Robust Dynamic Hamiltonian Engineering of Many-Body Spin Systems

TL;DR

This paper presents a systematic method for designing robust control pulse sequences to engineer desired Hamiltonians in many-body quantum systems, improving coherence, sensing, and simulation capabilities.

Contribution

It introduces a matrix-based algebraic framework for creating robust, efficient pulse sequences applicable to complex many-body Hamiltonians, with experimental validation.

Findings

Designed pulse sequences enhance quantum coherence and robustness.

Framework accelerates numerical search for optimal sequences.

Experimental demonstration confirms robustness in general interaction systems.

Abstract

We introduce a new approach 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 an efficient framework 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, (iv) provide optimal sequence length within given constraints, and (v) substantially accelerate numerical searches of pulse sequences. 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 system with the most general interaction form.

Figures

Equations

$$H_{\text{avg}}=\frac{1}{T}\sum _{k=1}^n{\tau }_k{\tilde{H}}_k,$$

$$H_s$$

$$=\sum _i{h}_i{S}_i^z+\sum _{ij}[J_{ij}^I{S}_i^z{S}_j^z+J_{ij}^S\text{quantity}(S_i^x{S}_j^x+S_i^y{S}_j^y)$$

$$+J_{ij}^A\text{quantity}(S_i^x{S}_j^y-S_i^y{S}_j^x)],$$

$${\tilde{S}}^z(t)$$

$$=\sum _{\mu }{F}_{\mu ,k}{S}^{\mu },\text{for tk−1<t<tk}.$$

$$F_{\mu ,k}=2\text{Tr}[S^{\mu }{\tilde{S}}_k^z],$$

$${\left(\begin{array}{c}\mathbf{F}\\ \mathbf{\tau }\end{array}\right)}_{\text{CPMG}}=\left(\begin{array}{cc}0& 0\\ 0& 0\\ +1& -1\\ \tau & \tau \end{array}\right),$$

$${\left(\begin{array}{c}\mathbf{F}\\ \mathbf{\tau }\end{array}\right)}_{\text{CPMG}}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& +1& 0& -1\\ +1& 0& -1& 0\\ \tau & 0& \tau & 0\end{array}\right),$$

$$K_{\mu }$$

$$L_{\mu }$$

$$S_i^z{S}_j^z$$

$$S_i^x{S}_j^x+S_i^y{S}_j^y$$

$$S_i^x{S}_j^y-S_i^y{S}_j^x$$

$$H_{\text{avg}}^{\text{dis}}$$

$$H_{\text{avg}}^{\text{int}}$$

$$+\sum _{ij,\mu }{J}_{ij}^S{S}_i^{\mu }{S}_j^{\mu }\cdot \left(1-L_{\mu }\right)$$

$$+\sum _{ij,\mu }{J}_{ij}^A({\vec{S}}_i\times {\vec{S}}_j{)}^{\mu }\cdot K_{\mu }.$$

$$\sum _{k=1}^n{F}_{\mu ,k}{\tau }_k=0,\text{for every μ=1,2,3}.$$

$$\sum _{k=1}^n|F_{x,k}|{\tau }_k=\sum _{k=1}^n|F_{y,k}|{\tau }_k=\sum _{k=1}^n|F_{z,k}|{\tau }_k.$$

$$\delta H_{\text{avg}}=\frac{t_p}{T}\sum _{k=1}^n{\bar{H}}_{P_k}^{(0)},$$

$$H_{\text{eff}}=H_{\text{avg}}+\delta H_{\text{avg}}.$$

$${\bar{H}}_{P_k}^{(0)}=\frac{1}{t_p}{\int }_0^{t_p}(P_{k-1}\cdots P_1{)}^{†}{\tilde{H}}_{P_k}(t)(P_{k-1}\cdots P_1)dt,$$

$${\tilde{H}}_{P_k}(t)$$

$${\bar{H}}_{P_k}^{(0)}=$$

$$+\left[\left(\frac{{\tilde{H}}_k^I+{\tilde{H}}_{k+1}^I}{2}\right)+\left(\frac{{\tilde{H}}_k^S+{\tilde{H}}_{k+1}^S}{2}\right)+{\tilde{H}}_{k,k+1}^C\right],$$

$${\tilde{H}}_{k,k+1}^{\text{C}}$$

$$=\sum _{ij,\mu \nu }{C}_{ij}^{\mu \nu }\sum _{k=1}^n{\mathcal{P}}_k^{\mu \nu },$$

$${\mathcal{P}}_k^{\mu \nu }=F_{\mu ,k}{F}_{\nu ,k+1}+F_{\nu ,k}{F}_{\mu ,k+1}.$$

$$H_{\text{eff}}=\frac{1}{T}$$

$$+\sum _{k=1}^n\text{quantity}({\tau }_k+t_p)\text{quantity}({\tilde{H}}_k^{\text{I}}+{\tilde{H}}_k^{\text{S}})+\sum _{k=1}^n{t}_p{H}_{k,k+1}^{\text{C}}],$$

$${\vec{\beta }}_k={\vec{F}}_{k+1}\times {\vec{F}}_k,$$

$$\delta H_{\text{avg}}^{\text{rot}}=\frac{1}{T}\sum _{i,\mu }{ϵ}_i\text{quantity}(\sum _{k=1}^n{\beta }_{\mu ,k}{S}_i^{\mu }),$$

$$\sum _k{\beta }_{\mu ,k}=0\text{for each μ=x,y,z}$$

$${\bar{H}}^{(1)}$$

$${\Theta }_k$$

Full text

I Introduction and Motivation

The ability to control and manipulate the dynamics of a quantum system in a robust fashion is key to many quantum technologies. In particular, the use of periodic control pulses, also known as Floquet driving, has emerged as a ubiquitous tool for the control and engineering of quantum dynamics Haeberlen and Waugh (1968); Vandersypen and Chuang (2005); Eckardt (2017); Oka and Kitamura (2018); Bukov et al. (2015); Goldman and Dalibard (2014); Poudel et al. (2015), with applications in protecting quantum coherence from environmental noise Hahn (1950); Carr and Purcell (1954); Meiboom and Gill (1958); Gullion et al. (1990); Viola et al. (1999); Khodjasteh and Lidar (2005); Viola and Knill (2005); Uhrig (2007); Khodjasteh and Viola (2009a); Biercuk et al. (2009); Du et al. (2009); Álvarez et al. (2010); West et al. (2010); de Lange et al. (2010); Ryan et al. (2010); Khodjasteh et al. (2013); Suter and Álvarez (2016); Burum and Rhim (1979); Cory et al. (1990a); Iwamiya et al. (1993); Naydenov et al. (2011); Genov et al. (2017) and frequency-selective quantum sensing Schirhagl et al. (2014); Degen et al. (2017); Taylor et al. (2008); Degen (2008); Bylander et al. (2011); Hall et al. (2010); Naydenov et al. (2011); Álvarez and Suter (2011); de Lange et al. (2011); Pham et al. (2012); Norris et al. (2016); Frey et al. (2017); Rose et al. (2018); Fiderer and Braun (2018); Lang et al. (2015). Periodic control can also be employed to engineer the interactions between qubits, even when only global control is available, enabling the study of out-of-equilibrium phenomena, such as dynamical phase transitions and quantum chaos, and the observation of novel phases of matter such as discrete time crystals Lindner et al. (2011); Jiang et al. (2011); Heyl (2018); D’Alessio et al. (2016); Garttner et al. (2017); Khemani et al. (2016); Else et al. (2016); von Keyserlingk et al. (2016); Yao et al. (2017); Choi et al. (2017a); Zhang et al. (2017); Sacha and Zakrzewski (2018); Nandkishore and Huse (2015); Abanin et al. (2018); Álvarez et al. (2015); Wei et al. (2018a, b); Ho et al. (2017); Choi et al. (2019).

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) Slichter (2013); Mehring (2012); Levitt (2001); Sørensen et al. (1984); Rhim et al. (1971); Drobny et al. (1978); Burum and Rhim (1979); Shaka et al. (1983); Baum et al. (1985); Shaka et al. (1988); Tycko (1990); Cory et al. (1990a); Lee et al. (1995); Hohwy et al. (1999); Carravetta et al. (2000); Takegoshi and McDowell (1985); Rose et al. (2018); Lee and Goldburg (1965); Vinogradov et al. (1999); Iwamiya et al. (1993); Sakellariou et al. (2000).

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 Kucsko et al. (2018); Mohammady et al. (2018).

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 many pulse sequences that retain robustness to some of these control imperfections Rhim et al. (1974); Burum and Rhim (1979); Cory et al. (1990a), 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.

