Robust Dynamic Hamiltonian Engineering of Many-Body Spin Systems

We introduce a new framework for the robust control of quantum dynamics of strongly interacting many-body systems. Our approach involves the design of periodic global control pulse sequences to engineer desired target Hamiltonians that are robust against disorder, unwanted interactions and pulse imperfections. It utilizes a matrix representation of the Hamiltonian engineering protocol based on time-domain transformations of the Pauli spin operator along the quantization axis. This representation allows us to derive a concise set of algebraic conditions on the sequence matrix to engineer robust target Hamiltonians, enabling the simple yet systematic design of pulse sequences. We show that this approach provides a streamlined way to (i) treat any secular many-body Hamiltonian and engineer it into a desired form, (ii) target dominant disorder and interaction characteristics of a given system, (iii) achieve robustness against imperfections, and (iv) provide optimal sequence length within given constraints. Using this systematic approach, we develop novel sets of pulse sequences for the protection of quantum coherence, optimal quantum sensing and quantum simulation. Finally, we experimentally demonstrate the robust operation of these sequences in a dense ensemble of nitrogen-vacancy centers in diamond.

The key tool to engineer the dynamics of periodically driven systems is average Hamiltonian theory (AHT). This technique has been particularly successful in nuclear magnetic resonance (NMR), where periodic driving protocols enable the suppression of unwanted evolution due to both disorder and interactions, effectively preserving quantum coherence and enabling high-resolution NMR spectroscopy and magnetic resonance imaging (MRI) [25-27, 40, 61-77].
However, conventional control pulse sequences are generally optimized for solid-state nuclear spin systems where dipolar interactions dominate. In particular, these sequences are often not applicable to other quantum systems, such as electronic spin ensembles or arrays of coupled qubits, where either on-site disorder dominates or interactions have a more general form [78,79].
Furthermore, periodic driving schemes are often vulnerable to perturbations caused by inhomogeneities of individual spins in the system, non-ideal finite pulse duration effects, as well as imperfect spin state manipulation. While there exist * These authors contributed equally to this work † lukin@physics.harvard.edu many pulse sequences that retain robustness to some of these control imperfections [25,26,80], a systematic framework to treat these errors in a general setting of interest is still lacking, hindering the customized design of pulse sequences optimally adapted for various applications across different experimental platforms.
Our approach is based on a simple matrix representation of pulse sequences that allows for their analysis and design in a straightforward fashion, using intuitive algebraic conditions. This matrix describes the interaction-picture transformations of the S z operator, the Pauli spin operator along the quantization axis, which can also be visualized in a very intuitive way. Crucially, we show that by decomposing all pulses into π 2-pulse building blocks, this representation not only gives the effective leading-order average Hamiltonian describing the driven spin dynamics, but also provides a concise description of dominant imperfections arising from non-ideal, finiteduration pulses and rotation angle errors. More specifically, we show that (i) the suppression or tuning of on-site potential disorder, (ii) the decoupling or engineering of spin-spin interactions, and (iii) the robustness of the pulse sequence against control imperfections, can all be extracted directly from our representation and algebraic conditions. The algebraic conditions also analytically provide the minimum number of pulses required to realize a target application, thus ensuring minimal sequence length under given constraints. This approach thus allows the incorporation of Hamiltonian engineering requirements in the presence of imperfections, enabling the versatile construction of sequences designed for a particular quantum application and tailored to the detailed properties of the exper-imental system at hand (see Fig. 1).
Specifically, we use our formalism to protect quantum information and benchmark the performance of two sequences with different design considerations, each suited for systems in different regimes of competing disorder and interaction energy scales. We also utilize our framework to design pulse sequences for robust and optimal quantum sensing, where our method provides a generalized picture of AC field sensing protocols in which an external AC field in the lab frame translates into an effective vector DC field in the driven spin frame. We show how to achieve highest sensitivity by engineering the orientation and magnitude of the effective DC sensing field and choosing optimal initialization and readout directions while extending coherence times through disorder and interaction suppression. We then further apply our framework to quantum simulation and engineer the bare system Hamiltonian to a desired target form, providing a new avenue to study many-body dynamics over a wide range of tunable parameters with different types of interactions and disorder. Finally, we experimentally demonstrate our results in a disordered, dipolar-interacting nitrogen-vacancy (NV) center ensemble in diamond.
The main advances enabled by our approach include: 1. Robustness: We show that all types of average Hamiltonian effects, including errors resulting from pulse imperfections [Sec. III], can be readily incorporated as concise algebraic conditions on the transformation properties of the S z Pauli spin operator in the interaction picture. This leads to a simple recipe for sequence robustness by design.
2. Generality: Our approach is applicable to generic twolevel spin ensembles in a strong quantizing field, as typically found in most experimental quantum manybody platforms such as solid-state electronic and nuclear spin ensembles, trapped ions, molecules, neutral atoms, or superconducting qubits. Our framework covers on-site disorder, various two-body interaction types such as Ising interactions and spin-exchange interactions, as well as complex three-body interactions.
3. Flexibility: The flexibility of our approach allows Hamiltonian engineering that takes the energy hierarchy into account, which can be tailored to specific physical systems exhibiting different relative strengths between disorder, interactions, and control errors. This enables the development of pulse sequences designed for disorder-dominated systems, beyond the typical NMR setting. 4. Efficiency: Using simple algebraic conditions, we can find the shortest possible sequence length required to achieve a target Hamiltonian [Sec. IV.2]. In addition, we use combinatorial analysis to demonstrate the necessity of composite pulse structures for efficient sensing, and provide optimized sequences that achieve maximum sensitivity to external signals [Sec. V.3].
The paper is organized as follows: In Secs. II and III, we provide the theoretical framework for systematic pulse design, in Secs. IV, V and VI, we present system-targeted sequence design for the applications of dynamical decoupling, quantum sensing and quantum simulation, respectively. Finally, Sec. VII presents the experimental demonstration of our results to dynamical decoupling and quantum sensing. We conclude with a further discussion and outlook of the framework in Sec. VIII.

II. GENERAL FORMALISM AND FRAME REPRESENTATION
We start by introducing a simple representation of pulse sequences based on the rotations of the spin frame in the inter-action picture: Instead of illustrating a sequence by the applied spin-rotation pulses, we describe it by specifying how the S z spin operator is rotated by the applied pulses in the interaction picture (also known as the toggling-frame picture). This method provides a one-to-one correspondence with the average Hamiltonian of the system and is an extension of the method presented in Ref. [94], also closely related to control matrices [95,96] and vector modulation functions [97][98][99]. We efficiently depict this S z operator evolution using a simple matrix, and show that this direct link to the average Hamiltonian is valid for any system under a strong quantizing field. In addition to its simplicity in describing the decoupling performance for the case of ideal, instantaneous pulses, this representation also allows the formulation of concise criteria to treat pulse imperfections, as will be discussed in Sec. III.

II.1. Frame Representation
The dynamics of periodically driven systems can be described and analyzed using AHT [1] (See Appx. A for a detailed review of AHT). In particular, for a pulse sequence consisting of n spin-rotation pulses {P 1,⋯,n }, the leading-order average Hamiltonian, H avg , is a simple weighted average of the toggling-frame Hamiltonians where τ k is the pulse spacing between the k − 1th and k-th control pulses P k−1 and P k ,H k = (P k−1 ⋯P 1 ) † H s (P k−1 ⋯P 1 ) is the toggling-frame Hamiltonian that governs spin dynamics during the k-th evolution period, τ k , in the interaction picture, and H s is the internal system Hamiltonian.
Here, we present a convenient alternative method to calculate the leading-order average Hamiltonian, utilizing our toggling-frame sequence representation [94]. Our representation is based on the time-domain transformations of a singlebody S z spin operator in the interaction picture. As shown in Sec. S1A [100], the representation works for general twolevel system Hamiltonians under the rotating wave approximation in a strong quantizing field (secular approximation). Physically, this corresponds to the common situation, realized in almost all experimental platforms, in which energetic considerations require all interaction terms to conserve the total magnetization along the quantization axisẑ, which can also be written as [H s , S z tot ] = 0, where S z tot = ∑ i S z i is the total spin projection operator along theẑ-axis. Thus, our framework is widely applicable to different experimental systems, including both ordered and disordered systems, and systems with different types of interactions, including Ising [52,93], spin-exchange [78,79], dipolar [87], and even exotic threebody interactions [101][102][103].
To introduce our framework in detail, let us first consider two-body interaction Hamiltonians with on-site disorder. The most general form of such a Hamiltonian is (Sec. S1A [100]) where the first term H dis is the on-site disorder Hamiltonian and the second term H int is a generic two-body interaction Hamiltonian, h i is a random on-site disorder strength, {S x i , S y i , S z i } are spin-1/2 operators, and J I ij , J S ij , J A ij are arbitrary interaction strengths for the Ising interaction and the symmetric and anti-symmetric spin-exchange interactions, respectively. According to AHT, a pulse sequence periodically applied to the system can engineer this into a new Hamiltonian, dictated by the control field that dynamically manipulates the spins (see Appx. A).
In our framework, the control field is assumed to be a time-periodic sequence of short pulses, with each pulse constructed out of π 2-rotation building blocks around thê x,ŷ axes [ Fig. 1(c)]. In this setting, let us consider the interaction-picture transformations of the S z operator: is the global unitary spin rotation due to the control field. We will assume in this section that the pulses are perfect and infinitely short. In such a case, the +S z operator transforms into ±S x,y,z operators, depending on the rotation angles and axes of the pulses. Hence, the effect of the pulse sequence is a rotation of S z in time in a discrete fashion, and the transformation trajectory in the toggling frame can be identified as Here, P k is the global spin rotation performed right after the k-th toggling frame, t k = ∑ k j=1 τ j with t 0 = 0, and is a 3 × n matrix containing elements 0 and ±1. The matrix elements F µ,k can be explicitly calculated as Intuitively, a nonzero element F µ,k indicates that the initial S z operator transforms into S µ for the duration of the free evolution interval τ k , with its sign determined by F µ,k . Additionally, the time duration of each toggling frame is specified by the frame-duration vector This representation is illustrated for three pulse sequences in Fig. 2. The CPMG sequence [9][10][11], consisting of equidistant π pulses to suppress on-site disorder, can be represented as  since the first +S z is flipped to −S z by a π pulse [ Fig. 2(a,d)]. Similarly, the WAHUHA sequence, consisting of four π 2 pulses [87], is shown in Fig. 2(b,e). The matrix clearly shows how the spin operator rotates over time, cycling through all three axes to cancel dipole-dipole interactions [104]. Finally, we present a sequence that combines the ideas of WAHUHA and CPMG to echo out disorder while symmetrizing interactions, as depicted in Fig. 2(c,f). The representation thus far uniquely specifies the toggling-frameS z orientation after each instantaneous pulse. However, the rotation axis of π pulses is not yet uniquely specified. To address this, we decompose all pulses into π 2 building blocks, and specify intermediate frames for pulses of rotation angles larger than π 2. Such a π 2-pulse decomposition also simplifies the analysis of finite pulse duration effects, which we discuss in Sec. III. As shown in Fig. 3(a), a π pulse is then split into two π 2 pulses with zero time separation, where the first π 2 pulse along thex-axis rotates +S z into the intermediate frame +S y and the second π 2 pulse along the same axis brings +S y into −S z . For example, the CPMG sequence involving π pulses can be represented as which now unambiguously specifies the rotation axes for the π pulses. Note that the sequence contains zeros in its frameduration vector, serving to indicate the intermediate frames, and in the following we shall indicate them with narrow lines in the pictorial representation (see Fig. 3). The use of such intermediate frames also allows a natural description of composite π 2 pulse structures, which will play an important role in robust quantum sensing sequences [Sec. V].
One key advantage of our pulse sequence representation {F, τ } is that we can now conveniently obtain the engineered H avg [Eq. (1)] from any many-body Hamiltonian of the form Eq. (2). More specifically, the weighted row-sums and weighted absolute row-sums (which are equivalent to row square sums, since each element F µ,k takes on values {0, ±1}) of the sequence matrix, fully specify the generalized formula for the average Hamil- Here, the three terms in H int avg [Eqs. (8)(9)(10)] correspond to the transformations of the Ising interaction and the symmetric and anti-symmetric exchange interactions, denoted as H I avg , H S avg , and H A avg , respectively (see Sec. S1B [100] for more details).
For equidistant pulse sequences where τ = τ I 1×n , the above condition can be further simplified to ∑ n k=1 F µ,k = 0. This suggests that each row µ of the matrix F should have an equal number of positive and negative elements, such that their sum is 0, resulting in H dis avg = H A avg = 0. Physically, this corresponds to guaranteeing a spin-echo-type structure, in which each precession period around a positive axis is compensated by an equal precession period in the opposite direction. Applying this criterion to the sequences in Fig. 2, we see that as expected, the CPMG [ Fig. 2(a,d)] and echo+WAHUHA [ Fig. 2(c,f)] sequences cancel on-site disorder. The WAHUHA sequence, however, produces a residual on-site disorder term, also known as the chemical shift, given by H dis Fig. 2(b,e)]. In contrast, the terms with quadratic dependence on F µ,k cannot be fully suppressed in general, as the isotropic (Heisenberg) component of the interaction is invariant under global rotations [83], leading to H int avg ≠ 0. However, it is still possible to fully symmetrize these interactions into a Heisenberg Hamiltonian, H int avg = 1 3 ∑ ij (J I ij + 2J S ij ) ⃗ S i ⋅ ⃗ S j , which is in fact sufficient to preserve spin coherence in many situations; In particular, globally polarized initial states that are typically prepared in experiments constitute an eigenstate of the Heisenberg interaction, and consequently do not dephase under the Heisenberg Hamiltonian. Such interaction symmetrization is satisfied when L x = L y = L z [Eqs. (8)(9)], giving Again, for equidistant pulses, this condition is simplified to the statement that the sum ∑ n k=1 F µ,k should be the same for each µ. Based on this analysis, we verify that the CPMG sequence [ Fig. 2(a,d)] does not symmetrize interactions, since it only employs π pulses, while the two sequences that incorporate π 2 pulses to switch between all axes in the toggling frame [ Fig. 2(b,c,e,f)] do indeed symmetrize the interaction Hamiltonian into the Heisenberg form. The spin-1/2 dipolar interaction with J I ij = −2J S ij is a special case where the WAHUHA sequence fully cancels interactions to leading order, giving H int avg = 0. While we have focused on single-body and two-body interactions, the above analysis can be extended to interactions involving more spins. In particular, in Appx. B, we utilize results from unitary t-designs [105][106][107] to prove that the conditions described above also guarantee decoupling of general three-body interactions for polarized initial states.

