Engineering dynamically decoupled quantum simulations with trapped ions

An external drive can improve the coherence of a quantum many-body system by averaging out noise sources. It can also be used to realize models that are inaccessible in the static limit, through Floquet Hamiltonian engineering. The full possibilities for combining these tools remain unexplored. We develop the requirements needed for a pulse sequence to decouple a quantum many-body system from an external field without altering the intended dynamics. Demonstrating this technique experimentally in an ion-trap platform, we show that it can provide a large improvement to coherence in real-world applications. Finally, we engineer an approximate quantum simulation of the Haldane-Shastry model, an exactly solvable paradigm for long-range interacting spins. Our results expand and unify the quantum simulation toolbox.


I. INTRODUCTION
The past few decades have seen a steady increase in the ability to control quantum states of simple systems, driven by improvements in maintaining the coherence of these systems. The most direct way to do this is to reduce their coupling to the environment through better engineering, which may involve materials developments or active stabilization of key parameters. However, since the development of NMR spin manipulation, a second strategy has also been known: improving the coherence of a system not in a fundamental way, but by strategically driving it with global operations that cancel out known and relatively slowly-varying noise sources [1][2][3]. This strategy, known as dynamical decoupling, has been extensively developed as a technique to reduce the coupling of quantum systems to external fields and thereby improve their internal coherences, which may allow them to serve as quantum memories [4][5][6]. Similarly, modern quantum computers often apply dynamical decoupling sequences to qubits that remain idle for part of a computation [7]. In both cases, one generally wants to completely cancel out any influence on the system. More recently, the dynamics of interacting quantum systems under a periodic external drive have themselves become a subject of extensive investigation. A primary example of this is the study of discrete time crystals [8][9][10][11][12][13][14]. Other recent works have explored Floquet Hamiltonian engineering: the ability to use a periodic driving sequence to design a desired effective Hamiltonian that the system does not natively realize [15][16][17][18][19].
Here, we combine these subjects to apply dynamical decoupling to a tunable quantum simulator. Although * current address: AWS Center for Quantum Computing, Pasadena, California 91125, USA. Work done prior to joining Amazon. † ksc64@umd.edu this technique is broadly applicable, we focus on its implementation with trapped ions, which are a leading quantum simulation platform [20,21]. This extends previous work in two ways. First, dynamical decoupling is applied not to preserve a quantum state, but to preserve evolution under a target Hamiltonian (or Floquet unitary operator). This adds the constraint, typically not present, that the decoupling pulse sequence not alter the desired dynamics. Demonstrating this on a trappedion simulator, we show that the coherence time can be substantially improved, in some cases by an order of magnitude. Second, while dynamical decoupling is often applied to a system with fixed intrinsic interactions, in a trapped-ion simulator there is significant freedom in the form of the Hamiltonian. The ability to periodically change the Hamiltonian in sync with global pulses opens up new possibilities for decoupling sequences, and for the combination of these with Hamiltonian engineering techniques. As a demonstration, we engineer an approximate decoupled quantum simulation of the Haldane-Shastry model, with dynamics reflecting key properties of this model including integrability [22,23]. This paper is organized as follows: Section II describes the general requirements for a decoupling pulse sequence to eliminate certain noise terms without affecting overall dynamics. Section III gives a brief overview of trappedion quantum simulation, with a particular focus on the primary sources of decoherence. With these preliminaries in hand, Section IV develops example dynamical decoupling and Hamiltonian engineering sequences for this trapped ion system and Section V presents tests with two ions validating the effectiveness of these sequences. Section VI extends this to larger systems, including an approximate realization of the Haldane-Shastry model incorporating dynamical decoupling and Floquet engineering. Finally, Section VII gives a summary and outlook.

II. PRINCIPLES OF DYNAMICAL DECOUPLING
We analyze dynamical decoupling sequences using the framework of Average Hamiltonian Theory, which concerns the average Hamiltonian H that approximately describes slow evolution under some pulse sequence [24]. For a periodic (Floquet) sequence, this is equivalent to the static Hamiltonian resulting from a series expansion (the Floquet-Magnus expansion) in the drive frequency [25]. Our goal will be to design a pulse sequence that makes H as similar as possible to the target Hamiltonian for the simulation, H t . Following previous expositions of Average Hamiltonian Theory, we consider the Hamiltonian in an interaction picture in which the operators act in a rotated frame determined by the previous pulses, which is often called the "toggling frame". For a pulse sequence consisting of evolution under Hamiltonians H 1 to H n for times t 1 to t n , which are each followed by pulses P 1 to P n (see Fig. 1), the unitary evolution operator ( = 1) is resulting in the toggling-frame Hamiltonians H n : A basic result of Average Hamiltonian Theory is that the resulting dynamics are approximately described by evolution according to the following static Hamiltonian: with T the total evolution time including the pulses. In words, the average Hamiltonian is, to lowest order in |H n t n |, nothing more than the weighted average of the different rotated Hamiltonians that make up the sequence. Within this framework, we define a dynamically decoupled sequence for some operator O, often an undesired noise term, as a pulse sequence which results in an average Hamiltonian that does not include O, despite O appearing in at least one of the individual H n of the sequence.
To design dynamical decoupling sequences that are suitable for a quantum simulation experiment, we usually impose two additional requirements on the pulse sequences. Most importantly, we require that, in the absence of noise (for O = 0), the unitary evolution is completely unchanged from that of H t , up to a known overall rescaling. To satisfy this, we also demand that the pulse sequence rotate the frame back to the original frame at the end of a cycle. These conditions are: P n · · · P 1 = I.
As a result, in the absence of noise, the Average Hamiltonian Theory description of the dynamics is exact. Higher-order corrections to H are determined by the commutators [H n , H m ], which all vanish. Consequently, provided that the pulse sequence is fast enough to effectively implement decoupling, the dynamics are guaranteed to follow the target Hamiltonian. While this is highly advantageous for many applications, in Sec. IV C we explore the possibility of removing this constraint to engineer average Hamiltonians that are not available in the static limit.
As a lesser constraint, we limit ourselves to decoupling sequences that are periodic in time. Accordingly, U 1 becomes the unitary evolution under one cycle, with the total evolution after m steps being U = Π m (U 1 ) m . Combined with the previous condition, this allows us to vary the total evolution time by integer cycles and observe time evolution matching the target Hamiltonian, as long as decoupling is effective.
While simple dynamical decoupling sequences can often be designed by inspection, more elaborate sequences can be systematized in several ways. Previous work [16] has developed an approach to determining H based on a pulse-sequence matrix representation, focusing on the case of a fixed H n . The Floquet-Magnus expansion can also be used to calculate higher-order terms in H [25] (see Appendix IX C), which can be suppressed at the cost of longer, symmetrized pulse sequences. Finally, we note that allowing non-periodic pulse sequences can be advantageous for cancelling certain types of noise [4], although such sequences may complicate the interpretation of dynamics with a varying evolution time.
We now provide a brief description of how the H n and P n operations are generated in our trapped-ion platform, and apply this framework to the specifics of this system.