In this work, we introduce a novel framework to systematically address these challenges and efficiently design robust, self-correcting pulse sequences for dynamic Hamiltonian engineering in interacting spin ensembles using only global control Hayes et al. (2014); Ajoy and Cappellaro (2013); Choi et al. (2017b); O’Keeffe et al. (2019); Haas et al. (2019); ’Attar et al. (2019). Such globally controlled spin ensembles are naturally realized in various systems Waugh et al. (1968); Wei et al. (2018a); Kucsko et al. (2018); Tyryshkin et al. (2003); Blatt and Roos (2012); Jurcevic et al. (2014); Bohnet et al. (2016); Zhang et al. (2017); Labuhn et al. (2016); Bernien et al. (2017). We demonstrate both theoretically and experimentally that our approach has immediate applications ranging from dynamical decoupling and quantum metrology to quantum simulation.

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 $\pi/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, finite-duration 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 experimental 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. Combining optimal choices of the effective DC sensing field and initialization/readout directions with coherence time extensions through disorder and interaction suppression, this approach can lead to high sensitivity magnetometry beyond the limit imposed by spin-spin interactions, as we show in a separate manuscript Zhou et al. (2019). 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 with the most general form of interactions.

The main advances enabled by our approach include:

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. 2.

Generality: Our approach is applicable to generic two-level spin ensembles in a strong quantizing field, as typically found in most experimental quantum many-body 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 [Sec. IV.3]. 3. 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 [Sec. V]. This enables the development of pulse sequences designed for disorder-dominated systems, beyond the typical NMR setting. 4. 4.

Efficiency: Using simple algebraic conditions, we can find the shortest possible sequence length required to achieve a target Hamiltonian [Sec. V.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. VI.2]. The algebraic conditions also substantially improve numerical searches of pulse sequences by constraining the search space to a set of good pulse sequences [Sec. IV.2].

The paper is organized as follows: In Secs. II and III, we provide the theoretical framework for systematic pulse design. This is extended to higher-order and more complex, multi-body interacting Hamiltonians via analytical and numerical approaches in Sec. IV. In Secs. V, VI and VII, we present system-targeted sequence design for the applications of dynamical decoupling, quantum sensing and quantum simulation, respectively. Finally, Sec. VIII presents the experimental demonstration of our results to dynamical decoupling, with a particular focus on the broad applicability under different forms of the Hamiltonian. We conclude with a further discussion and outlook of the framework in Sec. IX.

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 interaction 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. Mansfield (1971), also closely related to control matrices Green et al. (2013); Paz-Silva and Viola (2014) and vector modulation functions Lang et al. (2017, 2019); Schwartz et al. (2018); Wang et al. (2019). Note however that the form of the Hamiltonian is limited in these existing papers, and robust decoupling rules in the interacting regime have not been derived. 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 Haeberlen and Waugh (1968) (see Appx. A for a detailed review of AHT). In particular, for a pulse sequence consisting of $n$ spin-rotation pulses {$P_{1,\cdots,n}$}, the leading-order average Hamiltonian, $H_{\text{avg}}$, is a simple weighted average of the toggling-frame Hamiltonians

[TABLE]

where $\tau_{k}$ is the pulse spacing between the $k-1$-th and $k$-th control pulses $P_{k-1}$ and $P_{k}$, $\tilde{H}_{k}=(P_{k-1}\cdots P_{1})^{\dagger}H_{s}(P_{k-1}\cdots P_{1})$ is the toggling-frame Hamiltonian that governs spin dynamics during the $k$-th evolution period, $\tau_{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 Mansfield (1971). Our representation is based on the time-domain transformations of a single-body $S^{z}$ spin operator in the interaction picture. As shown in Sec. S1A SM , the representation works for general two-level 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 $\hat{z}$, which can also be written as $[H_{s},S^{z}_{\text{tot}}]=0$, where $S^{z}_{\text{tot}}=\sum_{i}S_{i}^{z}$ is the total spin projection operator along the $\hat{z}$-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 Bernien et al. (2017); Zhang et al. (2017), spin-exchange Kucsko et al. (2018); Mohammady et al. (2018), dipolar Waugh et al. (1968), and even exotic three-body interactions Büchler et al. (2007); Mezzacapo et al. (2014); Chancellor et al. (2017).

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 SM )

[TABLE]