III. ROBUST PULSE SEQUENCE DESIGN
For pulses of finite duration, on-site disorder and interaction effects acting during the pulses cause additional dynamics in the quantum system [25,26,80]. In addition, the spin rotations can also suffer from experimental control errors, such as over-or under-rotations. Both of these imperfections can contribute to an error Hamiltonian δH avg , which can be estimated to leading order using AHT as where t p is the duration of a π 2 pulse andH (0) P k is the zerothorder average Hamiltonian acting during the k-th π 2 pulse building block. Thus, the total leading-order effective Hamiltonian describing the driven spin dynamics is given by The goal of robust Hamiltonian engineering is to suppress the error δH avg = 0 by designing leading-order fault-tolerant, selfcorrecting pulse sequences.

III.1. Average Hamiltonian Analysis for Finite Pulse Duration
Turning to analyze finite pulse duration effects, we now provide an efficient method to understand and correct all pulse-related control errors. Here, the key insight is that our matrix representation F directly provides a simple way to obtain δH avg [Eq. (13)]. Intuitively, the form ofH (0) P k in δH avg is Table 1. Summary of robust dynamical decoupling conditions. A periodic pulse sequence consisting of n π 2 pulse building blocks is represented by a 1-by-n frame duration vector τ = [τ k ] and a 3-by-n frame matrix F = [F µ,k ] = [ ⃗ Fx; ⃗ Fy; ⃗ Fz], corresponding to the toggling-frame S z operator in different evolution blocks [Sec. II.1]. Based on this representation, the conditions for dynamical decoupling and fault-tolerance against leading-order imperfections can be phrased as intuitive statements on F, with precise algebraic conditions as listed. Above, tp is the duration of a π 2 pulse,êµ is the unit vector along axis µ, and µ, ν = x, y, z. For more details, see Sec. II.2, III.1 for conditions 1, 2, Sec. III.1 for condition 3 and definition of the "parity" of frame changes, Sec. III.2 for condition 4 and definition of the "chirality" of frame changes. expected to be the average of the neighboring toggling-frame Hamiltonians,H k andH k+1 , since the finite-duration pulse P k smoothly changes the spin frame from the k-th to the k +1th toggling frame [ Fig. 3(b)]. However, as shown below, detailed calculations reveal nontrivial prefactors, as well as an additional interaction cross-term that can be expressed as a parity condition on neighboring matrix columns.
Since our pulse sequences are constructed out of π 2-pulse building blocks, we can analytically calculateH (0) P k originating from the finite-duration pulse P k as where {P 1 , ⋯, P k−1 } denotes the preceding k−1 rotations and the effect of the pulse P k is given bỹ Here, Q k (t) = exp [−i ∑ i Ω i tS ν i ] is the time-dependent unitary operator due to the pulse P k that globally rotate spins along the ν-axis over the finite duration 0 ≤ t ≤ t p and Ω i is the Rabi frequency for a spin at site i, producing the π 2 rotation. For now, we assume no rotation angle errors: the treatment of them will be discussed in Sec. III.2.
Physically, the role of π 2 pulses is to smoothly interpolate the toggling-frame spin operatorS z (θ) = cos θS z k + sin θS z k+1 during the finite pulse duration, where θ = Ωt. Using this to evaluate the integral of Eq. (15), as detailed in Sec. S1C [100], we obtain whereH dis k ,H I k ,H S k andH A k are the disorder, Ising, symmetric and anti-symmetric spin-exchange interaction Hamiltonians at the k-th toggling frame, respectively, which can be obtained from replacing the original spin operators S µ in H s [Eq. (2)] to the toggling-frame onesS µ k , as shown in the individual terms in the summations of Eqs. (7)- (10). While most terms in Eq. (17) are the weighted average of the neighboring toggling-frame Hamiltonians, consistent with the original intuition, there is an additional average Hamiltonian contribu-tionH C k,k+1 resulting from two-body interaction cross-terms S z kS z k+1 acting during the pulse and given bỹ is the cross-interaction operator with J C ij = 1 π (J I ij − J S ij ) (see Sec. S1C) and P µν k defines the "parity" of neighboring k-th and k + 1-th frames, given as To cancel them, the sequence matrix should contain an equal number of even and odd parities for every pair of neighboring frames between two axes; if the spin frame maintains (changes) its sign when switching to a different axis, the parity is defined as even (odd). (f) Spin-rotation angle error effects. To cancel them, the matrix should contain an equal number of positive and negative chiralities for every axis; the rotation axis and sign of the spin frame rotation in the toggling frame define positive/negative chirality along that axis.
To cancel the interaction cross-terms, the parity should vanish when summed over one Floquet period for each pair (µ, ν): ∑ n k=1 P µν k = 0 [ Fig. 3(e)]. Intuitively, the parity can be understood as checking whether the signs of neighboring frames are the same (P µν k = +1, even parity) or different (P µν k = −1, odd parity).
Taking into account the distinct weighting factors of 4 π and 1 for the different interaction types as well as the additional interaction cross-term, as identified in Eq. (17), the effective Hamiltonian in the presence of finite pulse duration t p , H eff = H avg + δH avg , becomes where the base pulse sequence length, T = (∑ n k=1 τ k ) + nt p , includes the total length of n π 2 pulses. As described above in the discussion following Eq. (17), each of the terms in the toggling-frame Hamiltonian,H dis ,H A ,H I ,H S , andH C , can be readily computed using our sequence representation.

III.2. Analysis of Rotation Angle Error
We now analyze the effects of rotation angle errors in control pulses, resulting from imperfect and inhomogeneous global spin manipulation. At first glance, this may seem challenging, since the average Hamiltonian corresponding to a rotation angle error around a given axis in the lab frame depends on the transformations by previous pulses. However, we find that there is a simple intuition for their average Hamiltonian contribution in the toggling frame, whereby an imperfect rotation around the +μ direction (positive chirality) in the toggling frame can be compensated with another rotation around the −μ direction (negative chirality). Moreover, the rotation axis in the toggling frame can be readily described using our frame matrix F, allowing concise conditions based on rotation chirality to achieve self-correction of rotation angle errors in a pulse sequence.
More specifically, the rotation axis in the k-th toggling frame, ⃗ β k = ∑ µ β µ,kêµ , can be obtained by taking the cross product of the frame vectors before and after the pulse where ⃗ F k = ∑ µ F µ,kêµ . Physically, this cross product structure can be thought of as characterizing the chirality of the toggling-frame rotation from ⃗ F k to ⃗ F k+1 . As derived in more detail in Sec. S1D [100], the average Hamiltonian contribution corresponding to the rotation angle error is then given by where i is the static rotation-angle deviation from the target π 2 angle for a spin at site i. This allows us to identify the cancellation condition for rotation angle errors, δH rot avg = 0, as k β µ,k = 0 for each µ = x, y, z which corresponds to condition 4 in Tab. 1.

III.3. Decoupling Conditions for Finite Duration Pulses
By incorporating the dominant effects arising from finite pulse durations, we have obtained a more realistic form of the effective Hamiltonian. In the leading-order average Hamiltonian, we found that the finite pulse duration simply introduces corrections to the effective free evolution intervals in Eqs. (11,12), plus additional terms that are well-described by the parity and chirality associated with each toggling frame change. Combining these, we arrive at the conditions to achieve robust dynamical decoupling for a polarized initial state, as summarized in Tab. 1.
Remarkably, we have found that our matrix representation F provides a systematic treatment of all such imperfections in a simple, pictorial fashion. As an example, for the common case of equidistant pulses, condition 1 in Tab. 1 is satisfied when there is an equal number of yellow (F µ,k = +1) and green (F µ,k = −1) squares (lines) in each row, while condition 2 is satisfied when different rows have an equal number of squares, and an equal number of lines. Moreover, in Appx. B, we further show that our representations and methods can be applicable to more complex three-body interactions, even for finite pulse durations, highlighting the broad applicability of our sequence design framework.
Additional effects such as waveform transients and pulse shape imperfections can also be analyzed by a similar approach [80] (see for example Sec. S1E [100] for the treatment of rotation axis imperfections).