III. QUANTUM SIMULATION WITH TRAPPED IONS
Here we give an overview of the generation of longrange Ising-like Hamiltonians using the Mølmer-Sørensen scheme [26], following previous treatments [20,21,27]. We consider a linear chain of N ions in a radiofrequency (rf) Paul trap [28], which is uniformly illuminated by multiple optical tones nearly resonant with a series of joint spin and motional transitions (often referred to as red and blue sideband transitions). Tone λ generates a Hamiltonian term of the following form: FIG. 1. Schematic of the dynamical decoupling scheme for quantum simulations. a) We consider a spin system that can be subjected to two classes of controllable operations (illustrated as driven by a pair of laser beams): overall global rotation pulses (Pn), and Hamiltonian evolution for some period of time (Hn). The Hamiltonian evolution induces interactions (represented by fluctuating spin angles), but can also introduce a coherent error (represented by a tilt of in the average spin axis). b) The goal of the dynamical decoupling sequence is to alternate periods of Hn and Pn operations, so that the final evolution after the pulse sequence corresponds to evolution (up to an overall rescaling) under a target Hamiltonian Ht while the error is cancelled out. c) Level schematic for spin-spin interactions generated by the Mølmer-Sørensen technique described in Sec. III. For simplicity, the two-ion case is shown and only transitions from the absolute ground state are drawn.
Here j is the ion index and ν is the normal mode index, a ν is the destruction operator of a phonon of motion for a normal mode of the ion chain, σ + j is a raising operator for the two-state pseudospin (or qubit), η is the Lamb-Dicke parameter, b ν,j is the mode ν amplitude for ion j, ω ν is the mode frequency, and Ω λ , µ λ , and φ λ are the Rabi rate, detuning from the qubit transition, and phase of tone λ, respectively. We have taken the Lamb-Dicke limit, implying that η 1. These tones can be used to generate a family of Hamiltonians, of which we briefly describe the most relevant cases.
First, with two tones that are symmetrically detuned near the motional transitions, µ 1 = −µ 2 = µ, the evolution can be analyzed using the Magnus expansion and shown over long times to approximately follow Hamiltonian evolution with the following form [21]: where φ s = (φ 1 + φ 2 + π)/2 (see Fig. 1c). The first term, describing an effective spin-spin interaction, can be approximately written (for µ > 0) as an Ising-like term with a power-law spatial distribution that is tuned by µ, while the second term contains the AC Stark shifts (or light shifts) of the qubit levels from the two frequencies that vanish when they have balanced strengths. Second, with one tone of µ = 0, this Hamiltonian (Eq. 8) may be simplified to an effective magnetic field term: Such a tone may be applied in isolation, to realize global spin rotations, or in combination with an interaction term.
In addition to these directly generated Hamiltonian terms, an effective B z field can be implemented by means of a rotating frame transformation that redefines the qubit frequency. In Appendix IX B we describe how to combine this continuous rotating-frame transformation with the toggling-frame from the pulses to use this technique in dynamical decoupling sequences.
A common target Hamiltonian for quantum simulations is a long-range transverse Ising model, such as: However, several errors can cause the observed dynamics to deviate from evolution according to this Hamiltonian. These errors include: • Fluctuations in the ratio of Ω 1 /Ω 2 generating the interactions, causing the AC Stark term to drift away from zero. This is of special concern because the AC Stark shifts, unlike the spin-spin couplings, are not suppressed by η 2 (see Eq. 9). As a result, they are often much larger than the intended terms in the Hamiltonian (typically about an order of magnitude larger than J j,j+1 for our experimental parameters) and therefore require a high degree of fractional stability.
• More generally, fluctuations in any of the laser intensities generating the Hamiltonian terms, causing the values of J j,j or B y to vary over time.
• Errors arising from finite population of the motional states. This can occur from direct excitation of these transitions by the tones, or due to heating processes that these tones couple into the spin dynamics (see Appendix IX E) [21].
In the following sections we demonstrate dynamical decoupling sequences that suppress the AC Stark noise, and show that this leads to a large increase in coherence for sufficiently large detuning µ.