where the first term $H_{\text{dis}}$ is the on-site disorder Hamiltonian and the second term $H_{\text{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 $\pi/2$-rotation building blocks around the $\hat{x}$, $\hat{y}$ axes [Fig. 1(d)]. In this setting, let us consider the interaction-picture transformations of the $S^{z}$ operator: $\tilde{S}^{z}(t)=U_{c}^{\dagger}(t)S^{z}U_{c}(t)$, where $U_{c}(t)$ 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 $\pm 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

[TABLE]

Here, $P_{k}$ is the global spin rotation performed right after the $k$-th toggling frame, $t_{k}=\sum_{j=1}^{k}\tau_{j}$ with $t_{0}$ = 0, and ${\bf F}=[F_{\mu,k}]=[\vec{F}_{x};\vec{F}_{y};\vec{F}_{z}]$ is a $3\times n$ matrix containing elements 0 and $\pm 1$. The matrix elements $F_{\mu,k}$ can be explicitly calculated as

[TABLE]

with $\tilde{S}^{z}_{k}=\tilde{S}^{z}(t)$ for $t_{k-1}<t<t_{k}$. Intuitively, a nonzero element $F_{\mu,k}$ indicates that the initial $S^{z}$ operator transforms into $S^{\mu}$ for the duration of the free evolution interval $\tau_{k}$, with its sign determined by $F_{\mu,k}$. Additionally, the time duration of each toggling frame is specified by the frame-duration vector $\bm{\tau}=[\tau_{1},\tau_{2},\cdots,\tau_{n}]$.

This representation is illustrated for three pulse sequences in Fig. 2. The CPMG sequence Carr and Purcell (1954); Meiboom and Gill (1958); Gullion et al. (1990), consisting of equidistant $\pi$ pulses to suppress on-site disorder, can be represented as

[TABLE]

since the first $+S_{z}$ is flipped to $-S_{z}$ by a $\pi$ pulse [Fig. 2(a,d)]. Similarly, the WAHUHA sequence, consisting of four $\pi/2$ pulses Waugh et al. (1968), 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 111 In our algebraic conditions, we use the convention where each free evolution time is immediately followed by a pulsed rotation. For base pulse sequences that end with a free evolution in the original $\hat{z}$-axis without any following pulse, as in the WAHUHA sequence case [Fig. 2(b,e)], we move the final free-evolution block to the beginning of the pulse sequence representation and combine it with the first frame, before applying our algebraic conditions for robustness [Tab. 1].. 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-frame $\tilde{S}^{z}$ orientation after each instantaneous pulse. However, the rotation axis of $\pi$ pulses is not yet uniquely specified. To address this, we decompose all pulses into $\pi/2$ building blocks, and specify intermediate frames for pulses of rotation angles larger than $\pi/2$. Such a $\pi/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 $\pi$ pulse is then split into two $\pi/2$ pulses with zero time separation, where the first $\pi/2$ pulse along the $\hat{x}$-axis rotates $+S^{z}$ into the intermediate frame $+S^{y}$ and the second $\pi/2$ pulse along the same axis brings $+S^{y}$ into $-S^{z}$. For example, the CPMG sequence involving $\pi$ pulses can be represented as

[TABLE]

which now unambiguously specifies the rotation axes for the $\pi$ pulses. Note that the sequence contains zeros in its frame-duration 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 $\pi/2$ pulse structures, which will play an important role in robust quantum sensing sequences [Sec. VI].

One key advantage of our pulse sequence representation $\{{\bf F},\bm{\tau}\}$ is that we can now conveniently obtain the engineered $H_{\text{avg}}$ [Eq. (1)] from any many-body Hamiltonian of the form Eq. (2). More specifically, we find that the weighted row-sums and weighted absolute row-sums (which are equivalent to row square sums, since each element $F_{\mu,k}$ takes on values $\{0,\pm 1\}$) of the sequence matrix,

[TABLE]

fully specify the generalized formula for the average Hamiltonian $H_{\text{avg}}$ as a result of the following toggling-frame transformations of two-body Ising, symmetric exchange and anti-symmetric exchange interaction Hamiltonians:

[TABLE]

These expressions can be intuitively understood by examining how the $S^{z}$ operator is transformed, and using an analogy between the anti-symmetric interaction form and cross products (see Appx. B for detailed derivations). Using $K_{\mu},L_{\mu}$ defined above, we can thus write the leading-order average Hamiltonian, $H_{\text{avg}}=H_{\text{avg}}^{\text{dis}}+H_{\text{avg}}^{\text{int}}$, as

[TABLE]

II.2 Decoupling Conditions for Ideal Pulses

Our goal here is to perform dynamical decoupling and suppress both disorder and interaction effects, by generating a pulse sequence with a vanishing $H_{\text{avg}}=0$ Slichter (2013); Mehring (2012); Levitt (2001); Rhim et al. (1971); Burum and Rhim (1979); Cory et al. (1990a); Rose et al. (2018). Examining the above expressions in Eqs. (10-13), we observe that there are two types of functional dependencies on $F_{\mu,k}$: the disorder [Eq. (10)] and anti-symmetric spin-exchange [Eq. (13)] Hamiltonians involve terms linear in $F_{\mu,k}$, while the Ising [Eq. (11)] and symmetric spin-exchange [Eq. (12)] Hamiltonians involve terms quadratic in $F_{\mu,k}$.

The first type of contribution, which depends linearly on $F_{\mu,k}$, can be cancelled if $K_{\mu}$ = 0 for all $\mu=x,y,z$ axes [see Eqs. (10, 13)], giving

[TABLE]

For equidistant pulse sequences where $\bm{\tau}=\tau{\bf I}_{1\times n}$, the above condition can be further simplified to $\sum_{k=1}^{n}F_{\mu,k}=0$. This suggests that each row $\mu$ of the matrix ${\bf F}$ should have an equal number of positive and negative elements, such that their sum is 0, resulting in $H_{\text{avg}}^{\text{dis}}=H_{\text{avg}}^{\text{A}}=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_{\text{avg}}^{\text{dis}}=\sum_{i}\frac{h_{i}}{3}(S_{i}^{x}+S_{i}^{y}+S_{i}^{z})$ [Fig. 2(b,e)].

In contrast, the terms with quadratic dependence on $F_{\mu,k}$ cannot be fully suppressed in general, as the isotropic (Heisenberg) component of the interaction is invariant under global rotations Choi et al. (2017b), leading to $H_{\text{avg}}^{\text{int}}\neq 0$. However, it is still possible to fully symmetrize these interactions into a Heisenberg Hamiltonian, $H_{\text{avg}}^{\text{int}}=\frac{1}{3}\sum_{ij}(J_{ij}^{I}+2J_{ij}^{S})\vec{S}_{i}\cdot\vec{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}$, giving

[TABLE]

Again, for equidistant pulses, this condition is simplified to the statement that the sum $\sum_{k=1}^{n}|F_{\mu,k}|$ should be the same for each $\mu$. Based on this analysis, we verify that the CPMG sequence [Fig. 2(a,d)] does not symmetrize interactions, since it only employs $\pi$ pulses, while the two sequences that incorporate $\pi/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_{\text{avg}}^{\text{int}}=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 Sec. IV.3 and Appx. C, we utilize results from unitary $t$-designs Dankert et al. (2009); Webb (2015); Zhu (2017) 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 Rhim et al. (1974); Burum and Rhim (1979); Cory et al. (1990a). 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 $\delta H_{\text{avg}}$, which can be estimated to leading order using AHT as

[TABLE]

where $t_{p}$ is the duration of a $\pi/2$ pulse and $\bar{H}^{(0)}_{P_{k}}$ is the zeroth-order average Hamiltonian acting during the $k$-th $\pi/2$ pulse building block. Thus, the total leading-order effective Hamiltonian describing the driven spin dynamics is given by

[TABLE]

The goal of robust Hamiltonian engineering is to suppress the error $\delta H_{\text{avg}}=0$ by designing leading-order fault-tolerant, self-correcting 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 ${\bf F}$ directly provides a simple way to obtain $\delta H_{\text{avg}}$ [Eq. (16)]. Intuitively, the form of $\bar{H}_{P_{k}}^{(0)}$ in $\delta H_{\text{avg}}$ is expected to be the average of the neighboring toggling-frame Hamiltonians, $\tilde{H}_{k}$ and $\tilde{H}_{k+1}$, since the finite-duration pulse $P_{k}$ smoothly changes the spin frame from the $k$-th to the $k+1$-th 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 $\pi/2$-pulse building blocks, we can analytically calculate $\bar{H}_{P_{k}}^{(0)}$ originating from the finite-duration pulse $P_{k}$ as

[TABLE]

where $\{P_{1},\cdots,P_{k-1}\}$ denotes the preceding $k-1$ rotations and the effect of the pulse $P_{k}$ is given by

[TABLE]

Here, $Q_{k}(t)=\exp\quantity[-i\sum_{i}\Omega_{i}tS_{i}^{\nu}]$ is the time-dependent unitary operator due to the pulse $P_{k}$ that globally rotates spins along the $\nu$-axis over the finite duration $0\leq t\leq t_{p}$ and $\Omega_{i}$ is the Rabi frequency for a spin at site $i$, producing the $\pi/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 $\pi/2$ pulses is to smoothly interpolate the toggling-frame spin operator $\tilde{S}^{z}(\theta)=\cos\theta\tilde{S}_{k}^{z}+\sin\theta\tilde{S}_{k+1}^{z}$ during the finite pulse duration, where $\theta=\Omega t$. Using this to evaluate the integral of Eq. (18), as detailed in Sec. S1B SM , we obtain

[TABLE]

where $\tilde{H}_{k}^{\text{dis}}$, $\tilde{H}_{k}^{\text{I}}$, $\tilde{H}_{k}^{\text{S}}$ and $\tilde{H}_{k}^{\text{A}}$ 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^{\mu}$ in $H_{s}$ [Eq. (2)] to the toggling-frame ones $\tilde{S}^{\mu}_{k}$, as shown in the individual terms in the summations of Eqs. (10)-(13). While most terms in Eq. (20) are the weighted average of the neighboring toggling-frame Hamiltonians, consistent with the original intuition, there is an additional average Hamiltonian contribution $\tilde{H}_{k,k+1}^{\text{C}}$ resulting from two-body interaction cross-terms $\tilde{S}^{z}_{k}\tilde{S}^{z}_{k+1}$ acting during the pulse and given by

[TABLE]

where $C_{ij}^{\mu\nu}=\sum_{ij}J^{C}_{ij}(S^{\mu}_{i}S^{\nu}_{j}+S^{\nu}_{i}S^{\mu}_{j})$ is the cross-interaction operator with $J_{ij}^{C}=\frac{1}{\pi}(J^{I}_{ij}-J^{S}_{ij})$ (see Sec. S1B SM ) and $\mathcal{P}_{k}^{\mu\nu}$ defines the “parity” of neighboring $k$-th and $k+1$-th frames, given as

[TABLE]

To cancel the interaction cross-terms, the parity should vanish when summed over one Floquet period for each pair ($\mu,\nu$): $\sum_{k=1}^{n}\mathcal{P}_{k}^{\mu\nu}=0$ [Fig. 3(e)]. Intuitively, the parity can be understood as checking whether the signs of neighboring frames are the same ($\mathcal{P}_{k}^{\mu\nu}=+1$, even parity) or different ($\mathcal{P}_{k}^{\mu\nu}=-1$, odd parity).

Taking into account the distinct weighting factors of $4/\pi$ and 1 for the different interaction types as well as the additional interaction cross-term, as identified in Eq. (20), the effective Hamiltonian in the presence of finite pulse duration $t_{p}$, $H_{\text{eff}}=H_{\text{avg}}+\delta H_{\text{avg}}$, becomes

[TABLE]

where the base pulse sequence length, $T=(\sum_{k=1}^{n}\tau_{k})+nt_{p}$, includes the total length of $n$ $\pi/2$ pulses. As described above in the discussion following Eq. (20), each of the terms in the toggling-frame Hamiltonian, $\tilde{H}^{\text{dis}}$, $\tilde{H}^{\text{A}}$, $\tilde{H}^{\text{I}}$, $\tilde{H}^{\text{S}}$, and $\tilde{H}^{\text{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 $+\hat{\mu}$ direction (positive chirality) in the toggling frame can be compensated by another rotation around the $-\hat{\mu}$ direction (negative chirality). Moreover, the rotation axis in the toggling frame can be readily described using our frame matrix $\mathbf{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, $\vec{\beta}_{k}=\sum_{\mu}\beta_{\mu,k}\hat{e}_{\mu}$, can be obtained by taking the cross product of the frame vectors before and after the pulse

[TABLE]

where $\vec{F}_{k}=\sum_{\mu}F_{\mu,k}\hat{e}_{\mu}$. Physically, this cross product structure can be thought of as characterizing the chirality of the toggling-frame rotation from $\vec{F}_{k}$ to $\vec{F}_{k+1}$. As derived in more detail in Sec. S1C SM , the average Hamiltonian contribution corresponding to the rotation angle error is then given by

[TABLE]

where $\epsilon_{i}$ is the static rotation-angle deviation from the target $\pi/2$ angle for a spin at site $i$. This allows us to identify the cancellation condition for rotation angle errors, $\delta H_{\text{avg}}^{\text{rot}}=0$, as

[TABLE]

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. (14,15), 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 ${\bf 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_{\mu,k}$ = $+1$) and green ($F_{\mu,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 Sec. IV.3 and Appx. C, 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 Rhim et al. (1974) (see for example Sec. S1D SM for the treatment of rotation axis imperfections).

IV Extensions to Higher-order Average Hamiltonians and Multi-Body Interactions

IV.1 Suppression of Higher Order Effects

While our analysis thus far has focused on the zeroth-order average Hamiltonian, we can incorporate higher-order effects by considering the full Magnus expansion [Eq. (47)] to engineer effective Hamiltonians with higher accuracy. More specifically, our frame matrix representation allows us to readily evaluate the higher-order expansion terms, which consist of commutators between Hamiltonians at different times in the toggling frame (see Appx. A). For example, the first-order contribution [Eq. (49)] for a periodic pulse sequence, including finite pulse effects, can be expressed as

[TABLE]

where $n$ is the total number of evolution intervals and $\Theta_{l,k}$ are the time-weighted Hamiltonians in the $l$-th and $k$-th toggling frames, respectively, given as (including finite-pulse duration effects)

[TABLE]

and the $O(t_{p}^{2})$ term coming from commutations of a finite pulse Hamiltonian with itself will typically be small (there are $n$ such terms, compared to $n^{2}/2$ terms for the $\Theta$ terms, and $t_{p}$ is typically small). Recall from Eq. (23) that the zeroth-order Hamiltonian is simply $\bar{H}^{(0)}=\frac{1}{T}\sum_{k=1}^{n}\Theta_{k}$. Crucially, note that all $\Theta_{k=\{1,\cdots,n\}}$ are easily numerically computable with our matrix representation that specifies the toggling-frame evolution of the $S^{z}$ operator [Eqs. (10-13)]. This enables us to readily evaluate the contribution from the first-order term [Eq. (27)], which will result in algebraic conditions that involve second-order polynomials in $F_{\mu,k}$ and $|F_{\mu,k}|$. Analyzing the second-order or even higher-order terms is also straightforward as they can be obtained in a very similar fashion via recursive computation of nested commutators involving $\Theta_{k}$’s at different frames.

As an explicit example, the first-order term for the echo+WAHUHA sequence [Fig. 2(c,f)] can be analytically derived to be (assuming uniform $\tau_{k}$)

[TABLE]

where we have simplified the expression by dropping anti-symmetric exchange interactions that are typically not present, and assuming $t_{p}\ll\tau_{k}$.

In addition to explicitly evaluating the higher-order expansion terms, one can also use various heuristics to suppress higher-order terms and enhance the accuracy of Hamiltonian engineering. Developed primarily in the NMR community, there are several known approaches, such as reflection-symmetric pulse arrangements Mansfield (1971); Burum and Rhim (1979); Cory et al. (1990a); Li et al. (2007) and concatenated sequence symmetrization Khodjasteh and Lidar (2005); Souza et al. (2011); Wang et al. (2012); Farfurnik et al. (2015), 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 Mansfield (1971), that is, $\tilde{H}_{n+1-k}=\tilde{H}_{k}$, all odd-order terms in the Magnus expansion vanish, $\bar{H}^{(2l-1)}$ = 0 with integer $l$ [Eq. (47)]. In our framework, this imposes an additional condition on the sequence ${\bf 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. As an example, we can apply the reflection symmetry to the echo+WAHUHA sequence [Fig. 2(c,f)] to cancel higher-order effects, as shown in Fig. 4(a,b). 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. VI].

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 Burum and Rhim (1979) and concatenated symmetrization schemes Khodjasteh and Lidar (2005); Souza et al. (2011); Wang et al. (2012); Farfurnik et al. (2015), where a long pulse sequence is constructed from the repetition of short pulse sequences, symmetrized in a systematic pattern to suppress higher-order effects. Our method can facilitate the robust implementation of such concatenation schemes by providing both an intuitive visualization and precise algebraic conditions to analyze the error robustness of concatenated pulse sequences.

Second averaging: The technique of second-averaging has been developed in NMR to suppress dominant error terms in $\delta H_{\text{avg}}$ that do not commute with the leading contribution in $H_{\text{avg}}$ Haeberlen et al. (1971); Cory et al. (1990b); Cory (1996). Such methods can be readily incorporated in our framework by alternating the rotation axes of control pulses periodically every Floquet cycle, or by using off-resonant driving.

IV.2 Enhanced Numerical Search of Pulse Sequences

Our formalism not only enables efficient pen-and-paper pulse sequence design and provides important analytical insights, but can also greatly enhance the numerical search of pulse sequences. More specifically, the concise decoupling rules we have derived above provide a rapid means to narrow the search space down to pulse sequences that may have good performance, as a starting point for in-depth numerical simulations that capture the full dynamics of the system to all orders.

To illustrate this with a concrete example, we consider pulse sequences with 12 free evolution intervals that aim to efficiently decouple the effects of interactions and disorder. An exhaustive search of just such pulse sequences ignoring finite pulse duration effects would already require an enumeration of $6^{12}\approx 10^{9}$ possibilities (6 possible configurations for each toggling frame, i.e., $F_{\mu,k}=\pm 1$ for $\mu=x,y,z$), a prohibitively large number for numerical simulations. However, the application of our disorder- and interaction-decoupling rules on the generation of sequences can significantly narrow down the search space. In addition, depending on the target application, more constraints in the form of algebraic rules can be simultaneously applied to further reduce the size of the sequence space. For example, for efficient AC-field sensing we can impose a fast spin-echo structure whereby the signs of toggling frames are periodically flipped over the shortest possible period 2$\tau$ while maintaining a synchronized phase relation between different axes (see Fig. 6(b) as an example). Such a phase-locked, fast-echo structure acts as a bandwidth filter centered at a target frequency (see Sec. VI for detailed discussions). This allows us to find a total of 14,080 sequences that can be sorted out into four different categories according to their error robustness: Class I satisfies all decoupling rules (448 sequences), Class II does not fully decouple interaction cross-terms [violation of Condition 3 in Tab. 1], Class III does not fully suppress rotation angle errors [violation of Condition 4 in Tab. 1], and Class IV does not suppress interaction cross terms and rotation angle errors [violation of both Conditions 3, 4 in Tab. 1].

To evaluate the performance of the sequences, we numerically solve the exact Floquet spin dynamics for a disordered, interacting 8-spin system and monitor global spin polarization as a function of time.We choose Gaussian random on-site disorder ($\sigma_{W}=(2\pi)4.0$ MHz) and uniform random interactions ($J=(2\pi)0.2$ MHz), with pulse spacing $\tau=25$ ns and pulse duration $t_{p}=10$ ns, and extract the $1/e$ coherence decay times $T_{2}$ averaged over $\hat{x},\hat{y},\hat{z}$ initial states (with the average performed over decay rates). As shown in Fig. 4(c), Class I pulse sequences, satisfying all decoupling rules, perform considerably better than the other classes. In particular, we find that the top 10 sequences with longest coherence times all consistently belong to Class I. Note that the numerically-optimized sequences exhibit a broad distribution of coherence times due to different amounts of contributions from higher-order terms in the Magnus expansion. Indeed, using Eq. (27), we explicitly verify in Sec. S1F SM that the resulting coherence decay is strongly correlated with the first order contribution. These results further confirm that the analytical insights provided by our formalism can substantially improve the numerical search efficiency for optimal pulse sequences, allowing fast numerical optimization that can capture effects from all orders.

IV.3 Extensions to Multi-Body Interactions

Our discussion thus far has focused on the case of one- and two-body interactions. Interestingly, our versatile formalism can also be applied to 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. In particular, we show via a neat connection to unitary $t$-designs DiVincenzo et al. (2002); Dür et al. (2005); Emerson et al. (2003); Collins and Śniady (2006); Gross et al. (2007); Ambainis and Emerson (2007); Dankert et al. (2009); Webb (2015); Zhu (2017) that (i) 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, and (ii) we can extend the formalism that accounts for finite pulse duration effects to the case of three-body interactions, leading to new decoupling conditions beyond those discussed in Sec. III. Such interactions are important building blocks of exotic topological phenomena Moore and Read (1991); Fradkin et al. (1998); Moessner and Sondhi (2001); Levin and Wen (2005), and have been proposed to be realized in cold molecules Büchler et al. (2007) and superconducting qubits Mezzacapo et al. (2014); Chancellor et al. (2017). We sketch the main ideas of the derivation here, and detailed proofs can be found in Appendix. C.

First, let us consider the case of perfect, infinitely short pulses. We will show, via connections to unitary $t$ designs, that under the above decoupling conditions, 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

[TABLE]

Here, $\mathcal{U}(2)$ is the unitary group of dimension 2, used to describe two-level systems, $\mathcal{O}$ is an $N$-body operator with $N\leq t$ and $\mathcal{O}_{U}$ is the corresponding averaged observable. Intuitively, this expression 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. (30), where the average is taken over all elements of the unitary group over the Haar measure, imply that $\mathcal{O}_{U}$ must only contain terms proportional to elements of the symmetric group $S_{t}$ of order $t$ Collins and Śniady (2006). 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^{\otimes N}$ conjugates all spins identically.

It is known that the Clifford group forms a unitary 3-design Webb (2015); Zhu (2017). Combined with the fact that for interactions under the secular approximation, averaging over the Clifford group is equivalent to averaging over the six axis directions (see Appx. C), this implies that for any sequence that satisfies the above decoupling rules, 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_{\text{int}}=J\sum_{ijk}(S_{i}^{x}S_{j}^{y}S_{k}^{z}-S_{i}^{y}S_{j}^{x}S_{k}^{z})$. The symmetrized Hamiltonian can be calculated to be $\overline{H}_{\text{int}}=\frac{J}{3}\sum_{ijk}\sum_{\mu\nu\sigma}\epsilon_{\mu\nu\sigma}S_{i}^{\mu}S_{j}^{\nu}S_{k}^{\sigma}$, where $\epsilon_{\mu\nu\sigma}$ 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.

In fact, we can also extend this analysis to the case of finite pulse durations by expanding the Hamiltonian as a polynomial in $F_{\mu,k}$ and examining how different possible terms transform. As described in Appx. C, this gives rise to new decoupling conditions in the three-body case, as a generalization of the interaction cross-term decoupling condition [Condition 3 in Tab. 1].

V Application: Dynamical Decoupling

V.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 Hahn (1950); Carr and Purcell (1954); Meiboom and Gill (1958); Gullion et al. (1990); Viola et al. (1999); Khodjasteh and Lidar (2005); Viola and Knill (2005); Uhrig (2007); Khodjasteh and Viola (2009a); Biercuk et al. (2009); Du et al. (2009); Álvarez et al. (2010); West et al. (2010); de Lange et al. (2010); Ryan et al. (2010); Suter and Álvarez (2016); Burum and Rhim (1979); Naydenov et al. (2011); Rhim et al. (1971); Drobny et al. (1978); Burum and Rhim (1979); Shaka et al. (1983); Baum et al. (1985); Shaka et al. (1988); Tycko (1990); Cory et al. (1990a); Lee et al. (1995); Hohwy et al. (1999); Carravetta et al. (2000); Slichter (2013); Mehring (2012); Levitt (2001); Sørensen et al. (1984), magic echoes Takegoshi and McDowell (1985); Rose et al. (2018), frequency- and phase-modulated continuous driving Lee and Goldburg (1965); Vinogradov et al. (1999), and numerically-optimized control schemes Iwamiya et al. (1993); Sakellariou et al. (2000).

However, these sequences are optimized for interaction-dominated 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 Kucsko et al. (2018); Choi et al. (2017c), such that a naive application of the NMR pulse sequences performs poorly (see Sec. VIII). 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 and interactions, respectively (illustrated in Fig. 5(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\gg J$), we thus require $T_{W}\ll T_{J}$ to reduce the magnitude of higher-order contributions to $H_{\text{eff}}$ Cory (1991); Burum and Rhim (1979); Cory et al. (1990a) and maximize the leading-order approximation accuracy. If the control error magnitude $\epsilon$ 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 Cory et al. (1990a) developed for nuclear spin systems, where spin-spin interactions dominate over disorder. Indeed, as shown in Fig. 5(a), Seq. A symmetrizes interactions very rapidly and also cancels on-site disorder, but on a much slower timescale ($T_{J}\ll T_{W}$). It is also robust against leading-order imperfections resulting from finite pulse durations, and suppresses certain higher-order effects Cory et al. (1990a). For comparison, we design a new sequence, Seq. B in Fig. 5(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 $\pi$ pulses to echo out disorder on a rapid timescale while symmetrizing interactions on a slower timescale ($T_{W}\ll T_{J}$). We emphasize that Seq. B incorporates both $\pi$ pulses and composite $\pi/2$ pulses when switching between toggling frames [Fig. 5(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 disorder-dominated systems (e.g. for electronic spin ensembles). From numerical simulations, as shown in Fig. 5(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(\tau+t_{p})<0.09$ and $0.009<J(\tau+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.

V.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 $\pm 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 $\tau$ connected by composite $\pi/2$ pulses, satisfies all leading-order decoupling requirements:

[TABLE]

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. DD.1 for details.

VI 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 Naydenov et al. (2011); de Lange et al. (2011)), there are only a limited number of such demonstrations for strongly-interacting systems Cory et al. (1990a); McDonald and Tokarczuk (1989), not achieving optimal AC sensitivity, despite the pressing need for such protocols to further improve sensitivity in high density spin ensembles Acosta et al. (2009); Mitchell (2019); Mitchell and Alvarez (2019); oth . 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.

VI.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

[TABLE]

where $\gamma$ is the gyromagnetic ratio of the spins, $B_{\text{AC}},f$ and $\phi$ are the amplitude, frequency, and phase of the target AC signal, respectively.

For a given pulse sequence represented as ${\bf F}=[F_{\mu,k}]$, we can apply the same average Hamiltonian analysis to understand how driven spins experience the sensing field in their effective frame, giving

[TABLE]

with

[TABLE]

for $\mu=x,y,z$, and $\real$ denotes the real part. Physically, the time-averaged sensing-field Hamiltonian [Eq. (37)] has a simple and elegant interpretation: in the toggling-frame picture, all driven spins will undergo a coherent precession around the effective magnetic field $\vec{B}_{\text{eff}}$. Additionally, as seen in Eq. (38), the orientation and strength of $\vec{B}_{\text{eff}}$ are determined by the frequency-domain resonance characteristics of the applied pulse sequence, $\tilde{F}_{\mu}=\mathcal{F}[F_{\mu}]$, where $\mathcal{F}$ denotes the Fourier transform [Eq. (39)].

The AC magnetic field sensitivity $\eta_{\text{ac}}(f)$, characterizing the minimum detectable signal strength for an AC signal at frequency $f$, scales as

[TABLE]

where $|\tilde{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, $|\tilde{F}_{t}(f)|$ can be understood as the effective signal strength experienced by the driven spins at resonance, namely $|\tilde{F}_{t}(f)|=|\vec{B}_{\text{eff}}(f)|/B_{\text{AC}}$, given by

[TABLE]

Here, $\tilde{\phi}_{\mu}(f)$ is the spectral phase of $\tilde{F}_{\mu}(f)$ along the $\mu$-axis, identified from $\tilde{F}_{\mu}(f)=|\tilde{F}_{\mu}(f)|e^{-i\tilde{\phi}_{\mu}(f)}$. We immediately see that $|\tilde{F}_{t}(f)|\leq\sqrt{\sum_{\mu}|\tilde{F}_{\mu}(f)|^{2}}$ from Eq. (41), with the equality saturated when $\phi=\tilde{\phi}_{x}(f)=\tilde{\phi}_{y}(f)=\tilde{\phi}_{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

[TABLE]

To optimally detect this effective sensing field, spins then need to be initialized perpendicular to $\vec{B}_{\text{eff}}$ to form the largest precession trajectory and maximize signal detection contrast. In addition, to optimize contrast for a projective measurement along the $\hat{z}$-axis, for readout the precession plane should be rotated to contain the $\hat{z}$-axis.

VI.2 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. V.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 the standard sensing sequence XY-8 as well as a new sequence Seq. C we designed for interacting spin ensembles (see Fig. 6). 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 $\pi/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 $\pi/2$ pulses to intentionally inject even parities to counteract the odd parities. An example of such a composite pulse is shown in Fig. 6(b).

The new sensing sequence, Seq. C, has identical spectral responses between different axes, $|\tilde{F}_{x}(f)|=|\tilde{F}_{y}(f)|=|\tilde{F}_{z}(f)|$, leading to the transformation of a bare $\hat{z}$-axis resonant sensing field into the [1,1,1]-directional effective field, $\vec{B}_{\text{eff}}$, in the average Hamiltonian picture. While $\vec{B}_{\text{eff}}\parallel[1,1,1]$ 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. (41)]. 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. DD.2 and Seq. I in Fig. 9 for details).

Utilizing these ideas, we demonstrate in Ref. Zhou et al. (2019) a solid-state AC magnetometer operating in a new regime by surpassing the sensitivity limit imposed by spin-spin interactions at high densities. In addition, our average Hamiltonian approach also helps to identify other undesired effects, such as spurious harmonics Loretz et al. (2015); Wang et al. (2019), 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.

VII 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 Hayes et al. (2014); Ajoy and Cappellaro (2013); Bookatz et al. (2014); Choi et al. (2017b); O’Keeffe et al. (2019); Lee (2016); Haas et al. (2019), such that one can explore a range of interesting phenomena in out-of-equilibrium quantum many-body dynamics, including dynamical phase transitions Lindner et al. (2011); Jiang et al. (2011); Heyl (2018), quantum chaos D’Alessio et al. (2016); Garttner et al. (2017) and thermalization dynamics Nandkishore and Huse (2015); Abanin et al. (2018); Álvarez et al. (2015); Wei et al. (2018a, b); Ho et al. (2017); Choi et al. (2019). 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 Khemani et al. (2016); Else et al. (2016); von Keyserlingk et al. (2016); Yao et al. (2017); Choi et al. (2017a); Zhang et al. (2017); Sacha and Zakrzewski (2018); Rovny et al. (2018); Pal et al. (2018).

Indeed, as can be seen in Eqs. (10-13), we can design the toggling-frame spin operators $\tilde{S}^{z}_{k}=\sum_{\mu}F_{\mu,k}S^{\mu}$ to achieve a nonzero target sum in Eqs. (5,6) and engineer the leading-order average Hamiltonian $H_{\text{avg}}$. More specifically, we show that for the common form of two-body interaction Hamiltonians $H_{\text{int}}=\sum_{ij}J_{ij}^{S}(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y})+J_{ij}^{I}S_{i}^{z}S_{j}^{z}$, the relative strength between Ising and spin-exchange interactions can be tuned by the single parameter $c$ as

[TABLE]

where $c$ captures the imbalanced time evolutions in the toggling frames, defined such that the system evolves under the $\hat{z}$- and $\hat{x},\hat{y}$- axes toggling frames for total durations $cT$ and $(1-c)T/2$ in one sequence cycle $T$. Thus, the Floquet-engineered interaction Hamiltonian $H_{\text{avg}}^{\text{int}}$ now exhibits modified Ising and exchange interaction strengths of $\tilde{J}^{I}_{ij}$ and $\tilde{J}^{S}_{ij}$, respectively. Taking the case of interacting NV spin ensembles as an example Kucsko et al. (2018), where $J_{ij}^{S}=-J_{ij}^{I}$, 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 $H_{\text{dis}}=\sum_{i}h_{i}S_{i}^{z}$ to the Floquet-engineered version $H_{\text{avg}}^{\text{dis}}=\sum_{i}(\vec{h}_{\text{eff}})_{i}\cdot\vec{S}_{i}$, where the effective disorder field $\vec{h}_{\text{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. 7(a). Note that the maximum effective disorder strength $W_{\text{eff}}$ is dependent on $c$ (dashed line in Fig. 7(a), see Sec. S1G SM ). Two representative examples of how to engineer such interaction Hamiltonians in a robust fashion are shown in Fig. 7(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 Nandkishore and Huse (2015); Abanin et al. (2018); Álvarez et al. (2015); Wei et al. (2018a, b); Ho et al. (2017); Choi et al. (2019). Thus, our framework will open up a new avenue for the robust Floquet engineering of many-body Hamiltonians.

VIII Experimental Demonstration

Our framework is generally applicable to many different quantum systems, including interacting electronic spin ensembles, such as NV centers in diamond Doherty et al. (2013); Schirhagl et al. (2014); Awschalom et al. (2013); Dobrovitski et al. (2013); Koehl et al. (2011), phosphorus donors in silicon Tyryshkin et al. (2003); Feher and Gere (1959), and rare earth ions Thiel et al. (2011), conventional NMR systems Álvarez et al. (2015); Wei et al. (2018a), trapped ions Blatt and Roos (2012); Jurcevic et al. (2014); Bohnet et al. (2016); Zhang et al. (2017), and even to emerging platforms of cold molecules Carr et al. (2009); Yan et al. (2013); Bohn et al. (2017) and Rydberg atom arrays Labuhn et al. (2016); Bernien et al. (2017). These different systems will have a variety of competing energy scales and distinct interaction types that determine their dynamics, and will thus benefit from the flexibility of our system-targeted design formalism.

Here, we focus on the experimental implementation and demonstration of our results in an interacting ensemble of NV centers in diamond [Fig. 8(a)], tuned to realize the most general form of interactions, see Sec. S2A SM and Refs. Choi et al. (2017c); Kucsko et al. (2018); oth 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\approx(2\pi)4$ MHz) with modest interaction strengths ($J\sim(2\pi)35$ kHz).

To demonstrate the wide applicability of our pulse sequence design formalism, we tune two groups of NV centers with different lattice orientations onto resonance with an external magnetic field. The corresponding Hamiltonian exhibits all different interaction types (Ising, symmetric and anti-symmetric spin exchange) with disordered, position-dependent coefficients, which represents the most general class of one- and two-body interaction Hamiltonians (see Ref. Kucsko et al. (2018) for more details). Despite the complex form of the interaction, system-targeted pulse sequences designed with our formalism enable a sizable extension of coherence times, as shown in Fig. 8(b). More specifically, we find that while a conventional XY-8 pulse sequence is limited by interactions to a coherence time of 0.9 $\mu$s, and Seq. A (Cory-48) performs even worse in this parameter regime, our pulse sequence Seq. B leads to an extension of the coherence time to 3.0 $\mu$s. This observation is consistent with the theoretical prediction [Fig. 5(c)] for a disorder-dominated system, and hence corroborates the importance of considering the energy hierarchy in designing dynamic decoupling pulse sequences.

In Fig. 8(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. 9 for more details of the sequence).

Fig. 8(c) illustrates the coherence decay profile of driven spins of a single, isolated NV group under these two sequences when the rotation angle is chosen to be 92.5$\%$ of the correct rotation angle (gray dots in Fig. 8(d)); Seq. B does not show any oscillations, while Seq. D shows pronounced oscillations over time, resulting from a residual error term $\delta H_{\text{avg}}$ (see Sec. III.2). This behavior is further confirmed in Fig. 8(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.

IX 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. Furthermore, this approach can be readily interfaced with optimal control to substantially speed up the search of pulse sequences and take higher-order effects into account. Using a dense ensemble of interacting electronic spins in diamond, we experimentally confirm the wide applicability of our framework in systems with the most general form of one- and two-body interactions, thus confirming the generality of our approach.

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 dipoles in those high-spin systems for more effective sensing Choi et al. (2017b); O’Keeffe et al. (2019); Fang et al. (2013); Mamin et al. (2014); Bauch et al. (2018). Higher-order contributions beyond the leading-order average Hamiltonian can also be systematically incorporated using the proposed framework. In addition, our formalism may also be extended to the synthesis of dynamically-corrected gates and other nontrivial quantum operations Khodjasteh and Viola (2009b, a); Khodjasteh et al. (2012). While we have focused on the case of $\pi/2$ and $\pi$ pulses around $\hat{x}$, $\hat{y}$ 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 Khaneja et al. (2005); Doria et al. (2011); Iwamiya et al. (1993); Rose et al. (2018), we may be able to robustly engineer many-body Hamiltonians to create macroscopically entangled states, such as spin-squeezed states or Schrödinger-cat-like states, to be used as a resource for interaction-enhanced metrology beyond the standard quantum limit Cappellaro and Lukin (2009); Choi et al. (2018).

Acknowledgements

We thank P. Cappellaro, W. W. Ho, C. Ramanathan, L. Viola, F. Machado for helpful discussions, and J. Isoya, F. Jelezko, S. Onoda, H. Sumiya for sample fabrication. We also thank A. M. Douglas for critical reading of the manuscript and assistance with numerical calculations. This work was supported in part by CUA, NSSEFF, ARO MURI, DARPA DRINQS, Moore Foundation GBMF-4306, Samsung Fellowship, Miller Institute for Basic Research in Science, NSF PHY-1506284.

Appendix A Average Hamiltonian Theory

Here, we introduce the basic principles of AHT and start by considering a generic time-dependent Hamiltonian for a driven quantum system

[TABLE]

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=NT$, the many-body state is given by $\ket{\psi_{t=NT}}=\mathcal{U}(T)^{N}\ket{\psi_{0}}$ with the interaction picture unitary evolution operator Haeberlen and Waugh (1968)

[TABLE]

where $\mathcal{T}$ denotes time-ordering. Here, $\tilde{H}_{s}(t)$ is the rotated system Hamiltonian in the interaction picture with respect to control fields, given by $\tilde{H}_{s}(t)=U_{c}(t)^{\dagger}H_{s}U_{c}(t)$ with the unitary rotation operator $U_{c}(t)=\mathcal{T}\exp[-i\int_{0}^{t}H_{c}(t_{1})\;dt_{1}]$. The control unitary rotation operator over one period is chosen to be identity $U_{c}(T)=\mathbb{I}$.

AHT allows the identification of a time-independent effective Hamiltonian $H_{\text{eff}}$ such that

[TABLE]

The Magnus expansion of $\mathcal{U}(T)$ with expansion parameter $T$ Magnus (1954) can be used to approximate this effective Hamiltonian as

[TABLE]

where $l$ is the truncation order and $\bar{H}^{(k)}$ is the $k$-th order contribution in the Magnus expansion. The first two terms in the series are

[TABLE]

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_{\text{eff}}$ and approximate the dynamics of the many-body state to an accuracy improving exponentially in $l$ Abanin et al. (2017a, b); Mori et al. (2016); Kuwahara et al. (2016). In the following, we focus on the leading order contribution, corresponding to only retaining $\bar{H}^{(0)}$ in the series.

A general control field consists of $n$ pulses $\{P_{k=1,\cdots,n}\}$ with nonuniform pulse spacing $\{\tau_{k=1,\cdots,n}\}$, as shown in Fig. 1(d). Each $P_{k}$ defines a pulsed unitary rotation, generating a discrete set of rotated Hamiltonians $\{\tilde{H}_{k=1,\cdots,n}\}$, where

[TABLE]

with $\tilde{H}_{1}=H_{s}$. As the interaction-picture Hamiltonian is rotated (toggled) at every pulse, $\tilde{H}_{k}$ are also referred to as the “toggling-frame Hamiltonians” and govern the spin dynamics in their respective free evolution intervals $\tau_{k}$. For infinitely short pulses, the zeroth-order average Hamiltonian, $\bar{H}^{(0)}=H_{\text{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 Details of Average Hamiltonian During Free Evolution Time

In this section, we provide a detailed derivation of Eqs. (10-13) characterizing the various average Hamiltonian contributions. 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. Specifically, disorder and Ising interactions during the $k$-th free evolution period transform as

[TABLE]

in the toggling-frame picture. Here we have used the fact that $F_{\mu,k}F_{\nu,k}=\delta_{\mu\nu}F_{\mu,k}^{2}$, since each column of the matrix $\mathbf{F}=[F_{\mu,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_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}=\vec{S}\cdot\vec{S}-S_{i}^{z}S_{j}^{z}=(\sum_{\mu}S_{i}^{\mu}S_{j}^{\mu})-S_{i}^{z}S_{j}^{z}$, which gives:

[TABLE]

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 to $\tilde{S}^{x}(t)=\sum_{\mu}G_{\mu,k}S^{\mu}$, where $G_{\mu,k}$ satisfies the identity $G_{\mu,k}G_{\nu,k}=\delta_{\mu\nu}G_{\nu,k}^{2}$ and takes on values of $\{0,\pm 1\}$ (see Ref. SM for details of how one can explicitly construct $G_{\mu,k}$). Since the commutation relations between spin operators are conserved under frame transformations, the transformed $\tilde{S}^{z}(t)$ and $\tilde{S}^{x}(t)$ operators uniquely specify the $\tilde{S}^{y}(t)$ operator as:

[TABLE]

where $\epsilon_{\mu\nu\lambda}$ is the Levi-Civita symbol. Based on this, we can write the transformation of the anti-symmetric spin-exchange interaction term as:

[TABLE]

where we have used the identity presented above for $G_{\nu,k}$, as well as the fact that $G_{\nu,k}$ has only one nonzero element, squaring to 1. Combining these expressions gives the average Hamiltonian terms in Eqs. (10-13).

Appendix C Analysis of Three-Body Interactions

In this appendix, we analyze the decoupling conditions for spin-1/2 three-body interactions in more detail. 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.

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 Moore and Read (1991); Fradkin et al. (1998), and many other topological phases and spin liquids are ground states of such many-spin Hamiltonians Moessner and Sondhi (2001); Levin and Wen (2005). There have also been various proposals for the direct realization of three-body interactions in experimental platforms ranging from cold molecules Büchler et al. (2007) to superconducting qubits Mezzacapo et al. (2014); Chancellor et al. (2017). 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 $\hat{z}$-direction.

C.1 Ideal Pulse Limit

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 of $\tilde{S}^{z}$ that cover the six axis directions $\pm S^{x,y,z}$.

Proof: Consider a generic $N$-body interaction Hamiltonian $H$ and a set of unitary operators $U_{k}$ ($k=1,2,\cdots,K$). The average Hamiltonian over this set is given by

[TABLE]

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 satisfy $\tilde{S}^{z}=U^{\dagger}S^{z}U=(-1)^{\nu}S^{\mu}$ with $\nu=0,1$, while the rotated $\hat{x}$-axis spin operator, $\tilde{S}^{x}=U^{\dagger}S^{x}U$, can take four distinct values that are orthogonal to the $\tilde{S}^{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 by $\tilde{S}^{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 $\hat{z}$ axis such that $U_{1}=U_{z}U_{2}$ (this rotation leaves the interaction picture $\tilde{S}^{z}$ invariant, but changes $\tilde{S}^{x}$). The Hamiltonian under conjugation by $U_{1}$ is then given by

[TABLE]

where we use $(U_{z}^{\dagger})^{\otimes N}HU_{z}^{\otimes N}=H$. This holds because a rotation around the $\hat{z}$-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. $\blacksquare$

With this lemma in hand, we can utilize mathematical results from unitary $t$-designs DiVincenzo et al. (2002); Dür et al. (2005); Emerson et al. (2003); Collins and Śniady (2006); Gross et al. (2007); Ambainis and Emerson (2007); Dankert et al. (2009); Webb (2015); Zhu (2017) to show that a polarized initial state will be an eigenstate of the three-body interacting Hamiltonian after symmetrization along the six axis directions (i.e. the $\pm\hat{x},\pm\hat{y},\pm\hat{z}$-axes), as described in Sec. IV.3 of the main text. This is a consequence of the fact that the Clifford group is a unitary 3-design. However, 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})^{\otimes 4}+(S^{y})^{\otimes 4}+(S^{z})^{\otimes 4}$, which is found to act nontrivially on a generic polarized initial state.

C.2 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. (10-13), we write the $k$-th toggling-frame Hamiltonian as a polynomial in $F_{\mu,k}$:

[TABLE]

where $\mathcal{O}_{\mu,l}$ describes the generic operator form of interactions that transform as the $l$-th power of $F_{\mu,k}$. More specifically, $\mathcal{O}_{\mu,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 $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_{\mu,k}$. For example, we can rewrite the following three-body interaction as

[TABLE]

This is a three-body interaction with $l=2$, since it is proportional to the square of $F_{\mu,k}$. This suggests that during the finite pulse duration between free evolution blocks $k$ and $k+1$, where the interaction-picture operator $\tilde{S}^{z}(\theta)=\sum_{\mu}(\cos\theta F_{\mu,k}+\sin\theta F_{\mu,k+1})S^{\mu}$ with $\theta$ evolving from 0 to $\pi/2$, the corresponding average Hamiltonian can be written as

[TABLE]

where $\tilde{\mathcal{O}}$ contains operators acting during the rotation pulse and is implicitly dependent on $l$ and the indices in the bracket (for example, when $\theta=0$, $\tilde{\mathcal{O}}=\mathcal{O}_{\mu,l}$).

Expanding the polynomial in the bracket of Eq. (60), we shall find contributions corresponding to terms of degree $u$ in $F_{\mu,k}$ and $v$ in $F_{\mu,k+1}$, with $u+v=l$, since $F_{\mu,k}$ and $F_{\mu,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_{\mu,k}$, the decoupling conditions presented in Tab. 1 of the main text can be directly applied. Similarly, for $l=3$, the terms $F_{\mu,k}^{3}$ and $F_{\mu,k+1}^{3}$ 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_{\mu,k}^{2}F_{\nu,k+1}$ and $F_{\mu,k}F_{\nu,k+1}^{2}$ correspond to a generalization of the interaction cross-term decoupling condition [Condition 3 in Tab. 1] described in the main text:

[TABLE]

for each pair of $(\mu,\nu)$. As an example, for the pair of directions $\hat{x}$, $\hat{z}$, we consider all instances in which the $\hat{x}$ and $\hat{z}$ frames appear in the free evolution frames immediately preceding and following a $\pi/2$ pulse, and the above decoupling condition requires that the signs of all $\hat{x}$ 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 D Efficient Sequence Design Strategies

D.1 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 ${\bf F}$. Consequently, at least 3 free evolution blocks are required. A possible realization is (see also Seq. E in Fig. 9)

[TABLE]

However, according to the discussion in Sec. S1E SM , a single cycle evolution under this sequence does not return the transverse spin operators to their original configuration, i.e. $U_{c}(t)\neq\mathbb{I}$. To achieve $U_{c}(t)=\mathbb{I}$, we can use

[TABLE]

with composite $\pi/2$ pulses inserted at the interfaces of free evolution intervals (see Seq. F in Fig. 9). However, neither of the above sequences are robust against on-site disorder or finite pulse duration effects.

As discussed in Sec. V.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. V.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 in Tab. 1] 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 odd-parity frame changes. The odd number of odd-parity frame changes, however, results in a toggling-frame operator at the beginning of the next cycle $\tilde{S}^{z}(T)$ that has opposite sign from $\tilde{S}^{z}(0)$, violating the periodic condition $\tilde{S}^{z}(0)=\tilde{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)=\mathbb{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. 9)

[TABLE]

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. V. This necessitates a frame matrix structure in which the frames along each axis always consist of pairs with opposite sign. When permitting only a single $\pi/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. 9.

D.2 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 sign-modulation structure along each axis (for instance, the sequences in Fig. 6 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 $\pi/2$-pulse implementations to suppress interaction cross-terms [Condition 3].

If only single $\pi/2$ and $\pi$ 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. DD.1). 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.