IV.1. System-Targeted Dynamical Decoupling
The goal of dynamical decoupling is to extend the coherence time by cancelling the effects of disorder and interactions, and by suppressing pulse imperfections in a robust fashion. In the field of NMR, a number of decoupling methods have been developed, such as multiple-pulse sequences [8-22, 24-26, 28, 61-73], magic echoes [40,74], frequency-and phase-modulated continuous driving [75,76], and numerically-optimized control schemes [27,77].
However, these sequences are optimized for interactiondominated dipolar-interacting spin systems only, making it difficult to extend them to other Hamiltonians exhibiting different energy scales and interaction forms. For example, electronic spin ensembles typically display strong on-site disorder with weak interactions [78,108], such that a naive application of the NMR pulse sequences performs poorly (see Sec. VII). In the following, we show that our framework allows the design of system-targeted dynamical decoupling sequences that tackle the dominant effects on a faster timescale to achieve better performance.
As an example, for disorder-dominated systems, disorder cancellation needs to be prioritized and performed on a shorter timescale compared to interaction symmetrization and control error suppression. Our representation directly reveals the individual decoupling timescales T W and T J for disorder FIG. 4. System-targeted dynamical decoupling. (a,b) Sequence details and their matrix representations. The inset illustrates the pulse legend used to denote π 2 and π pulses. Decoupling sequences can be characterized by their different timescales of disorder and interaction decoupling, denoted as T W and T J , respectively. In (a), Seq. A (Cory-48) is designed for interaction-dominated systems, where axis permutation is performed on the fastest timescale to prioritize interaction symmetrization over disorder cancellation (T W ≫ T J ). In (b), Seq. B is designed for disorder-dominated systems where echolike operations are performed on the fastest timescale to prioritize disorder suppression over interaction symmetrization (T W ≪ T J ). (c) Numerical simulation results of the dynamical decoupling performance of Seq. A and Seq. B for a wide range of disorder W and interaction J strengths. The figure of merit for the comparison is the coherence time T2, extracted as the 1 e decay time of initially polarized spin states. More specifically, we perform an exact diagonalization simulation of a disordered, interacting ensemble of 12 spins. In the simulations, we initialize all spins along thex-axis, periodically drive them with the pulse sequences, and measure the spin decoherence profile at stroboscopic times t = N T . Each parameter set is averaged over 100 disorder realizations and the averaged profile is used to identify T2. The π pulse duration and spacing are chosen to be tp = 25 ns and τ = 20 ns, respectively. We have verified that choosing different tp values does not qualitatively change the results. and interactions, respectively (illustrated in Fig. 4(a)): They can be quantified as the minimum length of toggling-frame time evolution that fulfills their respective decoupling conditions [Condition 1,2]. With W and J describing the characteristic disorder and interaction scales of the driven system, in the disorder-dominated case (W ≫ J), we thus require T W ≪ T J to reduce the magnitude of higher-order contributions to H eff [25,26,109] and maximize the leading-order approximation accuracy. If the control error magnitude is comparable to W, J, then chirality cancellation [Condition 4] associated with pulse imperfections should also be performed at a relatively fast rate to suppress higher-order errors.
To illustrate the importance of system-targeted design, we provide two periodic pulse sequences (Seq. A, B) both designed to robustly decouple disorder and interactions, but with Seq. A(B) better suited for systems characterized by stronger interactions(disorder). For Seq. A, we adopt the Cory-48 sequence [26] developed for nuclear spin systems, where spinspin interactions dominate over disorder. Indeed, as shown in Fig. 4(a), Seq. A symmetrizes interactions very rapidly and also cancels on-site disorder, but on a much slower timescale (T J ≪ T W ). It is also robust against leading-order imperfections resulting from finite pulse durations, and suppresses certain higher-order effects [26]. For comparison, we design a new sequence, Seq. B in Fig. 4(b), based on the conditions in Tab. 1, to make the sequence operate better in the opposite, disorder-dominated regime. Specifically, Seq. B incorporates frequent π pulses to echo out disorder on a rapid timescale while symmetrizing interactions on a slower timescale (T W ≪ T J ). We emphasize that Seq. B incorporates both π pulses and composite π 2 pulses when switching between toggling frames [ Fig. 4(b)], to accomplish fast spin-echo operations and retain robustness to control imperfections, and thus lies beyond the design capabilities of previous approaches.
Given these design considerations, we expect Seq. A to perform better in the regime of large interaction strengths (e.g. for NMR), and Seq. B to perform better for disorderdominated systems (e.g. for electronic spin ensembles). From numerical simulations, as shown in Fig. 4(c), we indeed see a crossover in performance as the disorder and interaction strengths are tuned in the system within the range 0.009 < W (τ + t p ) < 0.09 and 0.009 < J(τ + t p ) < 0.023. At small disorder values, Seq. A shows a longer coherence time T 2 than Seq. B. However, as we increase the disorder strength, we observe a crossover beyond which Seq. B outperforms Seq. A.
Overall, Seq. B shows a stable performance within the range of parameters studied, while Seq. A shows a strong susceptibility to disorder. This example illustrates how our formalism enables the systematic design of pulse sequences adapted to the dominant energy scales of different systems.

IV.2. Shortest Sequence for Robust Dynamical Decoupling
The algebraic conditions introduced in Tab. 1 greatly simplify the design procedure, thereby allowing one to not only design pulse sequences that are robust against certain imperfections, but also guarantee via analytical arguments the shortest sequence length to achieve a set of target Hamiltonian engineering requirements.
The conditions to cancel disorder [Condition 1] and symmetrize interactions [Condition 2] at the average Hamiltonian level require an equal number of ±1 frames along each axis, resulting in at least 6 distinct free evolution intervals when neglecting pulse imperfections, which implies that the echo+WAHUHA sequence in Fig. 2(c,f) has shortest length for ideal pulses. When incorporating finite pulse durations and rotation angle errors, the conclusion may be less obvious; however, using the algebraic conditions in Tab. 1, we find that the following pulse sequence, consisting of 6 free evolution periods of duration τ connected by composite π 2 pulses, satisfies all leading-order decoupling requirements: To the best of our knowledge, this is the first pulse sequence that decouples all leading-order imperfections and achieves pure Heisenberg interactions with only 6 free evolution intervals, illustrating the power of our formalism. We note that the minimum achievable length of the pulse sequence may be modified by experimental considerations: for example, it may be challenging to apply pulses with different phases in close succession due to finite transient times of the experimental apparatus, and in such cases we can show that the minimum pulse length increases to 12 pulses, see Appx. C1 for details.

V. APPLICATION: QUANTUM SENSING WITH INTERACTING SPIN ENSEMBLES
Quantum sensing presents additional challenges beyond the simple decoupling of the effects that cause decoherence. Here, in addition to decoupling disorder and interactions to extend coherence time, one also needs to recouple the target signal to perform effective sensing. While there has been extensive research for quantum sensing with non-interacting systems (e.g. the XY-8 sequence [28,36]), there are only a limited number of such demonstrations for strongly-interacting systems [26,110], despite the pressing need for such protocols to further improve sensitivity in high density spin ensembles [111][112][113][114]. Here, we show that our framework addresses these challenges, by designing robust AC field-sensing pulse sequences that achieve maximal sensitivity to the target signal, while decoupling on-site disorder and spin-spin interactions. Furthermore, our formalism also provides a systematic approach to attain optimal sensitivity under given constraints, allowing diverse sensing strategies optimized for different scenarios.

V.1. General Formalism for AC Magnetometry
To achieve AC field sensing using interacting spin ensembles, we first incorporate external AC signals into the average Hamiltonian analysis. Specifically, the extra Hamiltonian due to the external AC signal can be modeled as for each sequence are shown with resonant AC signals (blue curves) at the frequency f0 to be detected. The sequence length, T , includes both free evolution periods and finite pulse durations. The π 2 and π pulses are depicted in the same way as in the pulse legend of Fig. 4. (c) Effective magnetic field ⃗ B eff generated at the resonance frequency f0 for XY-8 (top) and Seq. C (bottom). The initial state, ψ0⟩, needs to be prepared perpendicular to the direction of ⃗ B eff to achieve the best sensitivity. (d) Resonance spectra for XY-8 (left) and Seq. C (right), with their total spectral intensities (green) and individual spectral intensities for each axisx (blue),ŷ (red) andẑ (yellow). For Seq. C, the largest peak giving optimal sensitivity is identified at f = f0 [resonance 1 ], where all three axes act in phase. If the three-axis phases are mismatched, then the corresponding total peak height reduces, as seen at f = (2 3)f0 [resonance 2 ], despite a higher individual axis response. (e) Illustration of three-axis spectral-phase synchronization at resonance 1 , and out-of-synchrony at resonance 2 , resulting in a stronger and weaker effective magnetic field strength (black arrows), respectively.
where γ is the gyromagnetic ratio of the spins, B AC , f and φ are the amplitude, frequency, and phase of the target AC signal, respectively.
For a given pulse sequence represented as F = [F µ,k ], we can apply the same average Hamiltonian analysis to understand how driven spins experience the sensing field in their effective frame, giving for µ = x, y, z, and Re denotes the real part. Physically, the time-averaged sensing-field Hamiltonian [Eq. (25)] has a simple and elegant interpretation: in the toggling-frame picture, all driven spins will undergo a coherent precession around the effective magnetic field ⃗ B eff . Additionally, as seen in Eq. (26), the orientation and strength of ⃗ B eff are determined by the frequency-domain resonance characteristics of the applied pulse sequence, The AC magnetic field sensitivity η ac (f ), characterizing the minimum detectable signal strength for an AC signal at frequency f , scales as where F t (f ) is the total spectral response at the resonance frequency f under the pulse sequence and T 2 is the coherence time of the spin ensemble. Physically, F t (f ) can be understood as the effective signal strength experienced by the driven spins at resonance, namely F t (f ) = ⃗ B eff (f ) B AC , given by for a more general expression), with the equality saturated when φ =φ x (f ) =φ y (f ) =φ z (f ). Thus, it is crucial to align and synchronize the spectral phases of the pulse sequence at the target frequency f to the phase of the sensing signal, in order to achieve the best sensitivity. In such a phase-synchronized case, the effective magnetic field becomes To optimally detect this effective sensing field, spins then need to be initialized perpendicular to ⃗ B eff to form the largest precession trajectory and maximize signal detection contrast. In addition, to optimize contrast for a projective measurement along theẑ-axis, for readout the precession plane should be rotated to contain theẑ-axis.