IV. APPLICATION TO TRAPPED IONS AND EXAMPLE SEQUENCES
We now apply dynamical decoupling to an experimental trapped ion simulator, in which noise enters as a (possibly site-dependent) effective B z field: O = j j (t)σ z j . Several features particular to this experimental platform prove to be especially convenient for this goal. First, the noise varies relatively slowly compared to the duration of pulses. Second, both the pulse operations and Hamiltonian terms are generated by laser tones, and can be turned on and off or modified individually or in concert. Finally, the noise is only present when the interactions are on, and is approximately independent of the details of the Hamiltonian. Two consequences of this control are worth highlighting. First, there is essentially no unwanted evolution during the pulses. This can be contrasted with many natural systems with fixed interactions, in which there is some undesired evolution from the interactions during the finite lengths of the decoupling pulses. Second, while the average Hamiltonian that results from a pulse sequence applied to a fixed Hamiltonian is constrained by the symmetry group associated with the pulses [2], the ability in this system to turn on and off different Hamiltonians in sync with the pulses allows us to circumvent this limitation and realize average Hamiltonians with arbitrary symmetries. This is demonstrated explicitly in the following section.
It is useful to explicitly provide the transformations of the Pauli operators under global π pulses. These may be summarized as: where k, k = {x, y, z} index the Pauli operators (excluding the identity) and R k (θ) = e −i j θσ k j /2 . Geometrically, a π rotation about axis k flips the sign of a Pauli operator unless the axis of rotation and the Pauli operator direction coincide. This allows for a natural categorization of operators based on their parity under a given π rotation. For example, under the rotation R y (π), σ x j is odd but σ x j σ x j is even. Our strategy, similar to earlier proposals for high-fidelity quantum gates [29][30][31][32][33], is to use this difference to approximately invert the part of the unitary operator corresponding to evolution from the noise, undoing this evolution over two interaction periods separated by a π pulse, while the desired evolution is unchanged. Looking back at the fundamental Hamiltonian (Eq. 8), this is possible because although it generically has single spin operators with odd parities, the overall time dependence means that changing the sign is not equivalent to inverting the unitary. Instead, the Mølmer-Sørensen interaction relies on a geometric phase accumulated over the ions' spin-motion trajectory [34], which is invariant under a π rotation.
A. Example 1: CPMG sequence As a minimal example of a decoupling sequence satisfying our requirements, we consider a typical asymmetric Carr-Purcell-Meiboom-Gill (CPMG) decoupling sequence [3]. The target Hamiltonian includes spin-spin interactions and uniform fields along each direction: For this sequence n = 2 (of Eq. 2), t 1 = t 2 , and both pulses are π pulses along one axis taken to be y. A single unit of the evolution is: Note that we have explicitly included the noise term O = j (t)σ z j in both H 1 and H 2 . This sequence satisfies our requirements, with the terms in the target Hamiltonian remaining unchanged for different reasons. The spin-spin interaction term has even parity under any π pulse, and is therefore unchanged by R y (π), the B y term has even parity under R y (π) because their axes are aligned, and finally the B x and B z terms have odd parity under R y (π), but are correspondingly flipped between H 1 and H 2 to undo their transformation. As a result, the evolution under the target Hamiltonian is unchanged by the pulses, only occurring in a rotated frame that returns to the initial frame after both π pulses. However, the σ z j noise, which has odd parity and does not rotate with the pulse frame, reverses sign after each π pulse and is approximately cancelled. The exact degree of suppression depends on the power spectrum of (t) compared to the pulse frequency [4], as well as the commutator of the noise with the target Hamiltonian (see Appendix IX C).
As promised in the introduction to this section, this example demonstrates that the ability to vary the Hamiltonian during the pulse sequence has expanded the range of possibilities for H. If we were limited to a single Hamiltonian, for example H 1 , it would not be possible to have terms proportional to σ x j or σ z j in H, because they do not respect the symmetry of the pulse sequence. Thus, controllability of the Hamiltonian has significantly expanded the range of models amenable to dynamical decoupling.

B. Example 2: XY sequence
The second sequence we consider is a version of the XY dynamical decoupling sequence [3] for the same H t (Eq. 13): This sequence traverses a more complex rotating frame that only returns to itself after four global rotations. As a result, the signs of the B field terms must be cycled in the following way from H 1 to , while the spinspin term is unchanged for each Hamiltonian. While this has additional complexity relative to the CPMG sequence, it has the advantage of being insensitive to pulselength errors [3], due to a chirality condition satisfied by the pulse sequence [16], and also decouples noise along the x and y directions in addition to z. A similar robustness to pulse-length errors can be achieved in the CPMG sequence by alternating the sign of the rotations between R y (π) and R y (−π).

C. Example 3: Hamiltonian engineering a decoupled Heisenberg model
Until now, we have imposed the stringent requirement that the Hamiltonian be exactly the same in each rotated frame. This ensures that the dynamics are unchanged by the pulse sequence, removing any higher-order terms in the Floquet-Magnus expansion. However, one can also relax this condition to engineer a desired average Hamiltonian that is not accessible in the static limit. For example, in a toggling frame set by R y (π/2) pulses, σ x j σ x j is rotated to −σ z j σ z j , which is not otherwise accessible. Recently, this concept of Hamiltonian engineering has been used to extend the capability of quantum simulators [16][17][18][19]. While this imposes the extra condition that the cycle time (which can now be considered as a Trotter step [35,36]) be sufficiently small, within this limit it is natural to combine this technique with dynamical decoupling. For example, here is a pulse sequence that realizes, in the limit of a high-frequency drive, a dynamically decoupled long-range Heisenberg model (H t = j<j J j,j σ j · σ j /3) that is insensitive to pulse-length errors: where In Section VI B we apply this sequence to experimentally engineer the Haldane-Shastry model.