V.2. Examples
To better understand the above abstract analysis, we analyze detailed examples of different sensing pulse sequences designed for quantum sensing with non-interacting and interacting systems.
For non-interacting systems, the XY-8 sequence, consisting of periodic π pulses, is widely employed for AC magnetic field sensing [ Fig. 5(a)]. In our representation, the XY-8 sequence is written as F x = F y = 0, F z = [+1, −1, +1, −1, ⋯], and τ = [τ, τ, τ, τ, ⋯], with intermediate zero-duration frames more explicitly shown in Fig. 5(a). Due to the periodic sign modulation present in F z , its Fourier transformF z has a peak at the frequency f 0 [ Fig. 5(d), left], indicating that an AC signal oscillating at f can be efficiently detected if the pulse spacing τ is tuned to satisfy the resonance condition f 0 = f . According to Eq. (30), the resonant AC signal is then transformed to an effective DC field along theẑ-axis where . This is the theoretically maximal achievable effective field strength in periodically driven systems for an AC signal. The initial spin state should be prepared on the xy plane, orthogonal to theẑ-axis, for optimal sensitivity [ Fig. 5(c), top]. For interacting systems, however, the XY-8 pulse sequence is not ideal as it does not address interaction effects, resulting in rapid decoherence. Instead, we design a new sequence, Seq. C as shown in Fig. 5(b), optimized for quantum sensing based on interacting ensembles. Using our framework, Seq. C is also designed to satisfy all robust decoupling conditions provided in Tab. 1, indicating that it can significantly extend the coherence time of polarized initial states, especially for interacting systems with strong on-site disorder (T W < T J in Seq. C). More importantly, the frame-matrix representation of Seq. C reveals that it also has the same periodic sign modulation in F, this time along all three axes, suggesting that an external AC signal at f 0 can again be detected efficiently [ Fig. 5(b)]. Moreover, Seq. C shows identical resonance characteristics between all three axes, provid- For the individual spectral functions of F x,y,z , resonance 2 shows a higher peak intensity than resonance 1 , while in the total spectral strength F t , the resonance at f = f 0 becomes the strongest peak, giving maximum sensitivity [ Fig. 5(d), top right]. The significant intensity reduction in F t at resonance 2 originates from asynchronous phasesφ x ,φ y and φ z [ Fig. 5(e)], leading to a smaller total field strength as a result of destructive field addition [Eq. (29)]. However, at the resonance 1 , the spectral phases of all three-axes are aligned (φ x =φ y =φ z , Fig. 5(e)) and thus generate the largest effective sensing field along the [1, This new effective field direction also necessitates an unconventional sensor initialization and readout direction, as illustrated in the lower panel of Fig. 5(c): instead of preparing the spin in the usual initial states along thex-orŷ-axes, achieving maximal sensitivity requires preparing the spin in a plane orthogonal to the effective field direction. Moreover, to maximize signal contrast during readout, it is necessary to rotate the precession plane to a position that contains theẑ-axis. The requirements of equal evolution time along each axis for interaction symmetrization [Condition 2] imply that the single-axis effective field strength of Seq. C, B Seq. C , is three times smaller than that of XY-8, B XY-8 [ Fig. 5(b)]. This is partially compensated by aligning the resonance frequency and phase along each axis, as shown in Seq. C, leading to a total effective vector field strength that is only √ 3 reduced. Despite this penalty, we still expect Seq. C to outperform XY-8 in AC magnetometry for interacting spin ensembles, since the extension in coherence time T 2 due to interaction suppression (recall η ac ∝ 1 √ T 2 ) can more than counteract the √ 3 sensitivity loss [114].
In addition to the scenarios that we have described here, our average Hamiltonian approach also helps to identify other undesired effects, such as spurious harmonics [115], which appear as additional spectral resonances in the total modulation function for finite pulse duration. This clearly demonstrates the utility of our formalism for the design of quantum sensing pulse sequences in the presence of interactions, disorder, and control imperfections.

V.3. Design Considerations for Efficient Quantum Sensing
The additional requirements of optimizing magnetic field sensitivity impose new algebraic constraints within our framework. Here, we discuss the implications of these new constraints on the structure of sensing pulse sequences by utilizing the techniques described in Sec. IV.2.
For efficient quantum sensing and to decouple on-site disorder as rapidly as possible, it is desirable to maintain a periodic structure in which the free evolution periods have frame directions that alternate between +1, -1, as adopted in XY-8 as well as Seq. C [ Fig. 5(a,b)]. However, for interacting ensembles where interaction-symmetrization is performed, this implies that any interface between two frame orientations will always have a fixed odd parity, and will thus violate condition 3 in Tab. 1 if single π 2 pulses are used for the frame transformations. Thus, to preserve sequence robustness, it is necessary to use composite pulse structures in which each frame-switching rotation is realized by a combination of two π 2 pulses to intentionally inject even parities to counteract the odd parities. An example of such a composite pulse is shown in Fig. 5(b).
Moreover, while the effective [1,1,1]-directional field presented in the previous section is close to optimal for interacting ensembles, its strength can be further improved by adding an imbalance in the effective phase accumulation along each axis. Although the sum of the phase accumulation along all axes is fixed, the effective field strength depends on the sum of squares of the phase accumulation [Eq. (29)]. Thus, due to this nonlinearity, the effective field strength can be increased when the phase accumulation is different along the three axes, which is achieved by choosing the frame along one of the axes to be at the maxima of the sinusoidal sensing signal, resulting in enhanced phase accumulation (see Appx. C2 and Seq. I in Fig. 8 for details).

VI. APPLICATION: QUANTUM SIMULATION WITH TUNABLE DISORDER AND INTERACTIONS
Our framework can also be readily adapted to engineer various Hamiltonians in the context of quantum simulation. Here, the goal is to realize different types of interactions with tunable on-site disorder via periodic driving [81-85, 116, 117], such that one can explore a range of interesting phenomena in out-of-equilibrium quantum many-body dynamics, including dynamical phase transitions [42][43][44], quantum chaos [45,46] and thermalization dynamics [54][55][56][57][58][59][60]. Moreover, the interplay of disorder, interactions and periodic driving can also lead to novel nonequilibrium phases of matter, such as the recently-discovered discrete time crystals [47-53, 118, 119].
Indeed, as can be seen in Eqs. (7-10), we can design the toggling-frame spin operatorsS z k = ∑ µ F µ,k S µ to achieve a nonzero target sum in Eqs. (5,6) and engineer the leadingorder average Hamiltonian H avg . More specifically, we show that for the common form of two-body interaction Hamilto-  1 3), XY (c = 1 2) and dipolar interactions (c = 1), can be realized by varying the single parameter c which is defined as the relative fraction of the total time evolution along theẑ-axis in the toggling frame. In this example, the system is chosen to be a disordered dipolar interacting spin ensemble with J S ij = −J I ij and J A ij = 0, naturally realized with NV centers in diamond [78]. The dashed line indicates the maximum effective disorder strength W eff achievable for a given interaction type, tuned by the single parameter c as where c captures the imbalanced time evolutions in the toggling frames, defined such that the system evolves under thê z-andx,ŷ-axes toggling frames for total durations cT and (1 − c)T 2 in one sequence cycle T . Thus, the Floquetengineered interaction Hamiltonian H int avg now exhibits modified Ising and exchange interaction strengths ofJ I ij andJ S ij , respectively. Taking the case of interacting NV spin ensembles as an example [78], where J S ij = −J I ij , we can continuously interpolate between Ising (c = 0), Heisenberg (c = 1 3), XY (c = 1 2), and dipolar-like (c = 1) interactions by tuning the proportion c of the sequence. Moreover, on-site disorder can also be independently controlled by introducing an additional sign imbalance along each axis in the toggling frame, changing the original disorder Hamiltonian where the effective disorder field ⃗ h eff can now have both longitudinal and transverse field components.
We illustrate the accessible range of disorder and interaction Hamiltonians with this scheme in Fig. 6(a). Note that the maximum effective disorder strength W eff is dependent on c (dashed line in Fig. 6(a), see Sec. S1H [100]). Two representative examples of how to engineer such interaction Hamiltonians in a robust fashion are shown in Fig. 6(b).
Combined with the techniques for robust engineering of other terms in the Hamiltonian, such that imperfections are suppressed, this allows access to a broad range of interacting, disordered Hamiltonians that potentially exhibit very different thermalization properties [54][55][56][57][58][59][60]. Thus, our framework will open up a new avenue for the robust Floquet engineering of many-body Hamiltonians.
Here, we focus on the experimental implementation and demonstration of our results in an interacting ensemble of NV centers in diamond [ Fig. 7(a)], see Sec. S2A [100] and Refs. [78,108,114] for more details of our sample and experiments. It is characterized to be a disorder-dominated system [ Fig. 1(a)], exhibiting large on-site disorder (W ≈ (2π)4 MHz) with modest interaction strengths (J ∼ (2π)35 kHz), corresponding to the regime in Fig. 4(c) where Seq. B is expected to perform significantly better than Seq. A . We experimental confirm this in Fig. 7(b), where Seq. B shows a significantly longer coherence time compared to Seq. A.
In Fig. 7(c,d), to illustrate the importance of fulfilling the robust decoupling criteria [Tab. 1], we compare the robustness of two different sequences to systematic rotation angle deviations. Here, we design a non-robust Seq. D, which is almost identical to the robust Seq. B, but does not suppress spin-rotation angle errors; In Seq. D, intermediate frames are intentionally chosen to violate the suppression condition for rotation angle errors [Condition 4], while satisfying the rest of the conditions (see Fig. 8 for more details of the sequence). Fig. 7(c) illustrates the coherence decay profile of driven spins under these two sequences when the rotation angle is chosen to be 92.5% of the correct rotation angle (gray dots in Fig. 7(d)); Seq. B does not show any oscillations, while Seq. D shows pronounced oscillations over time, resulting from a residual error term δH avg (see Sec. III.2). This behavior is further confirmed in Fig. 7(d), where we extract the effective modulation frequency of the spin coherence as a function of the systematic rotation angle deviation. While Seq. B does not show any oscillations, Seq. D shows a linear dependence of oscillation frequency with the rotation angle error, indicating that it is not robust against perturbations.
Finally, we illustrate the quantum sensing applications of our framework, focusing on the importance of relative spectral phase alignment of the resonances in the three axis directions. In Fig. 7(e), we show the change in signal detection contrast under the presence of an AC magnetic field with varying frequency, for the robust sensing sequence Seq. C described in Sec. V and Fig. 5. As discussed in Sec. V, the resonances 1 and 2 have comparable field strengths along each individual axis (as an example,ŷ initialization is plotted in Fig. 7(e), red trace), but the total field strength is larger for resonance 1 , where the individual axis responses are phase-synchronized with each other (Fig. 7(e), green trace). We confirm this phase-sensitive response in Fig. 7(f), where the contrast is measured as a function of the sensing signal phase φ for three different initializations along thex-,ŷ-andẑ-axes. For resonance 1 , we see that the responses of all three initializations become phase-synchronized, while for 2 , the curves are outof-phase with respect to each other. As a result of the spectral phase alignmentφ x =φ y =φ z at the resonance 1 , a large effective sensing field is generated along the [1,1,1]-direction. This can yield a high magnetic field sensitivity, surpassing the interaction-limited value of the conventional XY-8 sequence, as we demonstrate in Ref. [114]. These results corroborate the quantum sensing capabilities of interacting spin ensembles as well as the robustness of the phase-aligned sensing sequence, which maximizes AC sensitivity while suppressing spin-spin interactions and on-site disorder.