V. TWO-ION TESTS OF PULSE SEQUENCE PARAMETERS
In the following sections, we experimentally benchmark these pulse sequences. The experimental apparatus (see Appendix IX A) consists of 2-10 171 Yb + ions trapped in an rf Paul trap. We use two magnetic fieldinsensitive ground hyperfine states as the pseudospin |↑ z and |↓ z and perform coherent operations consisting of individual spin rotations, global spin rotations, and global long-range spin-spin interactions.
We initially optimize and benchmark the effectiveness of dynamical decoupling with minimal tests using two ions. The sequence is a simplified version of the XY sequence [37]: Starting with both ions in | ↓↓ z , the ideal state after application of H XX (t) has a simple form: where J 0 is the spin-spin coupling of the averaged Hamiltonian H. Compared to the general form of the XY sequence in the previous section, both H XX and the average ion z magnetization σ z = j σ z j /N used here are symmetric under sign changes in x and y, allowing us to use a fundamental time step with only two π pulses, rather than four, without changing any measurement outcomes. Besides representing a basic building block of our quantum simulations, this sequence is also highly sensitive to the dephasing noise that we are trying to decouple. Fig. 2 shows a typical example of this technique. The interspersed π pulses effectively cancel out the slow AC Stark shift noise that is a primary cause of decoherence, while evolution under the target Hamiltonian is unaffected up to an overall scaling factor. Using data of this type, we quantify the success of the decoupling with the scaled coherence time J 0 τ , the product of the coherence time τ with J 0 . More specifically, we fit the average magnetization σ z to the function g(t) = A(1 − e t/τ cos (2J 0 t)) − 1. AC Stark shift noise leads both to finite τ , due to dephasing, and to A < 1, due to a nonzero average magnitude of the energy difference between the states |↓↓ z and |↑↑ z .
For the data shown in Fig. 2, the fit values (fits not shown) for J 0 /(2π) are 0.52 kHz for bare evolution and 0.33 kHz with dynamical decoupling. The fit values for τ are 1.2 ms and 7.5 ms for bare evolution and with dynamical decoupling, respectively. This illustrates the central trade-off of these sequences: overall slower dynamics, but the potential for an increase in coherence time such that the product of the two improves (from J 0 τ = 3.9 to 15.6).
We now describe a few parameters of the dynamical decoupling sequence which must be chosen for optimal performance: the pulse shaping properties, the drive rate, and the detuning µ of the tones used to generate H XX .

A. Pulse shaping
We smoothly turn on and off both the interactions and the global rotations, to reduce unwanted motional transitions resulting from spectral broadening of the pulses. This is especially crucial for the interactions, which are relatively near-detuned from motional transitions. We use a Tukey window [10] with shaping parameter α = 2t p /T , where t p , the pulse shaping time, describes the amount of time that the pulse is either being ramped up or down compared to total pulse length T . For the interactions, the most relevant quantity is t p δ, the product of the pulse shaping time, which determines the degree of broadening, and the detuning from the nearest motional resonance δ = (µ − ω 1 )/2π. We experimentally find that setting t p δ ≥ 3 makes the induced error from motional excitations negligible, leading to a typical value of t p = 20 µs for our parameters. Global rotations, while less sensitive to pulse shape, are similarly given an α of 0.4.
The pulse shaping of the interactions leads to a decrease in the overall scale of the averaged Hamiltonian, which depends on the precise way that the rf waveform is imprinted onto the laser intensity. This can be theoretically predicted or determined with measurements of the ex-situ beam power, but we choose to characterize it with the in-situ ion dynamics. We find that the averaged spinspin coupling for a single pulse of the interactions obeys where J(t) ramps from 0 to J 0 . For our normal drive parameters, β is near 0.8, meaning that the pulse shaping removes 20% of the area under J(t). The pulse-averaged βJ 0 can then be used the standard formulas of Average Hamiltonian Theory. For example, in the XY sequence of Fig. 2, the average Hamiltonian for two ions is: For sufficiently fast drives, we see a decrease in J0τ , which is caused by a decreasing J0 without any compensating improvement in τ . b): Dependence of J0τ on detuning δ, for both decoupled (black) and regular (gray) sequences. Points are experimental data, while solid lines are numerics incorporating the dominant noise sources (see Appendix IX E). Decoupling is most effective in the regime of large detuning, δ/ηΩ > 5, in which the evolution without decoupling is nearly completely destroyed. For this data, the drive was fixed at t1 = 120 µs, and pulse shaping was fixed at tp = 20 µs throughout. Experimental data in both plots is averaged over 300 shots with error bars from the fit uncertainty.
The drive rate, or the value of 1/t 1 , is set by two considerations. The drive rate determines the bandwidth of noise that is suppressed, and must be chosen to be fast enough for the observed noise. The drive also sets the minimum time step, so it should be chosen so that the dynamics of interest are resolved. On the other hand, increasing the drive rate does not lead indefinitely to better performance because it slows down the experimental timescale, increasing sensitivity to any slow drifts that are not decoupled.
Varying our experimental drive rate over the range of 40-200 µs (Fig. 3), we observe nearly no decrease in J 0 τ with increasing t 1 , indicating that the noise we are cancelling out is mostly slower than the kHz scale. At sufficiently fast rates of driving, we see a decrease in J 0 τ because the fast drive leads to a reduction in the average Hamiltonian parameter J 0 (see Eq. 20) without a corresponding improvement in τ .

C. Detuning
The Mølmer-Sørensen detuning δ is a key parameter in ion trap quantum simulations that determines the power-law range of the spin-spin coupling. It also controls the relative strengths of various sources of error. When δ is relatively small, coupling to the intermediate motional states is strong, and therefore errors associated with this motional coupling, such as fluctuations in the motional resonance frequencies, determine the coherence time. When δ is large, Stark shift noise becomes dominant, with the two being roughly equal near our normal point of operation of δ/ηΩ ≈ 3. Fig. 3 shows the dependence of J 0 τ on δ, for sequences with and without dynamical decoupling. At small δ the two are similar, while at large δ the dynamical decoupling begins to perform better. At the optimum value of δ/ηΩ ≈ 5, dynamical decoupling improves J 0 τ by over an order of magnitude. The dependence is captured well by a numerical model incorporating the relevant error sources (see Appendix IX E and Fig. 6). This suggests that this technique is especially valuable for quantum simulations requiring shorter-range interactions, which are reached with large δ.

A. Dephasing test
It is desirable to check that these techniques extend to larger ion chains. This requires a change in the methodology for determining J 0 and τ , which were previously found by a fit to an analytical form for the two-ion dynamics. For the following multi-ion studies, we calibrate J 0 using the Mølmer-Sørensen formula (Eq. 9) and Average Hamiltonian Theory, while τ is estimated from an exponential fit to an observable of the system that is conserved by the ideal Hamiltonian but modified by decoherence. As a minimal demonstration, we have prepared a chain of eleven ions, whose spin-spin interactions obey an approximate long-range power law with a similar nearestneighbor spin-spin coupling as in the two-ion tests (see Appendix IX A). The spins are individually initialized (with a single-spin addressing beam [38]) in a product FIG. 4. Performance of dynamical decoupling to preserve trivial evolution under a many-body Hamiltonian. Ions are initialized as either up or down along x, and a long-range spin-spin Hamiltonian HXX is applied for variable time, with and without CPMG decoupling pulses, before the x magnetization is read out. In the absence of decoupling, stray AC Stark shifts induce dynamics and dephasing, while with decoupling these are suppressed. The parameters are δ/ηΩ = 7.5, tp = 20 µs, and t1 = 120 µs. Evolution time is scaled by the nearest-neighbor spin-spin coupling averaged over the pulse sequence, J0. Left panels: average magnetization of individual spins. Right panel: generalized imbalance I(t), measuring preservation of initial state (see text), for regular sequence (gray) and decoupled sequence (black). Points are the average of 500 experimental repetitions, with error bars smaller than the symbol size. state in which each spin has a definite value of σ x j , and evolved under a CPMG-type sequence: To avoid errors associated with the rotations, we replace each bare rotation R −y (π) with a composite BB1 pulse [39]. Ideally, this shows no dynamics, since the system is initialized in an eigenstate of H XX , making any effects of decoherence clearly identifiable. However, with the addition of a small B y field this system also exhibits domain-wall confinement [40], which can result in nontrivial, slow dynamics between the confined quasiparticles. As a result, we expect that this sequence provides a reasonable estimation of the coherence we could expect when studying an interesting many-body system. As shown in Fig. 4, dynamical decoupling again greatly helps to preserve this state. This can be quantified with the generalized imbalance I(t), a spin-spin correlator reflecting the memory of the initial spin configuration [41]: I is a constant under the target Hamiltonian, while dephasing or thermalizing dynamics will both drive it towards zero. Fitting I(t) to an exponential decay, the scaled coherence time without dynamical decoupling is J 0 τ = 4.2, while with dynamical decoupling it is extended to J 0 τ = 22.1, showing that the benefit of dynamical decoupling demonstrated in two-ion tests persists in a much larger system. Addition of the decoupling pulses significantly extends the timescales accessible in this simulation, potentially improving the ability to resolve slow dynamics.

B. Decoupled Floquet Hamiltonian Engineering:
Realizing the Haldane-Shastry model As a final demonstration of these tools, we apply Floquet Hamiltonian engineering to a long-range Heisenberg model in the Haldane-Shastry regime. The Haldane-Shastry model [22,23] can be formulated as a periodic spin-1/2 chain with long-range antiferromagnetic spinspin interactions: with sites that are evenly spaced around a circle. This model has attracted extensive interest because it is exactly solvable using the asymptotic Bethe ansatz, and features a spin-liquid ground state and fractionalized spinon quasiparticles [42]. Previous work proposed to experimentally create this model with trapped ions [43][44][45] or atoms in a photonic crystal waveguide [46]; however to our knowledge it has not been demonstrated. We approximately engineer the Haldane-Shastry model using the dynamical decoupling sequence for a Heisenberg model (Eq. 18 and Fig. 5a) with four ions whose spin-spin interactions obey an inverse square power law: J j,j ≈ J 0 /|j − j | 2 . Because the sequence now relies on a fast drive limit, we set t 1 J 0 = 0.05. This model is a convenient target for our technique because reaching inverse-square power law couplings requires a large detuning of δ/ηΩ = 9.9. At this large detuning, the scaled coherence time without dynamical decoupling is J 0 τ < 1, as suggested by Fig. 3b, making decoupling necessary to see any nontrivial dynamics.
Our experimental realization deviates from the ideal Haldane-Shastry model in several ways. At the level of the Hamiltonian, the experiment has open boundary conditions, unlike the original version of this model, although open generalizations have also been previously studied [47,48]. Furthermore, the spin-spin couplings in the ion chain do not precisely obey an inverse square power-law dependence [21]. For the data shown, the best-fit power law constant is 2.05, and the largest deviation of a coupling from the fit is 0.08 J 0 . Beyond the Hamiltonian form, the experiment also has decoherence, which in this regime is primarily due to motional heating (see Appendix IX E). Averaging over the pulse sequence, this may be approximated as depolarizing noise (see Appendix IX F for numerics incorporating decoherence). Despite these differences, we can experimentally and numerically observe key signs of proximity to the ideal Haldane-Shastry limit. Due to the Heisenberg symmetry, any finite initial magnetization along each direction is conserved up to a slow decay set by the decoherence (Fig. 5b). We estimate the coherence time by a fit to the average of these decays, giving J 0 τ ≈ 4.5. Meanwhile, preparing another initial state, such as a Nèel initial state along z, results in persisting and oscillatory nonthermalizing dynamics consistent with near-integrability (Fig. 5c). Similar oscillations and recurrences are seen in numeric simulations of the experiment, which are performed by solving the Schrödinger equation for H using the experimental spin-spin couplings (Eq. 9), and in simulations of the exact Haldane-Shastry dynamics (Fig. 5d). We note that the inversion of the first recur-rence, seen in experiment and experimental numerics but not the ideal Hamiltonian, is a boundary effect that is determined in an open chain by the parity of the number of ions. In additional data (see Appendix IX G and IX H), we show that this behavior persists for an alternate initial state, and numerically study the level statistics of the approximate and exact Hamiltonians. Furthermore, utilizing the flexibility of Floquet Hamiltonian Engineering, we modify the driving sequence to engineer a lesssymmetric thermalizing model, for which each of these signatures is lost: the experimental dynamics no longer conserve σ z and no longer exhibit revivals (see Fig. 8).
Approximate realization of the Haldane-Shastry model opens up a number of exciting possibilities. These include studies of prethermalization and relaxation in a nearintegrable system [49][50][51], and of the thermodynamic consequences of the expected fractionalized quasiparticles [52]. We leave a full theoretical and experimental study of the trapped-ion implementation of this model to future work.