VIII. DISCUSSION AND CONCLUSION
In this paper, we have introduced a novel framework for the efficient design and analysis of periodic pulse sequences to achieve dynamic Hamiltonian engineering that is robust against the main imperfections of the system. Our approach provides versatile means to design and adapt pulse sequences for a wide range of experimental platforms, by considering their system characteristics such as disorder, interactions and control inhomogeneities. Key to our approach is the adoption of a toggling frame description of the sequence and the resulting average Hamiltonian. Crucially, we find that various types of leading-order control errors can be systematically described by the time-domain transformations of a single interaction-picture Pauli spin operator during free evolution periods. This allows us to derive a simple set of algebraic conditions to fully describe all necessary conditions for specific target applications, significantly simplifying the design of pulse sequences. Remarkably, these algebraic conditions also allow the construction of efficient strategies and the proof of their optimality to enhance various figures of merit, such as sequence length and sensitivity. Using a dense ensemble of interacting electronic spins in diamond, we experimentally confirm the wide applicability of our framework in protecting quantum information as well as detecting external signals with high sensitivity.
In addition to its wide-reaching consequences on the systematic design and analysis of pulse sequences for various applications, our framework also opens up a number of intriguing directions for future studies. For example, we can extend our approach to higher-spin systems to investigate more complex quantum dynamics, such as quantum chaos and information scrambling, as well as utilize larger effective FIG. 7. Experimental demonstration. (a) Sample used in experiments, containing strongly disordered, interacting NV centers in black diamond. We isolate two levels from the spin-1 ground state, { 0⟩ , ±1⟩}, by a static Zeeman shift (red arrow), to define an effective spin-1/2 (qubit) system, { 0⟩ , −1⟩}. Resonant pulsed driving at frequency ω0 is applied to the spin to manipulate its quantum states. In the sample, NV centers at a distance interact with one another via magnetic dipolar interactions (circled diagram). (b) Comparison of coherence decay of driven spins under Seq. A (Cory-48, red) and Seq. B (blue). The pulse spacing is fixed to τ = 25 ns, and the 1 e-decay times T2 are extracted from a stretched exponential fit, giving T A 2 = 1.2 µs and T B 2 = 7.9 µs for Seq. A and Seq. B, respectively. (c) Sequence robustness against spin-manipulation error. In the presence of systematic spin-rotation angle deviation 0.925(π 2) (gray points in (d)), the non-robust Seq. D (right) shows modulations in the coherence profile for all three initial states polarized along thex (blue),ŷ (red) andẑ (yellow) axes; the robust Seq. B is insensitive to such errors, corroborated by the absence of the modulation. (d) Modulation frequency for Seq. B and Seq. D as a function of systematic rotation angle error relative to the perfect π 2 rotation, Ω0tp = π 2. The small lateral offset observed in Seq. D at zero modulation frequency is due to a slight Rabi frequency calibration error. (e) AC magnetic-field sensing using the robust Seq. C, presented in Fig. 5, for initial states prepared along theŷ-axis (red) and optimal [−1, −1, +2]-axis (green). Solid lines are theoretical fits from the spectral intensities obtained from Seq. C, showing the largest peak at f = f0 [ Fig. 5(b)]. (f) AC sensing response under Seq. C as a function of external AC signal phase φ for signal frequencies at f = f0 (resonance 1 ) and at f = (2 3)f0 (resonance 2 ). We observe in-phase and out-of-phase responses at resonance 1 and 2 , respectively. See Sec. S2B [100] for details of the fitting functional forms. dipoles in those high-spin systems for more effective sensing [83,84,[129][130][131]. Higher-order contributions beyond the leading-order average Hamiltonian can also be systematically incorporated using the proposed framework, as we discuss in Appx. D. In addition, our formalism may also be extended to the synthesis of dynamically-corrected gates and other nontrivial quantum operations [16,132,133]. While we have focused on the case of π 2 and π pulses around x,ŷ axes for simplicity, it will be interesting to extend the analysis to more general control pulses, which could enable shorter protocols for Hamiltonian engineering. Moreover, by employing optimal control techniques to further boost the performance of the pulse sequences [27,40,134,135], we may be able to robustly engineer many-body Hamiltonians to create macroscopically entangled states, such as spinsqueezed states or Schrödinger-cat-like states, to be used as a resource for interaction-enhanced metrology beyond the standard quantum limit [136,137].
Here, we introduce the basic principles of AHT and start by considering a generic time-dependent Hamiltonian for a driven quantum system where H s is the system Hamiltonian governing the internal dynamics and H c (t) describes the time-dependent control field used to coherently manipulate the spins (qubits). For a Floquet system, the control field is modulated in time with a periodicity of T , i.e., H c (t) = H c (t + T ). At times t = N T , the many-body state is given by ψ t=N T ⟩ = U(T ) N ψ 0 ⟩ with the interaction picture unitary evolution operator [1] where T denotes time-ordering. Here,H s (t) is the rotated system Hamiltonian in the interaction picture with respect to control fields, given byH The control unitary rotation operator over one period is chosen to be identity U c (T ) = I.
AHT allows the identification of a time-independent effective Hamiltonian H eff such that The Magnus expansion of U(T ) with expansion parameter T [138] can be used to approximate this effective Hamiltonian as where l is the truncation order andH (k) is the k-th order contribution in the Magnus expansion. The first two terms in the series arē s (t 1 )dt 1 and (A5) Although the accuracy of the average Hamiltonian approximation depends on the truncation order l, if the Floquet driving frequency 1 T is much faster than the local energy scales associated with the system Hamiltonian H s , then the first few terms are sufficient to model H eff and approximate the dynamics of the many-body state to an accuracy improving exponentially in l [139][140][141][142]. In the following, we focus on the leading order contribution, corresponding to only retainingH (0) in the series. A general control field consists of n pulses {P k=1,⋯,n } with nonuniform pulse spacing {τ k=1,⋯,n }, as shown in Fig. 1(c). Each P k defines a pulsed unitary rotation, generating a discrete set of rotated Hamiltonians {H k=1,⋯,n }, wherẽ withH 1 = H s . As the interaction-picture Hamiltonian is rotated (toggled) at every pulse,H k are also referred to as the "toggling-frame Hamiltonians" and govern the spin dynamics in their respective free evolution intervals τ k . For infinitely short pulses, the zeroth-order average Hamiltonian, H (0) = H avg , can be simplified from an integral to a weighted average of the toggling-frame Hamiltonians, as presented in Eq.
(1) of the main text.

Appendix B: Analysis of Three-Body Interactions
In this appendix, we analyze the decoupling conditions for spin-1/2 three-body interactions. Interestingly, our versatile formalism can be applied to these more complex scenarios, leading us to a new set of rules that allow for the implementation of robust protocols in the presence of three-body interactions. First, we show that in the limit of ideal pulses, the decoupling conditions described in Sec. II.2 are also sufficient to fully suppress dynamics under any secular three-body interaction for a polarized initial state. We then describe an extension of the formalism that accounts for finite pulse duration effects in the presence of three-body interactions, leading to new decoupling conditions beyond those discussed in Sec. III.
While most naturally occurring physical systems only involve two-body interactions, interactions involving more particles can lead to a number of exotic physical phenomena. For example, fractional quantum Hall state wavefunctions appear as the ground state of Hamiltonians involving three-body interactions [143,144], and many other topological phases and spin liquids are ground states of such many-spin Hamiltonians [145,146]. There have also been various proposals for the direct realization of three-body interactions in experimental platforms ranging from cold molecules [101] to superconducting qubits [102,103]. They may also emerge in the form of a higher-order term in the Magnus expansion of a system with only two-body interactions.
As a first step towards the control and engineering of such interactions, we analyze the conditions for dynamical decoupling for a polarized initial state. As in the main text, we will be focusing our attention on interactions under the secular approximation, where all terms in the Hamiltonian commute with a global magnetic field in theẑ-direction.

Ideal Pulse Limit
Before we discuss the details of three-spin coupling, we first prove a useful lemma for interactions under the secular approximation in the perfect, infinitely short pulse limit.
Lemma: For any interaction under the secular approximation, averaging over the spin-1/2 single qubit Clifford group is equivalent to averaging over toggling frames ofS z that cover the six axis directions ±S x,y,z .
Proof: Consider a generic N -body interaction Hamiltonian H and a set of unitary operators U k (k = 1, 2, ⋯, K). The average Hamiltonian over this set is given by Let us now group the elements of the Clifford group into sets defined by how the elements transform the S z operator. Each set contains elements that satisfyS z = U † S z U = (−1) ν S µ with ν = 0, 1, while the rotatedx-axis spin operator,S x = U † S x U , can take four distinct values that are orthogonal to theS z direction. In our toggling-frame representation, however, any of the four Clifford elements in the same set will correspond to a single term specified byS z . Thus, proving the lemma reduces to proving that the four Clifford elements above give identical Hamiltonians. We prove this by observing that for any two elements U 1 and U 2 in the same set, there exists a rotation U z around theẑ axis such that U 1 = U z U 2 (this rotation leaves the interaction pictureS z invariant, but changesS x ). The Hamiltonian under conjugation by U 1 is then given by where we use (U † z ) ⊗N HU ⊗N z = H. This holds because a rotation around theẑ-axis does not modify the secular Hamiltonian, which commutes with the global S z operator. Consequently, a conjugation of the average Hamiltonian above by U 1 will be equal to a conjugation by U 2 , and thus each set of Clifford elements that transform the S z operator in the same way will result in identical Hamiltonians. ∎ With this lemma in hand, we now utilize mathematical results from unitary t-designs [105][106][107][147][148][149][150][151][152] to prove that after symmetrizing along the six axis directions (i.e. the ±x, ±ŷ, ±ẑ-axes), a polarized initial state will be an eigenstate of the resulting symmetrized Hamiltonian. A unitary t-design is a set of unitary operators {U k }, such that for every polynomial P (t,t) (U ) of degree at most t in U and at most t in U † , the average over {U k } is equivalent to the average over the Haar measure of all unitaries of the same dimension [105]. More specifically, the following relation holds: Here, U(2) is the unitary group of dimension 2, used to describe two-level systems, O is an N -body operator with N ≤ t and O U is the corresponding averaged observable. This means that for observables up to order t, the effect of averaging over the finite set of unitary operators {U k } is equivalent to averaging over all unitaries of dimension 2. The symmetrizing properties of the right-hand-side of Eq. (B3), where the average is taken over all elements of the unitary group over the Haar measure, imply that O U must only contain terms proportional to elements of the symmetric group S t of order t [150]. This is because all other terms will be transformed and symmetrized out by the average, but elements of the symmetric group, which only permute the labels of the states, will be invariant, as the unitary operator U ⊗N conjugates all spins identically.
It is known that the Clifford group forms a unitary 3design [106,107]. Combined with the lemma presented above, this implies that for any sequence that equally covers all six axis directions, all interactions involving three particles or fewer will be symmetrized into a form that only contains terms proportional to elements of the symmetric group. Any initial state with all spins polarized in the same direction will then be an eigenstate of this symmetrized interaction, since this state is invariant under any permutation of the elements. Correspondingly, a polarized initial state does not experience decoherence under this interaction.
As a nontrivial example of this result, let us consider the interaction H int = J(S x S y S z − S y S x S z ). The symmetrized Hamiltonian can be calculated to be H int = J 3 ∑ µνσ µνσ S µ S ν S σ , where µνσ is the Levi-Civita symbol. One can explicitly verify that any globally polarized initial state is an eigenstate of the symmetrized Hamiltonian with eigenvalue 0.
As the Clifford group is not a unitary 4-design, four-body interactions will still induce dynamics after symmetrization. Indeed, we can explicitly verify this by considering the symmetrized interaction (S x ) ⊗4 + (S y ) ⊗4 + (S z ) ⊗4 , which is found to act nontrivially on a generic polarized initial state.

Finite Pulse Duration Effects
We now illustrate how to analyze finite pulse duration effects for three-body interactions using our sequence representation matrix F. We consider generic interaction Hamiltonians with up to three-body interactions and, in analogy to Eqs. (7-10), we write the k-th toggling-frame Hamiltonian as a polynomial in F µ,k :H where O µ,l describes the generic operator form of interactions that transform as the l-th power of F µ,k . More specifically, O µ,l can be written as a sum of terms, each composed of a product of Pauli operators that preserve the total magnetization and thus include an even number of S x or S y operators. Consequently, the interaction must either be of the form S z S z S z , or involve the tensor product of an S z operator and a polarization-conserving two-body operator, which can be S x S x + S y S y or i(S x S y − S y S x ). Each of these terms can thus be written as a product of individual components that transform as F µ,k . For example, we can rewrite the following three-body interaction as This is a three-body interaction with l = 2, since it is proportional to the square of F µ,k . This suggests that during the finite pulse duration between free evolution blocks k and k + 1, where the interaction-picture operatorS z (θ) = ∑ µ (cos θF µ,k +sin θF µ,k+1 )S µ with θ evolving from 0 to π 2, the corresponding average Hamiltonian can be written as whereÕ contains operators acting during the rotation pulse and is implicitly dependent on l and the indices in the bracket (for example, when θ = 0,Õ = O µ,l ). Expanding the polynomial in the bracket of Eq. (B6), we shall find contributions corresponding to terms of degree u in F µ,k and v in F µ,k+1 , with u + v = l, since F µ,k and F µ,k+1 each have only one nonzero element. This allows us to generalize the conditions for decoupling finite pulse imperfections discussed in the main text to three-body interactions, and also provides an alternative perspective to the conditions in the main text. For interaction terms exhibiting linear (l = 1) or quadratic (l = 2) dependence on F µ,k , the decoupling conditions presented in Tab. 1 of the main text can be directly applied. Similarly, for l = 3, the terms F 3 µ,k and F 3 µ,k+1 directly correspond to the three-body interactions appearing in the original Hamiltonian within the free evolution periods k and k +1, and can be easily incorporated into the sequence design by extending the effective duration of the free evolution time. Meanwhile, the decoupling of cross terms F 2 µ,k F ν,k+1 and F µ,k F 2 ν,k+1 correspond to a generalization of the interaction cross-term decoupling condition [Condition 3 in Tab. 1] described in the main text: for each pair of (µ, ν). As an example, for the pair of directionsx,ẑ, we consider all instances in which thex andẑ frames appear in the free evolution frames immediately preceding and following a π 2 pulse, and the above decoupling condition requires that the signs of allx frames appearing in such positions sum up to 0. Combining these results, we see that our formalism provides a systematic method to robustly decouple the effects of any three-body interaction under the secular approximation, on any polarized initial state, even in the presence of finite pulse durations.
Appendix C: Efficient Sequence Design Strategies

Minimal Length for Robust Dynamical Decoupling
In this section, we discuss the minimal sequence lengths required to satisfy different combinations of decoupling conditions, and provide examples of pulse sequences that achieve these minimal lengths.
We start by considering the minimal number of free evolution blocks to fully symmetrize the interaction Hamiltonian. From condition 2 of Tab. 1, we see that the minimal nontrivial solution requires nonzero elements in each of the three rows of F. Consequently, at least 3 free evolution blocks are required. A possible realization is (see also Seq. E in Fig. 8) However, according to the discussion in Sec. S1G [100], a single cycle evolution under this sequence does not return the transverse spin operators to their original configuration, i.e. U c (t) ≠ I. To achieve U c (t) = I, we can use with composite π 2 pulses inserted at the interfaces of free evolution intervals (see Seq. F in Fig. 8). However, neither of the above sequences are robust against on-site disorder or finite pulse duration effects. As discussed in Sec. IV.2 in the main text, to fully symmetrize interactions and cancel disorder, at least 6 free evolution periods are required. If composite pulses are allowed, there is a pulse sequence consisting of 6 free evolution periods that also satisfies all robustness requirements, as described in Sec. IV.2 the main text. In some scenarios, however, composite pulses may be undesirable due to the technical challenges of implementing independent pulses in quick succession; in this case, we can show that at least 12 free evolution intervals are required. Let us first examine why 6 intervals are not sufficient: The cross-interaction parity condition [Condition 3] requires the number of odd-parity frame changes to be equal to that of even-parity frame changes; a sequence with 6 free evolution intervals thus gives 3 even-parity and 3 oddparity frame changes. The odd number of odd-parity frame changes, however, results in a toggling-frame operator at the beginning of the next cycleS z (T ) that has opposite sign from S z (0), violating the periodic conditionS z (0) =S z (T ). To prevent this, we thus require an even number of odd-parity frame changes. Since the length of the sequence has to be an integer multiple of 6 to simultaneously accomplish disorder and interaction decoupling, the minimal sequence length is at least 12 frames, which now realizes a valid Floquet cycle with U c (T ) = I and satisfies all decoupling conditions.
One realization of such a pulse sequence, for example, can be written in our representation as (see also Seq. G in Fig. 8) The above argument can also be readily extended to other scenarios. For disorder-dominated systems it is desirable to maintain a fast spin-echo structure on the toggling-frame evolution of a sequence, as discussed in Sec. IV. This necessitates FIG. 8. Complete representation of periodic pulse sequences. All sequences are shown in the conventional pulse representation and our toggling-frame transformation-based representation. In the toggling frame representation, squares indicate free evolution times of length τ and narrow lines indicate short intermediate frames; yellow is positive and green is negative. Seq. A: The Cory-48 sequence, which decouples interactions on a faster timescale and disorder on a slower timescale, and is robust against leading-order imperfections; Seq. B: Pulse sequence designed to decouple disorder on a faster timescale and interactions on a slower timescale, and is robust against leading-order imperfections; XY-8: Standard pulse sequence for dynamical decoupling with non-interacting spins; Seq. C: Pulse sequence to illustrate robust dynamical decoupling of interactions and disorder, and well-aligned sensing resonances; Seq. D: Pulse sequence that has identical free evolution frames as Seq. B, but with intermediate frames in the second half permuted to remove the robustness against rotation angle errors. Seq. E: Simple sequence to fully symmetrize interactions; however, the transverse spin operators do not return to themselves at the end of this sequence; Seq. F: Modified version of Seq. E, in which the composite pulses ensure that the transverse spin operators return at the end of the sequence; Seq. G: Minimal pulse sequence that meets all dynamical decoupling and robustness requirements (Tab. 1 in the main text) without composite pulses; Seq. H: Minimal pulse sequence that meets all dynamical decoupling and robustness requirements without composite pulses, and has a fast spin echo structure; Seq. I: Minimal pulse sequence that achieves optimal vector sensitivity under interaction-decoupling constraints, and illustration of phase accumulation along each axis.
a frame matrix structure in which the frames along each axis always consist of pairs with opposite sign. When permitting only a single π 2 pulse to switch frames, the parity product of the last element of the first pair and the first element of the second pair will again have the same parity constraints as the preceding case. Thus, we can apply the same argument to show that the 12-frame sequence also needs to be doubled to satisfy the parity condition if composite pulses are not allowed and a spin-echo structure is required on a fast timescale. One example that achieves this is illustrated as Seq. H in Fig. 8.

Composite Pulses for Sensing
We now discuss the implications of the algebraic conditions in Tab. 1 in the context of quantum sensing. Here, we show that for pulse sequences that follow a periodic signmodulation structure along each axis (for instance, the sequences in Fig. 5 of the main text), to fully utilize the effective sensing field and satisfy the robust dynamical decoupling conditions, it is necessary to employ composite pulses for the effective π 2-pulse implementations to suppress interaction cross-terms [Condition 3].
If only single π 2 and π pulses are employed to connect different axes in the toggling frame, the fixed spin-echo-type sign modulation patterns inevitably result in fixed parities at every interface between two axis directions, leading to the violation of condition 3 in Tab. 1 (see also discussion in Appx. C1). To address this, composite pulses should be utilized instead to connect different axes and adjust the parities, allowing one to balance the number of even and odd parity interfaces.

Optimal Sensitivity for Interacting Spin Ensembles
Our formalism also allows us to easily determine the optimal achievable effective field strengths under the condition that interaction effects are decoupled.
As described in the main text, the condition for Ising and spin-exchange interaction-decoupling requires equal evolution times along each of thex,ŷ,ẑ axes in the toggling frame. In the typical scenario where a spin echo is performed along each axis, an effective sensing field ⃗ B eff = B XY-8 1 3 , 1 3 , 1 3 is generated at the dominant resonance, with a magnitude a factor of √ 3 smaller than the maximal field strength B XY-8 achievable for the XY-8 sensing sequence, where interactions are not decoupled. In fact, despite the relative reduction in the sensing field magnitude, our sequence is optimal for sensing with disorder-dominated interacting spin systems, as interaction-decoupling requires equal evolution time along different axes, and a fast spin-echo structure is required to not only effectively detect external AC signals but also rapidly suppress the effects of disorder.
If we broaden the scope to include pulse sequences that do not have an immediate echo structure, then the above pulse sequence can be further improved into a version that is provably optimal in the limit of infinitely short pulses. Here, the key observation is that the amount of phase accumulation along each axis can be imbalanced, even while spending equal amount of time along all axes to decouple interactions. We can theoretically prove that the effective field strength is maximized by making the individual phase accumulations along each axis as different from each other as possible, as shown in Seq. I in Fig. 8. The detailed derivation for its optimality is provided in Sec. S1I. In this case, the corresponding effective field is 1 2 , which gives an optimal field strength for interaction-decoupling sequences, resulting in a 10% improvement over fast-echo-based sensing sequences, such as Seq. B, that require an even distribution of phase accumulation between the three axes.

Appendix D: Suppression of Higher Order Effects
While our analysis has mainly focused on the zeroth-order average Hamiltonian, we can readily incorporate higher-order effects by considering the full Magnus expansions [Eq. (A4)] to engineer effective Hamiltonians with higher accuracy. Developed primarily in the NMR community, there are several known approaches, such as reflection-symmetric pulse arrangements [25,26,94,153] and concatenated sequence symmetrization [13,[154][155][156], to suppress higher-order contributions in driven spin dynamics. In the following, we discuss how these techniques can also be naturally incorporated into our sequence design framework.
Reflection symmetry: When Hamiltonians in the toggling frame respect reflection symmetry [94], that is,H n+1−k =H k , all odd-order terms in the Magnus expansion vanish,H (2l−1) = 0 with integer l [Eq. (A4)]. In our framework, this imposes an additional condition on the sequence F; Generically however, any sequence can be extended into a pulse sequence that respects reflection symmetry simply by doubling the length of the frame matrix and filling the second half with its own mirror image in time, taking care of pulse imperfections at the central interface. We note, however, that for certain applications such as quantum sensing, one needs to take additional care when performing such symmetrizations, since reflection symmetry (similar to a time-reversal operation) may accidentally cancel the desired sensing field contributions [Sec. V].
Concatenated sequence symmetrization: One way to understand and engineer higher-order Hamiltonian engineering properties of a sequence is to decompose it into smaller building blocks. A few techniques developed along these lines include pulse-cycle decoupling [25] and concatenated symmetrization schemes [13,[154][155][156], where a long pulse sequence is constructed from the repetition of short pulse sequences, symmetrized in a systematic pattern to suppress higher-order effects. Such constructions can be easily added to our framework, and we leave the detailed application of such techniques to a future work.
Second averaging: The technique of second-averaging has been developed in NMR to suppress dominant error terms in δH avg that do not commute with the leading contribution in H avg [157][158][159]. Such methods can be readily incorporated in our framework by alternating the rotation axes of control pulses periodically every Floquet cycle, or by using offresonant driving. CONTENTS S1. Formalism details S1 A. S2. Experimental Details S10 A. Experimental Setup S10 B. Data Analysis S10 References S11 S1. FORMALISM DETAILS

A. General Hamiltonian Under the Secular Approximation
In this section, we discuss the consequences of the secular approximation (i.e. rotating wave approximation under a strong quantizing field), and explain why under this approximation the transformation properties of the S z spin operator uniquely specify the Hamiltonian.
The secular approximation is widely used in a large variety of contexts. Since the strong quantizing field along theẑ axis sets a dominant energy scale in the system Hamiltonian, under the secular approximation we only retain Hamiltonian contributions that conserve the total spin polarization along theẑ axis; in other words, a change in the totalẑ-axis polarization, such as a single spin flip, is associated with a large energy cost and is thus energetically suppressed. In mathematical terms, this approximation corresponds to only keeping Hamiltonian contributions that commute with S z tot = i S z i . This allows us to easily write down all allowed terms in the Hamiltonian for one-and two-body operators in the Pauli operator basis. For one-body operators, the only nontrivial term is S z , since the other Pauli operators, S x,y , do not commute with the S z tot operator. Similarly, the only two-body operators commuting with S z tot are the following three types: S2 where S ± = S x ± iS y . Hence, a linear combination of the one-and two-body operators allowed in the secular approximation results in the following general form of many-body Hamiltonians which corresponds to Eq. (2) of the main text. More importantly, when global spin manipulation is performed, the S z tot -preserving nature of the interaction Hamiltonian under the secular approximation also greatly simplifies the representation of the Hamiltonian. To describe how a fully generic Hamiltonian is transformed after global rotations, we usually would need to specify how both S z and S x are transformed. However, due to the fact that the Hamiltonian commutes with the global S z tot operator, one can perform a rotation with S z to take the S x operator into any direction that is perpendicular to the transformed S z direction, and the Hamiltonian remains invariant. Consequently, the Hamiltonian is uniquely specified given the transformation properties of the S z operator alone. This justifies our approach in the main text, where, given a Hamiltonian under the secular approximation, we only need to examine the transformation properties of the S z operator.
While we have focused on the case of Hamiltonians satisfying the secular approximation, it would also be interesting to extend the techniques we have developed in this paper to situations where such an approximation does not hold, such as zero-field NMR [1].

B. Details of Average Hamiltonian During Free Evolution Time
In this section, we provide a detailed derivation of the various average Hamiltonian contributions Eqs. (5)(6)(7)(8)(9)(10) in the main text. The key idea is to express all Hamiltonian contributions in terms of rotationally-invariant terms and terms that only depend on the S z operator direction.
As defined in the main text, during the k-th free evolution periodS z (t) = µ F µ,k S µ . Correspondingly, we can directly substitute this into the expression for disorder and Ising interactions: where we have used the fact that F µ,k F ν,k = δ µν F 2 µ,k , since each column of the matrix F = [F µ,k ] has only one nonzero element. Using these expressions, we can also easily find the transformed interaction for symmetric spin-exchange interactions by making use of the identity S x i S x j + S y i S y j = S · S − S z i S z j = ( µ S µ i S µ j ) − S z i S z j , which gives: Finally, we derive the transformation of the anti-symmetric spin-exchange interaction. To this end, we assume that the S x Pauli spin operator is transformed toS x (t) = µ G µ,k S µ , where G µ,k satisfies the identity G µ,k G ν,k = δ µν G 2 ν,k and takes on values of {0, ±1} (See Sec. S1 G for details of how one can explicitly construct G µ,k ). Since the commutation relations between spin operators are conserved under frame transformations, the transformedS z (t) and S x (t) operators uniquely specify theS y (t) operator as: where µνλ is the Levi-Civita symbol. Based on this, we can write the transformation of the anti-symmetric spinexchange interaction term as: where we have used the identity presented above for G ν,k , as well as the fact that G ν,k has only one nonzero element, squaring to 1. Combining these expressions gives the average Hamiltonian terms in Eqs. (7)(8)(9)(10).

C. Details of Finite Pulse Width Analysis
Here, we provide more details for the derivation of the average Hamiltonian for a π/2 pulse with finite pulse width. For illustration purposes, we consider a concrete example in which the initial toggling frame withS x (0) = S x , S y (0) = S y ,S z (0) = S z is rotated, such that the rotation brings the frame from +ẑ into +ŷ. Correspondingly, in the interaction picture relative to the instantaneous control pulse, the operators are transformed asS x (θ) = S x , S y (θ) = S y cos θ −S z sin θ,S z (θ) = S z cos θ +S y sin θ, where θ smoothly changes from 0 to π/2 over the finite duration of the pulse. Plugging this into the form of the average Hamiltonian, we find that the instantaneous toggling-frame Hamiltonian during the rotation is Upon careful inspection, we see that the transformation properties of the Hamiltonian terms during the free evolution period (see Sec. II.3 in the main text) also dictate the response during the finite pulse width evolution. More specifically, the disorder and imaginary spin exchange terms, which transform linearly in F µ,k , are proportional to µ,ν (cos θF µ,k + sin θF ν,k+1 ) during the finite pulse duration. Subsequently, after integrating over the angle θ, we arrive at Hamiltonian contributions that are proportional to F µ,k and F ν,k+1 , with proportionality factor 2/π due to the integration. Meanwhile, the symmetric spin exchange and Ising terms, which transform quadratically in F µ,k , will produce contributions proportional to [ µ,ν (cos θF µ,k + sin θF ν,k+1 )] 2 during the finite pulse duration. After integrating over the angle θ, this will result in Hamiltonian contributions that are proportional to F 2 µ,k and F 2 ν,k+1 with factor 1/2, and a cross-term proportional to F µ,k F ν,k+1 with factor 1/π. Combining these considerations gives rise to the expression in Eq. (25) of the main text. While we have assumed that the Rabi frequency remains fixed between different pulses, our framework can be easily generalized to incorporate differing Rabi frequencies, simply by specifying a different pulse width for each pulsed rotation. The above derivation also provides a simple intuition for the different components of the finite pulse duration effect Hamiltonian: there will be terms that are proportional to the average Hamiltonian in the free evolution periods preceding and following the pulse, but for terms that transform quadratically, the square during the finite pulse duration will also generate cross terms. The same intuition will also prove to be useful when deriving conditions for finite pulse duration effects with three-body interactions (Appx. B).

D. Analysis of Rotation Angle Errors
We now consider the effects of control errors, in the form of a systematic rotation angle deviation. We will show that the effective Hamiltonian corresponding to the rotation angle error is well described by the chirality of frame changes, thus enabling a simple algebraic description of this decoupling condition.
To be more concrete, consider a π/2-pulse P k = i exp −i π 2 µ α µ,k S µ i , where α µ,k specifies the rotation axis of the k-th pulse (e.g. α x,k = +1 implies a +x-rotation). A rotation angle error corresponds to an actual rotation being applied, where i is the angle error for the spin at site i, assumed to be static and independent of the applied pulse P k ; equivalently, we can regard this as an error term in the engineered Hamiltonian acting during the rotation. The toggling-frame Hamiltonian corresponding to this rotation angle error can be calcu-S4 lated as where U k−1 = P k−1 · · · P 2 P 1 is the unitary rotation due to the preceding control pulses. We simplify Eq. (S10) by making use of the expressions for the initial and final frames connected via the pulse P k : Inserting U k−1 U † k−1 = I into the latter equation, we arrive at the expression While this equation looks complicated, upon closer inspection, it implies that U † k−1 P k U k−1 is nothing other than the rotation pulseP k = i exp −i π 2 µ β µ,k S µ i that in the toggling-frame representation, rotates the k-th toggling frame F k into the (k + 1)-th one F k+1 , where F k = µ F µ,kêµ . The rotation axis direction in the toggling frame specified by β µ,k can be simply obtained from a cross product between the toggling-frame orientations before and after the rotation: where β k characterizes the chirality of the frame change between neighboring k-th and k + 1-th toggling frames.
Since P k = i exp −i π 2 µ α µ,k S µ i , where µ α µ,k S µ is a rotation along thex orŷ axis and is thus proportional to a Pauli matrix, we can move the unitary U k−1 into the exponential This implies that α k = µ α µ,kêµ , the rotation axis in the lab frame, transforms to β k in the toggling frame: An alternative, more formal way to see this is that β µ,k satisfies the expression Using the Baker-Campbell-Hausdorff formula, the above expression can be simplified to

S5
which gives the same result as Eq. (S14). Summing the total contributions of all imperfect rotations, and using Eq. (S16), we obtain the error Hamiltonian, δH rot avg , due to the rotation angle errors which corresponds to Eq. (23) of the main text. We can also formulate a more intuitive geometric picture of the preceding derivation, by considering the Bloch sphere picture of interaction picture spin operators. To better understand this, let us return to the Bloch sphere representation shown in Fig. 2(a-c), lower panels of the main text. We shall call the coordinate system generated by the initial orientation of thex,ŷ,ẑ axes the fixed external coordinate system, and the coordinate system generated by the transformed axis directions as the toggling-frame coordinate system. Physically, the axis direction of the toggling-frame coordinate system that points along theμ-axis of the external coordinate system is the transformed S µ operator: for example, in Fig. 2(b), the second Bloch sphere has theŷ axis pointing up along theẑ direction, indicating thatS z (t) = U † k−1 S z U k−1 = S y at this time point. With this picture in mind, we can now pictorially understand the above technical derivation. At the beginning of the sequence, we start with a toggling-frame coordinate system that coincides with the external reference coordinate system, as shown in the first Bloch sphere of Fig. 2. We then apply rotations along the axes of the fixed external coordinate system, which transform the spin operators and rotate the toggling-frame coordinate system. To specify the axis of rotation in the toggling frame however, we should conjugate by all of the preceding pulses; this corresponds to following the rotations into the toggling-frame coordinate system, and asking which axis the rotation is being applied around in this coordinate system. Since we know the initial and final frame directions in the toggling-frame coordinate system, the rotation effect is uniquely specified and can be easily characterized by the chirality of the rotation connecting the initial and final frames.
With these considerations, we have explicitly shown that cancellation of rotation angle errors is well-captured by the chirality of each frame change, given by β µ,k . Note that this applies to both a uniform rotation angle error as well as cases with a field inhomogeneity.

E. Analysis of Rotation Axis Errors
While we have focused on a particular set of finite pulse effects and pulse imperfections that we believe to be dominant in most experimental systems, our analysis can be readily generalized to other types of imperfections [2]. As an example, here we illustrate how to incorporate into our framework the effects of a rotation axis error, in which the applied rotation axis of e.g. anx-rotation deviates from the designated direction.
To model this effect, we treat axis errors as a modification of the rotation exp[−iπS x /2] into exp[−iπ(S x ± ζS y )/2], and correspondingly exp[−iπS y /2] into exp[−iπ(S y ∓ ζS x )/2]. We assume that the rotation axis error results in shifts in the same direction for bothx andŷ axes. This is justified when microwave control pulses are directly synthesized using an arbitrary waveform generator (AWG): Here, while thex andŷ axes are 90 degrees phase shifted with respect to each other for pulses of the same length, different pulse transients for pulses of different length (e.g. π/2 and π pulses) may lead to different relative phases between the rotation axes of each pulse type. Other types of rotation axis errors can also arise in experiments due to e.g. phase distortions of IQ modulators or nonlinear transient behaviors of microwave amplifiers. We leave the treatment of these more general errors to future work.
We rewrite the rotation axis deviation as the cross product between theẑ-axis, corresponding to the disorder direction, and the ideal rotation axis direction. For example, for the rotation exp[−iπS x /2], we obtain the imperfection contributionẑ ×x = +ŷ. Directly working in the spin operator language, we can make use of the identity which resembles the form of the cross product. Thus, the rotation axis error Hamiltonian for the pulse P k can be written as where α k = ν α ν,k e ν is the rotation axis in the lab frame. In the toggling frame, this term becomes where we used Eq. (S11) and Eq. (S16). Utilizing this expression, the average Hamiltonian of the rotation axis error during the k-th pulse is calculated as where F µ (θ) = cos θF µ,k + sin θF µ,k+1 describes the continuous rotation of the frame direction during the pulse. To further simplify the expression, we make use of the fact that β k specifies a rotation from the frame F µ,k into F µ,k+1 , and can thus be written as This allows us to further simplify the expression as where in the third line, we used the identity in the fourth line, we performed the summation over λ and σ to eliminate the Kronecker deltas, and in the fifth line, we used the fact that at times k and k + 1, the frame orientation is different, and thus F µ,k · F µ,k+1 = 0 for all µ, as well as the fact that µ F 2 µ,k = 1. Thus, the condition for rotation axis deviations to be cancelled is that considering all π/2 pulses (i.e. neglecting π pulses), for each axis directionx,ŷ,ẑ, the sum of parities of the frames before each pulse is equal to the sum of parities of the frames after each pulse. Note that if we had included π pulses as well, then the sum would fully cancel and become 0, which is consistent with the fact that a global phase rotation of all pulses does not affect the decoupling performance, and can be counteracted by a change of axis.
Incorporating a rotation axis deviation between thex andŷ rotations requires modifications to the formalism. This is because the average Hamiltonian rotation axis error term will now depend on the specific pulse physically applied (x orŷ), which cannot be directly read out in our formalism using only the information of the free evolution frames immediately preceding and following the pulse. One possibility would be to keep track of the full configuration of frame orientations or make use of the history of rotations, as described in Sec. S1 G, though the analysis would be more involved. However, we emphasize that in the case where bothx andŷ rotations are directly synthesized using an AWG, we expect this type of rotation axis deviation to be very small.

F. AC Sensing Effective Field Strength
As discussed in the main text, the effective magnetic field strength as seen by the driven spins is SinceF µ (f ) 2 has the complex phase dependence e 2i(φ−φµ) , for a given target field frequency, the preceding field strength has maximum value In Eq. (S28), the equality is satisfied only when the sensing sequence is phase-matched with the target signal at the frequency f for all axes, i.e., φ =φ x (f ) =φ y (f ) =φ z (f ).

G. Obtaining S x Transformation Trajectories
In the main text, we have assumed that the periodic pulse sequence gives rise to an identity evolution operator over one Floquet cycle. However, this need not be the case, since our framework only specifies the transformations of the S z operator to be periodic, but the transverse S x,y operators may not return to their initial configurations. This may be particularly problematic for sensing purposes, since the phase accumulation in thex orŷ directions would then cancel over multiple cycles due to the changing frame directions ofx andŷ axes. In such situations, one should repeat the frame matrix multiple times until the total unitary rotation by the pulses is identity and use the enlarged frame matrix to define the Floquet period.
Here, we provide a simple method to obtain the evolution of the S x operator from the frame matrix F specifying the transformations of S z , such that one can easily keep track of the total unitary rotation up to any point. This is not only useful to determine whether the total unitary rotation is identity, but also helps to directly write down the rotation pulses that need to be applied in the lab frame to realize the transformations specified by F.
To intuitively understand how the criteria are obtained, let us consider a generic frame orientation in the free evolution period k, with the spin operators transformed into S z →S z k and S x →S x k respectively. We now consider the Bloch sphere intuition developed in Sec. S1 D. IfS z k+1 = ±S x k , this implies that the rotation P k should be applied along theŷ axis to bring the toggling frame operator pointing along thex direction into theẑ direction. Correspondingly, it will also rotate theẑ direction into thex direction, with an additional negative sign. Thus, wheneverS z k+1 = ±S x k , thex andẑ operators will switch in the subsequent frame, with the product of the signs of S x k ,S x k+1 ,S z k ,S z k+1 being −1. On the other hand, ifS z k+1 = ±S x k , then the rotation P k should be applied along thê x axis, andS x k+1 =S x k . With this information, we can sequentially construct the matrix G that corresponds to a similar representation as F, but for the transformed operatorS x . The rules are: 1. The first column is G 1 =   1 0 0   .
2. For each subsequent column k + 1, compare the location of the nonzero element in F k+1 to G k . If they are different, copy G k into G k+1 ; otherwise, copy F k into G k+1 , and multiply −1 if the nonzero elements of F k+1 and G k have the same sign.
Combining these techniques with the other rules in Tab. 1 for robust Hamiltonian engineering, we can engineer a wide range of Hamiltonians. Importantly, the disorder and interaction Hamiltonians [Eqs. (S31,S32)] can be independently engineered due to their different functional dependence on F µ,k (linear for H dis avg and quadratic for H int avg ). In addition, one can also engineer uniform single-body terms via intentional detuning or systematic rotation angle deviations (rule 4), and disordered XY-type interactions via incomplete cancellation of cross terms (rule 3). By utilizing higher-order Magnus expansions, one can also engineer effective three-body interactions generated by the commutators. This allows one to obtain the desired many-body Hamiltonian with tunable disorder and interaction type. As discussed in the main text, taking the case of NV centers as an example where J S ij = −J I ij and J A ij = 0, we can engineer Ising (c = 0), Heisenberg (c = 1/3), XY (c = 1/2), and dipolar-like (c = 1) interactions. These can exhibit very different thermalization behavior, as the Ising interaction limit is expected to be integrable, while the exchange interactions can facilitate spreading of excitations.
There are, however, a few constraints on the form of the engineered Hamiltonian: First, the Heisenberg component, S i · S j , of the interaction will remain invariant, since it is not transformed under global rotations; Second, some of the Hamiltonian terms transform in the same way (for instance disorder and anti-symmetric spin exchange), and consequently cannot be independently engineered; Finally, the resultant Hamiltonian may experience a rescaling in magnitude, due to the finite projection of the initial Hamiltonian onto the final Hamiltonian. Eventually, contributions from higher-order terms in the Magnus expansion may also become important.

I. Proof of Optimal Sensitivity Under Interaction Decoupling Constraints
To theoretically derive the maximum AC sensitivity under the constraint of interaction decoupling (where the evolution time along each axis is required to be equal), the sinusoidal amplitude modulation of the target sensing signal needs to be carefully taken into account. The effective field strengths along each axis at the resonance are proportional to the weighted integrals of the sinusoidal target signal (see Eq. (26) of the main text), and the sum of amplitudes is constrained by the identity where B XY-8 is the effective field strength obtained from the conventional XY-8 sequence. At the same time, the best AC sensitivity is achieved when the total effective field strength B eff = B 2 x + B 2 y + B 2 z is maximized. To optimize | B eff |, we first assume that the inequality Eq. (S36) is saturated, in which case |B x | = B XY-8 − |B y | − |B z |. Without loss of generality, let us also assume B x > B y > B z > 0. Taking the partial derivative, we obtain Since B z < B y < B XY-8 , we see that it is desirable to minimize the smallest element. We thus hold B z fixed at its minimal value, and optimize the values of B x and B y , in which case Similarly, the effective field strength is maximized when B y is minimized under the constraint that B z is minimal. Given this understanding, optimal sensitivity will be achieved when the phase accumulation under the sinusoidal target signal is distributed along the three axes to maximize their difference. This implies that it is desirable to distribute all phase accumulation along thex-axis to be near the anti-nodes of the sinusoid and phase accumulation along theẑ-axis to be near the nodes, as illustrated in Seq. I in Fig. 8 of the main text. Using Eqs. (26,27) presented in the main text, we estimate the effective sensing field of Seq. I to be B eff = B XY-8 1 − in the perfect pulse limit with zero pulse duration. The corresponding sensing field magnitude | B eff | ≈ 0.634B XY-8 is the theoretical upper limit of the effective field strength under the constraint of interaction decoupling. For realistic cases of finite pulse duration, we can also use the same strategy to optimize the effective field strength, by maximizing the phase accumulation imbalance between different axes. In such cases, we can choose to reduce the number of short intermediate frames along theẑ axis, and instead satisfy the interaction decoupling requirements by increasing the duration of free evolution time along theẑ axis. This will further increase the difference between the phase accumulation along different axes, leading to high sensitivity for the case of finite pulse duration. S10

S2. EXPERIMENTAL DETAILS
A. Experimental Setup Our experimental system consists of a dense electronic spin ensemble (∼15 ppm) of nitrogen-vacancy (NV) centers in diamond. We utilize a static external magnetic field to split the spin-1 ground states |0, ±1 and resonantly address the |0 and |−1 levels, effectively isolating an ensemble of two-level systems consisting of NVs with the same crystallographic axis orientation. This subgroup of NVs (∼4 ppm) is characterized to have large on-site disorder (W ∼ (2π) 4 MHz) originating from paramagnetic impurities and inhomogeneous strain, with modest interaction strengths through long-range magnetic dipolar interactions (J ∼ (2π) 35 kHz within a group of NVs with the same axis orientation). The intrinsic magnetic dipolar interaction strength is significantly larger than extrinsic decoherence rates due to coupling to the environment. Indeed, we have verified that the coherence time of the standard XY-8 sequence is limited by NV-NV interactions [3][4][5]. We utilize optical excitation and readout to initialize and detect the spin states and use resonant microwave driving to manipulate the spin state. The microwave pulses are directly synthesized using a high sampling rate AWG (Tektronix 7122B), amplified through a microwave amplifier (Mini Circuits ZHL-16W-43-S+), and delivered to the NV ensemble with an Ω-shaped coplanar waveguide. We cut out a nanobeam from a bulk sample using an angled etching procedure [6], to improve optical and microwave control homogeneity. For the sensing measurements, an additional radio-frequency signal is combined with the microwave pulses through a diplexer (Gowanda D3005001). See Refs. [3,4,7] for a more complete discussion of the experimental system and setup.

B. Data Analysis
We now discuss the detailed forms of fitting performed for the experimental data. For the spin coherence decay measurements in Fig. 7(b), we fix the free evolution time τ = 25 ns and π/2 pulse length t p = 10 ns, and sweep the number of times the sequence is repeated. We fit the experimental data to a stretched exponential profile y(t) = y 0 exp[−(t/T 2 ) α ], with fitting parameters being the 1/e coherence time T 2 , the exponent α and the maximum contrast y 0 . For ease of comparison, we plot the normalized decay curves defined as y(t)/y 0 .
In the rotation error robustness measurements in Fig. 7(c) of the main text, we fit the experimental data to a stretched exponential with a finite contrast sinusoidal modulation: where A and f characterize the contrast modulation depth and frequency, respectively. In Fig. 7(d), we extract f as a function of the Rabi frequency Ω, and we introduce the systematic rotation error = (Ω − Ω 0 )t p , where Ω 0 = π/2t p is the precalibrated Rabi frequency. Taking the mean of the extracted f values averaged over different initial states, we fit the profile of f to the functional form f (r) ∝ |r − r 0 |, where r = (Ω − Ω 0 )/Ω 0 and r 0 is a small lateral offset. In Fig. 7(e), we study the experimental contrast as a function of AC sensing signal frequency. We calculate the theoretical ESR response of the Floquet sequence based on the expected effective sensing field orientation in the toggling frame. Here, for a given pulse sequence, we first obtain the time-domain modulation functions F µ (t), taking into account the finite pulse durations. We then use Eqs. (26,27) of the main text to evaluate the effective magnetic field B eff for a given phase φ of the target AC sensing signal. For small sensing signal phase accumulation and for an initialization and readout direction alongx, the resulting contrast C due to the tilted precession can be written as where θ = γB t t/ is the rotation angle due to the AC sensing signal and B t is the total field strength. Experimentally, we optimize the phase φ of the target AC sensing signal to maximize contrast; correspondingly, we should choose φ to maximize the transverse magnetic field strength B ⊥ = B 2 y + B 2 z to capture the expected contrast theoretically [7]. We then calculate the expected contrast from Eq. (S40) using this optimal phase φ to obtain the theoretical curves shown in Fig. 7(e) for different spin initialization directions.
For Fig. 7(f), we study the contrast as a function of AC sensing signal phase φ for different initial states. The response can be calculated using the frequency-domain modulation function and effective magnetic field, as above, for each of the initialization directions. When the spectral modulation functions are phase-synchronized along all three axes, the oscillations in contrast C with respect to φ are aligned, and choosing the optimal initialization direction (perpendicular to the effective field direction) yields larger contrast (Fig. 7(f), top panel, resonance 1 ). The slight