VII. OUTLOOK
We have demonstrated a general strategy to extend dynamical decoupling to quantum simulators undergoing unitary many-body evolution. Applied to a trapped-ion experimental platform, we have shown that this technique may extend its coherence time and accessible parameter regime, and can be integrated with Floquet Hamiltonian engineering to further expand simulation possibilities. We expect that these decoupling strategies will be useful in systems with other noise sources as well. For example, our experimental qubit is insensitive to external magnetic fields, making decoherence from physical magnetic field fluctuations (as opposed to simulated magnetic fields from qubit energy shifts) negligible. However, these are a leading source of decoherence in other types of trapped ion experiments, and are amenable to the techniques developed here. These methods are also potentially relevant to other platforms with similar control, including quantum simulators based on neutral atoms [53][54][55] or superconducting circuits [56]. Extensions of this approach, in general including rotations on individual spins rather than global rotations [15,57], can be applied to decouple essentially any noise that can be modeled as arising from a slowly-varying Hamiltonian with undesired couplings, while preserving a target Hamiltonian of arbitrary symmetry properties.
Our results also help to clarify the limiting noise in various regimes of our simulator's operation, suggesting possible paths forward. With dynamical decoupling, our noise model (see Appendix IX E) suggests that the limiting noise over most regimes of interest is slow drifts in the motional trap frequencies, possibly due to temperature fluctuations that shift the resonant frequency of the rf trap voltages or small shifts in the electronic grounding. Therefore, future improvements should target this noise, which we have not investigated as extensively as other error sources. We can further predict that operating at δ/ηΩ = 3.5, if we can reduce these errors arising from motional frequency drifts by about a factor of 5 so that they are below other error sources, we will be able to increase our coherence to J 0 τ ≈ 30. In addition to hardware improvements, this will be facilitated by advances in automated and efficient re-calibration of experimental parameters, which are already in use in similar experiments [58]. Once achieved, assuming that an experimental signal typically persists up to time t = 2τ , this would correspond to a time of tJ 0 ≈ 60 available for quantum simulation.
While we have focused on demonstrations with small systems, the fundamental error sources we have studied couple to each ion equally, providing promising fundamentals for system size scaling. For sufficiently long ion chains, the decreasing axial frequency leads to an additional source of decoherence [59] that may ultimately require moving to a modular architecture [60]. However, even direct scaling has led to quantum simulations with up to about 50 ion chains [61,62]. Maintaining this noise up to comparable sizes would enable our simulator to reliably investigate dynamics that are highly challenging to access either for classical simulation or near-term gatebased quantum computers [63].

A. Experimental details
Our apparatus has been described in previous recent works [12,40,41,64]. For the data shown here, we create crystals of 2 to 11 171 Yb + ions in an rf Paul trap with anisotropic secular trapping frequencies of ω xy ≈ 2π× 4.8 MHz and ω z = 2π× 0.5 MHz. Pseudospin states are encoded in two ground-state hyperfine levels: |F = 0, m F = 0 = |↓ z and |F = 1, m F = 0 = |↑ z . State initialization via optical pumping prepares all the spins in | ↓ z with a fidelity greater than 0.99 per ion. For the data shown in Fig. 2, the average magnetization is read out using a photomultiplier tube (PMT), with a typical detection fidelity of 0.99 per ion. For all other data, individual spins are detected using fluorescent light imaged onto a charge-coupled device (CCD) sensor, with a typical detection fidelity of 0.97 per ion.
Coherent global operations are created using stimulated Raman transitions, driven between two optical frequency combs generated by a pulsed laser [65], with typical π times of t π = 5 µs. Interactions are generated with three frequency combs set up to drive two stimulated Raman transitions, creating the two tones for the Mølmer-Sørensen scheme. The resulting effective coupling parameters primarily depend on the laser power and detuning, while only weakly depending on system size. For the two-ion data presented, a typical spin-spin coupling rate is J 0 = 2π×400 Hz (δ/ηΩ ≈ 4). For the multi-ion data presented in Fig. 4, the typical energy scale is J 0 = 2π×204 Hz (δ/ηΩ = 7.5), and the spin-spin couplings approximately follow a power law dependence, J j,j ≈ J 0 /|j − j | p , with p = 1.32. For the four-ion Haldane Shastry data presented in Fig. 5, J 0 = 2π×84 Hz (δ/ηΩ = 9.9) and p = 2.05.
For the data in Sec. V, involving tests with two ions, we extract both τ and J 0 self-consistently from the oscillatory dynamics, as described in the main text. For the multi-ion data in Sec. VI and Appendix IX G, there is no longer a single oscillation frequency at 2J 0 . We therefore calculate J 0 from calibrations of the experimental parameters and Eq. 9.

B. Dynamical decoupling with rotating-frame B z terms
Our quantum simulator platform has the capability to implement an effective B z field through a rotatingframe transformation. To perform this transformation, the qubit transition is defined with an energy shift relative to the "true" energy splitting of the atomic states. Practically, this requires changing the frequencies of the Mølmer-Sørensen beams to ±µ−2B z , making them symmetrically detuned from the redefined transition, while also calculating the phase for all coherent operations, normally defined in the rotating frame of the qubit transitions, with a corresponding shift of −2B z t.
It is straightforward to combine this rotating frame with the toggling frame conditions for dynamical decoupling. As a concrete example, consider the CPMG sequence of Example 1, consisting of two periods of evolution with interactions separated by two R y (π) rotations. In the toggling frame of the rotations, the two periods of interactions alternate between B z and −B z . Accordingly, φ(t) should alternate between −2B z t and +2B z t, with boundary conditions that maintain phase continuity. Applying these conditions at each transition, we find that for the first part of each cycle, which begins at time t 0 + 2(n − 1)(t 1 + t π ), the phase should be φ(t) = −2B z t + 2B z (t 0 + 2(n − 1)(t 1 + t π )), while for the second half beginning at time t 0 + (2n − 1)(t 1 + t π ) it should be φ(t) = +2B z t − 2B z (t 0 + 2n(t 1 + t π )).

C. Error analysis for the CPMG sequence
Here we provide a more detailed error analysis for the CPMG sequence shown in Section IV A, focusing on the degree to which dynamically decoupling suppresses the noise. We first separate the Hamiltonian terms into ideal (static) terms and the error: The evolution operator for one cycle is: tc (H02+O(t))dt × R y (π)T e −i ta +t 1 ta (H01+O(t))dt (29) with T the time-ordering operator. The exact frame transformation allows this to be rewritten as: At this point, we invoke the Magnus approximation, in which the average Hamiltonian is calculated in progressive orders of the commutator between terms, which we explicitly write up to second order: (37) We can perform the same manipulations for the term beginning at t a : T e −i ta +t 1 ta Ht−O(t)dt = e −it1Ha , and combine exponentials with the Baker-Campbell-Hausdorff (BCH) theorem: Finally, the above expression reduces to: where we have defined variables related to the mean and fluctuations of O over a duration of t 1 : t 1 n = tn+t1 tn (t)dt and (t 1 ) 2 δ n = tn+t1 tn t tn ( (t) − (t ))dt dt. It can be seen from this expression that the lowestorder noise is caused by any drift over the experimental sequence, while higher-order contributions come from the commutator of the noise term with the target Hamiltonian. Both can be reduced by reducing t 1 , while additional symmetrization of the pulse sequence can be used to partially cancel higher-order terms [35].

D. Drive error
We benchmark our drive error by applying a given drive to a single ion initialized in | ↑ z 500 times (for 1000 total π pulses) and fitting the decay. For the XY sequence, the fit is consistent with no decay over this period, and the fit parameters naively imply a lower fidelity bound of F > 0.99997 per drive period (two π pulses). For the BB1 sequence, we measure the fidelity of a single driving period as F = 0.9998. As a typical experimental time scan uses 30 drive periods, the drive error is therefore negligible when using either robust sequence.

E. Noise model
For the noise model used in Fig. 3, we simulate the full Mølmer-Sørensen evolution (Eq. 8) for both the red and blue tones, truncating the maximum number of phonons in each motional mode to two. The time step for the numerical evolution (using a Krylov subspace algorithm) is fixed at 15 ns. In addition, we add three types of experimental noise: • Stark shift noise. This consists of a random, independent intensity fluctuation on each of the two tones. These are constructed to have a power spectrum of 1/f 2 , over the bandwidth of [100 Hz, 1 MHz], and to have an overall fractional variation of σ = 0.021 over the total sample time of 100 ms.
• Detuning noise. This consists of a Gaussian random variable, with σ = 4.5 kHz, this is applied as a static shift to the detuning δ for each experimental run.
• Motional heating. This consists of a displacement operator with a random phase that is applied to the highest motional mode (the center-of-mass mode) every 1.5 µs, representing fast fluctuations of the global electric fields holding the ions [66]. The amplitude of these random displacements is set at 0.01 to match the decoherence observed in Fig. 3. This results in a simulated heating rate of roughly 50 quanta/s. We can compare this to the rate observed in independent measurements, which is typi-cally between 100 and 200 quanta/s, pointing to future opportunities to further refine our noise models with independent characterizations of the various processes.
We note that all of these are technical noise sources. The fundamental limit to coherent Mølmer-Sørensen evolution, off-resonant Raman scattering, is several orders of magnitude below these for the chosen parameters. The strengths of the Stark shift noise and detuning noise were taken to match the data in Fig. 3, as were the exact power spectra over the experimental timeframe (choosing from possibilities of 1/f noise, 1/f 2 noise, and an ensemble average of constant values). The strengths are similar to previous estimations of our error sources [12,64], while the power spectra over experimentally relevant times are not strongly constrained by our previous measurements. Simulations were run for evolution times up to 0.5 ms, and averaged over 20 random realizations of all the error sources, before being fit to extract J 0 τ . For numeric simulations of the dynamically decoupled sequence, the form of the experimental pulse shaping of the interactions was included, while we took the simplification of approximating the π pulses as perfect and instantaneous. Fig. 6 shows the relative contribution of each of these noise sources to 1/J 0 τ , as a function of the detuning δ/ηΩ. Consistent with our expectations, the decoupled dynamics effectively suppress the Stark shift noise without similarly affecting the other sources.

F. Haldane-Shastry numerics with decoherence
To support the comparison of experimental data with numerics, in Fig. 7 we present a comparison of the experimental approximate Haldane-Shastry data (reproduced from Fig. 5) with numerics incorporating decoherence. From the noise model described in the previous section, it can be seen that in the regime of the Haldane-Shastry data, δ/ηΩ = 9.9, for a decoupled simulation the dominant noise source is expected to be motional heating. This noise source can be treated as a stochastic σ x term [64], but after considering the overall global rotations about the Bloch sphere in the Haldane-Shastry sequence the resulting decoherence is reasonably approximated as a simple depolarizing model. This leads the density matrix for the system to be approximated as ρ = (1 − p 1 (t))I + p 1 (t)|ψ(t) ψ(t)|, with I the diagonal fully mixed density matrix, |ψ(t) the desired evolution, and p 1 (t) = e −t/τ .
We use the decay of the polarized data to estimate the rate of this decoherence process. Fitting the decay of the polarized data to a single exponential, we extract a coherence time of J 0 τ = 4.5. After applying the resulting correction to the numerics, the resulting dynamics are shown in Fig. 7. Reassuringly, the numerics incorporating decoherence show close agreement with experiment.

G. Modified Haldane-Shastry sequence
In Fig. 8 we show additional data for the Haldane-Shastry Hamiltonian engineering sequence. An additional initial state shows dynamics consistent with exact numerics, including a revival of the magnetization at late times.
To contrast with these results, and show the flexibility of our Hamiltonian engineering scheme, we also present data for a modified sequence in which two pulses are removed. This retains the decoupling of the Haldane-Shastry sequence, but results in an target Hamiltonian that has an anisotropic XY form: H HSmod lacks the integrability of the Haldane-Shastry model, as well as the conservation laws for each spin projection. Correspondingly, we see dynamics that do not show simple oscillations with revivals or conserve total spin.

H. Energy level statistics of Haldane-Shastry Hamiltonian
As described in the main text, our experimental Hamiltonian, while approximately realizing the Haldane-Shastry Hamiltonian, differs from it in several ways. While the exact Haldane-Shastry Hamiltonian can be described as spins arranged on a circle with couplings scaling as 1/r 2 , our experimental realization has open boundary conditions and couplings that are close to, but not precisely, a 1/r 2 power-law dependence. In the main text, we have presented numerics for our specific experimental sequence to demonstrate that some of the key signatures of the Haldane-Shastry Hamiltonian, such as approximately integrable dynamics, are robust to these changes. However, this comparison is necessarily limited to the dynamics following certain initial states. To provide an alternate comparison reflecting the entire structure of the Hamiltonian, we show in Fig. 9 the distribution of energy level spacings for the different Hamiltonians. Level spacings are quantified with the ratio of adjacent energy level gaps, defined as The statistics of r are an established probe of the extent to which a Hamiltonian is chaotic, by capturing the structure of the eigenvalues. For example, a transition from r obeying a Wigner-Dyson distribution to a Poissonian distribution is commonly seen for models undergoing a generic thermalizing to localized phase transition [67].
In the current application, the non-generic Hamiltonians can have levels that are exactly degenerate to within numerical precision, so we add an infinitesimal term to all energy level differences that is reduced until the results converge. Thus, indeterminate values appear as r = 1. In Fig. 9 the distribution of r is shown for three different cases: • an exact realization of the Haldane-Shastry Hamiltonian (Eq. 23), • the approximate realization of the Haldane-Shastry Hamiltonian in our experimental ion chain. This consists of the average Hamiltonian H corresponding to the pulse sequence Eq. 18, which describes a Heisenberg Hamiltonian with approximate 1/r 2 couplings and open boundary conditions, • and H for the modified pulse sequence which breaks the Heisenberg symmetry (Eq. 41).
For each, we numerically diagonalize a system of twelve ions, rather than the experimental four, so that there are enough eigenvalues to resolve the level structure. While this does have the effect of reducing the impact of boundary terms relative to the experimental data, it is nonetheless useful for understanding of the degree to which the class of Hamiltonians that we can experimentally create capture Haldane-Shastry physics.
At a qualitative level, it is clear that the Haldane-Shastry Hamiltonian shows a non-generic structure with peaks at certain spacings and a high degree of degeneracy (r = 0 or 1). The experimental variant of this Hamiltonian shows a very similar structure, with shifted weights that reflect the boundary conditions and small integrability-breaking terms. Finally, the modified Hamiltonian is very different from the others, and has a more generic structure that is similar to other longrange spin Hamiltonians that are believed to be chaotic [41]. Thus, the level spacing structure is consistent with the dynamics in suggesting that key properties of the Haldane-Shastry model are resilient to the changes in our experimental approximation.