Hamiltonian Engineering with Constrained Optimization for Quantum Sensing and Control

While quantum devices rely on interactions between constituent subsystems and with their environment to operate, native interactions alone often fail to deliver targeted performance. Coherent pulsed control provides the ability to tailor effective interactions, known as Hamiltonian engineering. We propose a Hamiltonian engineering method that maximizes desired interactions while mitigating deleterious ones by conducting a pulse sequence search using constrained optimization. The optimization formulation incorporates pulse sequence length and cardinality penalties consistent with linear or integer programming. We apply the general technique to magnetometry with solid state spin ensembles in which inhomogeneous interactions between sensing spins limit coherence. Defining figures of merit for broadband Ramsey magnetometry, we present novel pulse sequences which outperform known techniques for homonuclear spin decoupling in both spin-1/2 and spin-1 systems. When applied to nitrogen vacancy (NV) centers in diamond, this scheme partially preserves the Zeeman interaction while zeroing dipolar coupling between negatively charged NV$^{\text -}$ centers. Such a scheme is of interest for NV$^\text{-}$ magnetometers which have reached the NV$^\text{-}$-NV$^\text{-}$ coupling limit. We discuss experimental implementation in NV ensembles, as well as applicability of the current approach to more general spin bath decoupling and superconducting qubit control.


Introduction
From quantum computing to sensing, both the strength and fragility of the underlying quantum systems stem from interactions. Coupling enables gate operations and metrology, but also introduces crosstalk and decoherence. As quantum technologies mature, the demands of tailoring properties and dynamics increase. Hamiltonian engineering describes a family of classical control techniques on quantum systems to achieve a desired state, process, or observable behavior [1]. This encompasses decoupling to eliminate or reduce unwanted interactions [2,3], optimal control using numerical optimization for maximizing the fidelity of quantum computing operations given experimental limitations [4][5][6][7][8], and Hamiltonian simulation [9][10][11] of models for exploring quantum phenomena not otherwise immediately accessible in experiment [12,13].
Solid-state quantum systems are emerging as viable candidates for high-performance sensing [14], both filling existing technology gaps, such as stable vector magnetometry [15,16], and opening new capabilities in areas like high-spatial-resolution sensing [17][18][19]. Ensemble magnetometry with negatively charged nitrogen vacancy (NV − ) centers in diamond [20] has demonstrated sensitivities of 15 pT Hz 1 2 -/ for DC [21] and 0.9 pT Hz 1 2 -/ for AC [22] sensing. In the latter, the application of dynamical decoupling consistent with narrowband magnetometry [23,24] suppresses low frequency fluctuations and inhomogeneity, yielding improved performance. Still, T 2 * dephasing times and T 2 coherence times currently remain orders of magnitude shorter than T 1 spin relaxation times [25], leaving much room for improvement to reach the spin relaxation limited coherence time T 2 =2T 1 .
The spin-projection-limited sensitivity of an ensemble magnetometer consisting of N spins can be intuitively understood by [26] C N where γ is the gyromagnetic ratio, C the measurement contrast, and τ the interrogation time during which the spins precess, typically on the order of the relevant coherence time, i.e. T 2 * for DC magnetometry and T 2 for AC magnetometry. To date, most spin environment engineering for solid-state ensembles has focused on extending NV − coherence time by synthesizing diamond to eliminate other spin impurities, which commonly include other electronic defects (substitutional nitrogen atoms known as P1 centers, neutrally charged NV 0 , divacancies, NVH − ) [27], as well as nuclear spin species ( 13 C) [28,29]. Although isotope engineering is efficient to remove nuclear spin dephasing [30], the interdependence between NV center formation and the incorporation of a number of spin impurity species makes eliminating all spin impurities via diamond growth intractable. Various broadening mechanisms from paramagnetic defects can be mitigated by driving the environment spins [31,32] while inhomogeneous broadening mechanisms like strain variation and temperature fluctuation can be mitigated by double quantum [33,34]. Used independently, these techniques help to identify sources of dephasing; in combination, they can result in over an order of magnitude improvement in T 2 * [35]. NV − -NV − dipolar interactions will likely ultimately limit NV − coherence time [36]. Inhomogeneous interactions between spins in an ensemble limit coherence [37], with each spin experiencing a different local magnetic field B r loc 3 m due to a neighboring magnetic moment μ a distance r away. A rough estimate for the dephasing time T B 1 2 loc * g  shows that it varies inversely to the density of spins n=N/V. While increasing the spin density for a fixed sensor volume V∼r 3 appears to improve (decrease) the sensitivity in(1), its effect on coherence time, and thus sensing time τ, must also be accounted for.
Ramsey magnetometry entails preparing a superposition state, letting it undergo free evolution for some time τ R , then rotating the state to a convenient basis for readout, with the typical scheme shown in figure 1(b). For an 'ideal' ensemble with no interactions between spins or with the environment, the Ramsey signal is a sinusoid with angular frequency ω 0 = γB, whereas inhomogeneous interactions or coupling to the environment, typical of solid state systems, produces free induction decay, or equivalently a broadened spectral line, shown in figures 1(c) and (d), respectively.
Solid state spins are subject to a large number of anisotropic interactions, including chemical shifts, dipolar, and quadrupolar, which typically leads to broad spectral lines in nuclear magnetic resonance (NMR) and electron spin resonance. A line narrowing workhorse of solid state NMR, the WHH-4 pulse sequence [38] decouples homonuclear dipolar interactions to leading order in spin-1/2 systems, and when repeated in reverse order (MREV-8) also compensates for pulse imperfections [39][40][41]. Choi et al recently discovered a six pulse sequence [3], which we refer to as CYL-6, that achieves dipolar decoupling in spin-1 ensembles, effectively generalizing WHH-4. Drawing motivation from this work, we evaluate both sequences in the context of Ramsey magnetometry, and find that neither sequence produces single frequency observable response. Such an exercise makes clear the need to identify both desirable (the magnetic field to be detected) and undesirable (inhomogeneous line broadening) interactions and the corresponding terms in the Hamiltonian.
Given a system's Hamiltonian with desirable and undesirable terms, the goal of Hamiltonian engineering is to find a set of unitary operators that transforms the original Hamiltonian to a target Hamiltonian containing only desirable terms with maximal strength. In this paper, we explicitly formulate Hamiltonian engineering as a quantum pulse sequence search problem with conditions for achievability and optimality (section 2). We cast the pulse sequence search as constrained optimization in which the search criteria are formulated as constraints and the objective function incorporates sequence length and cardinality regularization suitable for linear or integer programming (LP/IP) frameworks (section 3). As a concrete example, we choose the goal of decoupling dipolar interactions in a spin ensemble used for broadband magnetometry (section 4). The undesirable term is the inhomogeneous dipolar interaction, and the desirable term is a Zeeman term from a global magnetic field of interest. Success criteria include a clean, single frequency Ramsey magnetometry signal, and the strength, or frequency scaling, of the resulting oscillations. We present and analyze pulse sequences for qubit and qutrit systems that effectively average the dipolar interaction term in the Hamiltonian to zero, even with inhomogeneous interaction strengths, while preserving the desired Zeeman term with strength multiplied by 1/3. Because these sequences can be applied during the free evolution interval of a Ramsey experiment and preserve the sinusoidal dependence of the observed signal on the free evolution time, albeit with a frequency scaled by some constant, we name this family of sequences 'Homonuclear Ramsey decoupling' (HoRD), and adopt the naming convention HoRD-[qudit type]-[number of pulses]. Table 1 compares this new family of pulse sequences with existing dipolar decoupling sequences for qubits and qutrits. We discuss considerations for experimental demonstration of the HoRD sequence for spin-1 systems in NV ensembles and, with general techniques germane to a variety of systems, point towards a number of possibilities for future work (section 5) prior to concluding insection 6.

Hamiltonian engineering as a quantum pulse sequence search problem
Hamiltonian engineering entails finding a set of operations, realized with a sequence of pulses, that transform a given Hamiltonian into a target Hamiltonian. In this section, we review average Hamiltonian theory for periodically driven systems. We then make the quantum pulse sequence search problem definition concrete, first abstractly, then explicitly, by using the traditional linear algebraic framework employed in quantum mechanics, and then by defining a convenient 'Pauli projection' representation of the problem.

Average Hamiltonian theory
The state of quantum systems can be manipulated with pulsed control, through which some set of quantum operations can be performed. Average Hamiltonian theory [42] provides a mathematical formalism to evaluate the time evolution of a system subject to a pulse sequence. Originally developed for NMR, these techniques have found renewed interest in the context of quantum control and computation [43].
We capture the effect of intermittent driving with short 'delta function' pulses P k on the system separated by time intervals τ k of free evolution during the kth interval. The time evolution operator propagates the state of the system over a cycle time t k n k c  [42] shows that the time evolution operator can be expressed in terms of a time-independent average Hamiltonian H   ) !, where the operator norm of H gives the largest eigenvalue, so making the pulses and intervals short compared to the fastest timescale of the Hamiltonian generally yields better performance. Zeroing higher order terms, which is generally desirable, can also be achieved by symmetrizing a pulse sequence that gives the desired leading order behavior [39][40][41], or incorporating higher order terms in the Magnus expansion into the optimization goal. For periodically driven systems, the so-called Floquet time evolution operator satisfies ( ) a property we will use for efficient time dependent simulation of the pulse sequences discussed in section 4.

The search criteria and the search problem
This identification of H 0 and H 1 is the specification of  , the goal criterion. The zeroing out is accomplished by applying sequences of unitary rotations U k to the system, such that the effective first order 'average Hamiltonian' (see (6) and surrounding discussion) is: where w k , the weight for unitary U k , is proportional to the amount of time spent during evolution in the reference frame determined by U k , known as the toggling frame. The weights w k are positive and nonzero, and typically are integers, but in principle they could be real numbers. The quantum pulse sequence search problem is thus summarized as follows. Given H H H • Achievability: Does there exist any subset seq   Ì of available unitaries, and weights w k such that where β is bounded away from zero by a constant?
• Optimality: What is the shortest sequence such that the above equation holds, with the largest β (ideally β=1), smallest size seq  | |, and smallest total weight w k k å ?
Note that this formalism is completely general to Hamiltonian engineering, and not limited to interaction decoupling. For example, to obtain and maximize a desired interaction term H 2 that does not occur in the original Hamiltonian H, simply replace H 1 with H 2 on the right-hand side of (12). Reference [3] presents necessary and sufficient conditions for both decoupling and engineered interactions.

The Pauli projection representation
The operators H and U arise from elementary dynamics of constituent subsystems. Thus, it is highly convenient to represent all these operators in terms of the d Let these be indexed from k=1 to k=d, and define σ 0 as the identity. Each σ k is a d×d matrix; for qubits, they are the Pauli matrices, and for qutrits they are commonly taken to be the Gell-Mann matrices. Explicit representations as well as useful properties are given in appendix A. These are all Hermitian, i.e. k k T k . Except for σ 0 , these are traceless matrices, and they satisfy a trace orthogonality condition, We shall refer to the σ k matrices as the generalized Pauli matrices. In addition to the trace orthogonality condition, an exceptionally useful fact is that they form a linear basis for any d×d matrix, i.e. any such matrix M can be expressed as a linear combination of the σ k : and because of the trace orthogonality condition, where the normalization factor α d depends on how the σ k are defined. For Hermitian matrices such as H, Let the full quantum system comprise n constituent subsystems. Then the n-fold tensor products of the generalized Pauli matrices This isomorphic map represents H using a set of d 2 n real numbers, which we shall find convenient to vectorize as c H , and call the Pauli projection coefficients.

Matrix-vector formulation of the search problem
The first-order search problem, defined in (12), is conveniently re-expressed using the Pauli projection, by applying the map to each term in the sum, to define Pauli projection coefficients for each of the U k transformed Hamiltonian terms: Recall that the goal is to average away H 0 , but leave H 1 intact. Thus, in terms of these vectors, the search problem becomes two sets of equations: . The matrix-vector formulation relies on our choice of solving the first order achievability of an average Hamiltonian where linear dependence on the pulse sequence dominates. This is well justified if the duration of each pulse and intervals between pulses are sufficiently short.

Pulse sequence search as constrained optimization
Brute force search checks every possible solution, and is guaranteed to find a solution if it exists in the search space. The computational complexity is linear in the number of possible solutions. Because sequences made of composite pulses contain combinations of available pulses, the number of possible solutions scales combinatorially with the number of available pulses. We prune operators that transform the system Hamiltonian in the same way, and provide two relevant examples to establish a sense of scale. The 24 qubit Clifford operators, defined in section 4.3, only produce 6 unique mappings of the Hamiltonian describing dipolar interacting pins (spin-1/2) subject to a global magnetic field (39), allowing us to prune the search set from 24 to 6 operators. The 13 824 Clifford operators between sublevels in a qutrit system produce 558 unique mappings of the corresponding Hamiltonian for spin-1. While this reduces the number of elements in the search table by an order of magnitude, there are over 10 13 and 10 24 combinations for 6 and 12 pulse sequences, respectively. To overcome the potentially prohibitive cost of a brute-force search, we cast the pulse sequence search problem as constrained optimization. Concretely, we formulate the search problem, (23) and (24), as linear constraints, and impose sequence length and cardinality penalties in the linear objective function.

Constraint formulation
The pulse sequence search problem, as defined by (23) and (24) can immediately be cast in terms of the standard framework for constrained optimization using LP. In that framework, we wish to find a vector x such that: where A eq and A ub are matrices, and together with vectors b eq and b ub , specify the optimization problem. For our search problem, from (24), suppose  projects into the subspace of c H 1 which has nonzero coefficients, such that c 0 .
where N U  = | |. Note the overall minus sign on A ub , and the minus sign on b ub : these are present because we actually have a lower bound, but the standard form of the optimization equations expects an upper bound. Similarly, we may define  as the projector into the subspace of Pauli projection coefficients expected to be zero (e.g. Thus, from (23) we have that b 0 eq = , and A eq is the m eq ×N U matrix formed by taking c k 0 ( )as the kth column, i.e.
The solution x to this LP problem would give the weights w k of the unitaries U k used to achieve the desired goal. However, while this formulation brings our search problem into standard form for LP, this construction does not include additional terms which would minimize the number of unitaries, or minimize the total weight employed. These can usually be included using other constraint mechanisms available within an optimization package. Also, the weights w k are typically desired to be integers. With this constraint, the optimization problem becomes a mixed IP problem, which is harder to solve than the relaxed, LP formulation given above.

Optimization formulation
We now introduce an objective function into the formulation. With the goal to maximize the proportionality constant β of the desired Hamiltonian H 1 , one could choose to minimize A x ub . However, with this requirement encoded in a constraint, we use the objective function for cardinality and weight regularization. Incorporation of a cardinality penalty term offers control over the overall size of the set of pulse primitives (dictionary) seq  , where the parameter α may be tuned to weight cardinality versus total sequence length, and u i is an upper bound on | | is a nonlinear function, so to incorporate such a penalty term in a conventional LP or IP framework, we introduce a set of binaries, z i , satisfying 0x i ub i ·z i , where x i can either be general integer or continuous variable, ub i is the upper bound of x i . Consequently, a value of z i =0 implies x i =0. Thus, to add a cardinality penalty we can simply add a term z i i å to the objective.
x z min 30 If the right-hand side assignment of variables is not supported on the designated optimization solver (e.g. Matlab CPLEX version), it is advisable to define an extension of the inequality matrix A ub , with diagonal of ones (for x i ) and a shifted diagonal (by number of elements in x) with b ub as follows where e denotes the Hadamard product.

Dipolar decoupling for broadband magnetometry
We now consider the problem of decoupling inhomogeneous dipolar interactions in spin ensembles while retaining sensitivity to an external magnetic field. This section begins with a model of a spins in a solid subject to a global magnetic field, with dipole-dipole interactions between spins, followed by the definition of decoupling pulse sequence success criteria for ensemble magnetometry. A spin-1/2 qubit model provides geometrical intuition, and we present a pulse sequence that achieves the success criteria. Applying the previously developed optimization formulation to the spin-1 qutrit model generates a family of pulse sequences, and we examine one in particular with the most desirable characteristics for broadband magnetometry.

Spin ensemble model
The dynamics of a solid-state spin ensemble are governed by an external magnetic field B and interactions between spins J ij , shown schematically in figure 1(a). The spins in the system are generally subject to an external magnetic field where B z is the projection of the magnetic field onto the quantization axis z, with γ the spin gyromagnetic ratio, and a dipolar interaction term depends on the geometrical orientation and physical properties of the two spins, with r ij the distance between spins, θ ij the angle between the separation vector and the quantization axis, and J 0 the interaction strength.
For noninteracting spins, J ij =0, Ramsey fringes of constant amplitude oscillate sinusoidally in time at a Zeeman frequency ω=γB z , with the Fourier transform producing a spectral line at this frequency. Averaging over Gaussians with linewidths J ij drawn from a distribution of coupling strengths representative of an ensemble [45] results in free induction decay of the Ramsey fringes in time and broadening of the spectral line. Figures 1(c) and (d) compare the free induction decay and corresponding frequency spectrum, respectively, of these limiting cases. We take γB z = 2π, Γ = 2π × 10 −2 , and average over 10 000 values of pairwise interactions J ij for each case. Time dependent simulations obtained using QuTiP [46,47] calculate the state at time t=ndt using the time evolution operator U(dt), given by (3), along with the Floquet property, (8).

Success criteria
We define the success criteria for dipolar decoupling in spin ensemble magnetometry as follows.
• Clean Zeeman-the Ramsey signal of an experimental observable is a pure, single frequency sinusoidal oscillation.
• Zeeman Strength-the scaling factor β of the Ramsey signal frequency of an experimental observable with respect to the same observable under ideal Zeeman evolution.
These criteria express achievability and optimality conditions defined in section 2.2, including (12). The search succeeds when it finds a set of unitaries that transforms the system Hamiltonian, (39), to one proportional to the desired Zeeman terms H Z only, and averages the undesired dipolar coupling terms H dd to zero, which indeed results in an experimentally observable clean Zeeman signal. Optimality includes maximizing the scaling factor β of the resulting desired terms This is precisely the scaling of the Ramsey fringe frequency with respect to a given magnetic field.
We will find it useful to compare the Zeeman strength of decoupling sequences that achieve the clean Zeeman criterion with existing dipolar decoupling sequences that do not meet the clean Zeeman requirement. The Zeeman strength for a Hamiltonian that is not a clean Zeeman is calculated by normalizing the projection onto the generalized Pauli basis, and taking the inner product with S z .

Qubit ensembles
For spin-1/2 qubits, the spin matrices S k =σ k /2 are proportional to the Pauli matrices σ k for k=x, y, z. Rotation operators about each axis by an angle θ are and for π/2 rotations we adopt the following notation analogously defined for y and z. The set of rotations that maps Pauli operators to Pauli operators (up to a multiplicative factor ±1) forms a group, known as the Clifford group. Let us formalize this set by defining the operators, from which we may construct any of the 24 Clifford rotations from the product V i W j . Thinking of the Clifford group as the set of operations that maps the x, y, z axes, or corresponding Pauli operators, to distinct orientations along the cardinal axes ±x,±y,±z axes provides geometrical intuition. Consider, for example, that the six V i map σ x to unique cardinal directions, and the four W j orient σ y and σ z with respect to σ x . The WHH-4 sequence [38] X Y Y X 2 4 7 t t t t t¯( ) known from NMR decouples homonuclear dipolar interactions to leading order. The instantaneous pulses punctuate free evolution time intervals τ i =τ for i=0, 1, 3, 4 and τ 2 =2τ. The time ordering of the free evolution times and pulses goes from right to left. The corresponding unitary basis transformations Note the zero frequency component in addition to Ramsey oscillations scaled in amplitude by a factor of 2/3, and in frequency by 1 3 . The frequency scaling is consistent with calculation of the Zeeman strength 1 3 b = obtained by normalizing the projection of the Hamiltonian onto the Pauli matrices and taking the inner product with S z =σ z /2.
While WHH-4 works in practice for homonuclear decoupling in NMR, given that low frequencies are effectively filtered out, for magnetometry, it may prove suboptimal. Applying a decoupling pulse sequence that effectively decouples dephasing interactions increases the available measurement time τ, but also alters the effective gyromagnetic ratio γ and measurement contrast C must also be accounted for. The magnetic field sensitivity, (1) is inversely proportional to the maximum slope γC of the Ramsey signal. For WHH-4, C 0.385 where γ and C capture the scaling of the gyromagnetic ratio and contrast, respectively, due to the pulse sequence. In cases where the contrast noise is independent of the signal amplitude, e.g. photon shot noise limited optically detected magnetic resonance, reducing the contrast is generally undesirable. The zero frequency component may limit the use of WHH-4 in low or zero bias field operation, of current interest for certain applications [48,49].
The following dipolar decoupling for broadband magnetometry sequence, which we call HoRD-qubit-5, t t t t t t¯¯( ) found analytically using average Hamiltonian theory, yields a single frequency Ramsey oscillation ( figure 2). The corresponding basis transformations  [38] and no decoupling shown for comparison. The HoRD-qubit-5 pulse sequence decouples dipolar interactions in spin-1/2 ensembles, resulting in a constant amplitude single frequency sinusoid (clean Zeeman) with frequency scaled by a factor of β=1/3 (Zeeman strength) relative to ω/ω 0 =1, where ω 0 =γB z . By comparison, the WHH-4 sequence does not produce clean Zeeman behavior, as it reduces the amplitude of the sinusoid and produces a DC offset, resulting in a spectral component at ω=0 in addition a component scaled by a factor of with free evolution for intervals τ i =τ for all i averages H dd to zero while scaling the strength of H Z by a factor of 1/3. With no change in contrast, this gives a maximal slope C 0.333 1 3 g = » slightly lower yet comparable to WHH-4. The HoRD-qubit-5 pulse sequence also zeros the second order average Hamiltonian, (7). These unitary operators are in fact Clifford operators, V W i j , defined in (46), where relations such as R XR X z y  Figure 2(a) depicts the orientation of the Pauli operator axes during each interval, taking advantage of the mapping of SU(2) rotations in spin space to SO(3) rotations in real space. Note that any sequence that orients the +z axis along ±x, ±y, ±z, for equal intervals achieves the desired leading order decoupling for a static relative arrangement of dipoles. This freedom to tailor the pulse sequence may prove useful in canceling higher order terms or simplifying the experimental control pulses. Time-dependent simulations show that HoRD-qubit-5 results in a single frequency Ramsey signal S S cos 3 z z tot tot 0 R w t á ñ = ( )and corresponding spectrum with a single resonance at ω/ω 0 =1/3 (figures 2(b) and (c), respectively).
The HoRD-qubit-5 sequence would find application in magnetometry with dense spin-1/2 ensembles for which dipolar coupling between spins is a limiting factor for phase coherence. The negatively charged silicon vacancy center (SiV − ) in diamond is a spin-1/2 color center [50], for which coherence of a single SiV − is limited by coupling to a spin bath of substitutional nitrogen atoms at sufficiently low (millikelvin) temperatures [51], and by phonon-mediated interactions at higher temperatures [52]. Note that while the qubit approximation to higher spins, in which only a pair of levels is considered, can describe the Zeeman splitting (with the appropriate gyromagnetic ratio), it fails to capture physics of dipole-dipole coupling.

NV ensemble magnetometry
Demonstrations of sensitive broadband magnetometry [21,53], and recent progress elucidating and eliminating dominant dephasing mechanisms [35] in these systems, make dense NV spin ensembles particularly interesting for applying the dipolar decoupling techniques investigated in this work. Here we describe some essential physics of NV ensembles relevant to magnetometry, discuss how these systems relate to the spin ensemble model in . B z is the component of the magnetic field projected onto the NV − axis, D≈2.87 GHz is the zero-field splitting, and S z is a spin-1 operator. Figure 3(a) shows the NV energy levels. In an ensemble containing all four NV orientations, the magnetic field projects differently in principle onto each axis. The model considered here corresponds to the magnetic field oriented along a single NV − class. This can be realized by optically initializing all NV − orientations to the 0ñ | state, then preparing a superposition of only a single chosen orientation.
Tailoring the growth process to achieve preferential NV orientation along two axes [54] or a single axis [55][56][57][58] has resulted in increased measurement contrast in dilute [59] and dense [60] ensembles, yet the irradiation and annealing used to optimize NV − concentration destroys preferential orientation. Note that hyperfine interactions introduce beats in a typical NV − Ramsey signal; the current analysis considers only the electronic spin degrees of freedom. Other color centers to which the current result applies include neutrally charged silicon vacancies [61,62] in diamond and divacancies in silicon carbide [63].
To connect the NV model to the spin ensemble model, (39), we analyze the effect of the zero field splitting. Moving to the rotating frame of H D leaves H Z unchanged, and shows that the large zero-field splitting of the NV suppresses certain dipole-dipole matrix elements. This makes any sequence that decouples dipolar interactions in a qutrit ensemble model applicable to NV − ensembles, whereas the converse is not necessarily true.
The spin-1 property of the NV − provides three choices of superposition state for the Ramsey measurement, enabling both single quantum (SQ) and double quantum (DQ) magnetometry. SQ magnetometry uses a superposition of 0ñ | and 1 + ñ | , or 0ñ | and 1 -ñ | states. Because the zero-field splitting D has temperature dependence and is broadened by inhomogeneous strain, this basis is suboptimal for broadband magnetometry. DQ magnetometry employs a superposition of 1 + ñ | and 1 -ñ | states in order to cancel out common-mode shifts due to temperature fluctuations or strain inhomogeneity. The doubled effective gyromagnetic ratio results in a Ramsey frequency that exhibits a two-fold increase, and dephasing due to dipolar broadening increases by a factor of two for spin bath coupling and four for NV − -NV − coupling, which can be seen in figures 3(b)-(e). This amplifies the requirement for a pulse sequence to produce an average Hamiltonian containing the full S S z i z i tot = å Zeeman term, in order to produce a clean Zeeman Ramsey signal for any of the three superposition state bases.
The quantum states of a given class of NV − in an ensemble can be controlled by delivering resonant microwave pulses to the entire sample. Transitions between 0 and±1 are driven directly at frequencies D±γB z , whereas transitions between±1 states require composite microwave pulses. Current experimental setups include multiple microwave drive channels, each with independent phase and amplitude control, and well-tuned π/2 and π pulses suitable for SQ and DQ magnetometry [35,64].

Qutrit ensembles
With the physics of NV systems, relationship to the qutrit ensemble model, and available controls established, we now develop an explicit model for the spin-1 quantum pulse sequence search problem. The increased dimensionality of spin-1 qutrit ensembles both affords greater flexibility for Hamiltonian engineering and requires a more systematic approach to push beyond the geometrical intuition apparent for spin-1/2 systems.
Consider Clifford rotations between pairs of levels {+1, 0}, {0, −1}, {+1, −1} with k=1, 2, 3 corresponding to each qubit subsystem, respectively, The Gell-Mann matrices λ k , k=1, K, 8, given explicitly in appendix A, provide an orthonormal basis for traceless Hermitian 3×3 matrices, and are the generators of SU(3). The rotation operators  |  ) ) bases, respectively, of a Ramsey experiment with the HoRD-qutrit-8 pulse sequence applied during the free evolution time, with CYL-6 [3] and no decoupling shown for comparison. The HoRD-qutrit-8 pulse sequence decouples dipolar interaction in a spin-1 ensemble and produces a clean Zeeman signal with frequency scaled by a factor of the Zeeman strength β=1/3 relative to ω/ω 0 =1 (ω/ω 0 =2) for SQ (DQ) magnetometry, where ω 0 =γB z . The CYL-6 sequence also performs dipolar decoupling, but produces a Ramsey signal containing multiple frequency components which also depend on the choice of basis. Note that without decoupling, inhomogeneous interactions produce a faster decay and broader linewidth for DQ compared to SQ. k exp i 2 , 1, 2, 3, 51 for each subsystem are built from the corresponding SU(2) subalgebra, (A.6), where ζ and η are linear combinations of λ 7 and λ 8 . The Clifford group for subsystem k is defined by using the rotation operators analogously defined for y and z. In the definitions of {V i , W j } in (46), with the substitution V W VW i j i j k  ( ) . While we consider only Clifford rotations here given the prevalence of pulsed quantum control of NV ensembles [35], this restriction is not fundamental. Quantum control techniques that are robust against noise, such as fast holonomic gates [65], which have recently been demonstrated in NV systems [66], present an interesting avenue for future work.
With a clean Zeeman target Hamiltonian, we search for and find a set of unitary operators. We harness the IBM ILOG CPLEX [67] library that handles LP, IP, and constraint programming (CP). shows 558 unique mappings, allowing us to significantly prune the search set. We anticipate sequences containing at least six unitaries, such as CYL-6, needed to decouple dipolar interactions and possibly longer sequences to obtain a clean Zeeman term.
We discover a family of sequences for qutrit ensembles that decouples dipolar interactions with various resulting Zeeman terms, which we describe here. LP finds the set of twelve unitary operators, given explicitly in appendix B. While this set does zero dipolar interactions, it also undesirably zeros the Zeeman term exactly . This is a clean Zeeman with Zeeman strength 1 6 . To maximize the Zeeman strength, we systematically turn on terms to which the Zeeman Hamiltonian maps, then look for a unitary operator from the search space that maps the average Hamiltonian back to S z . The Zeeman strength resulting from a set of unitary operators should increase with fewer intervals canceling each other out, so we want to maximize the number of Zeeman interval pairs that are on.
We find that the sequence {U 0 , U 2 , U 4 , U 5 , U 6 , U 8 , U 10 , U 11 } with weights {2, 2, 1, 1, 2, 2, 1, 1} produces an average Hamiltonian with terms proportional to λ 1 +λ 2 +λ 4 +λ 5 . A brute force search of products of twolevel Cliffords that maps this to S z finds U V W V W V W   , with Zeeman strength β=1/3. We call this sequence HoRD-qutrit-8 and it is one of the main results of this work. The corresponding pulses may be found using P U U The Zeeman strength 1 6 0.408 b = » of CYL-6 is slightly larger than 1 3 0.333 b = » of HoRD-qutrit-8. However, resonances of the observable Ramsey signal of CYL-6 depend on the choice of initial state (see figure 3 for a comparison of SQ and DQ initial state), unlike those of HoRD-qutrit-8, hindering direct comparison of the Zeeman strength.
An optimal pulse sequence will achieve the largest Zeeman strength β1 with the fewest number of pulses seq  | | and smallest total weight w k k å . Proving optimality is difficult in practice. Brute force search can place a lower bound on sequence length, but becomes computationally expensive after a few pulses due to combinatorial explosion, as discussed in section 3. Analyzing known sequences that decouple dipolar interactions but do not produce a clean Zeeman term provides some insight. The CYL-6 sequence, and a second six pulse sequence which can be constructed from HoZD-qutrit-12, (B.1), by taking only even (or odd) terms, {U 0 , U 2 , U 4 , U 6 , U 8 , U 10 } all with equal weight, are the shortest currently known. As neither of these sequences produces a clean Zeeman term, sequence length optimality of the eight pulse HoRD-qutrit-8 seems plausible, if not provable.

Discussion
We now discuss experimental realization of the HoRD-qutrit-8 pulse sequence in NV ensembles, opportunities for future work on spin ensemble sensing, and potential application of the current Hamiltonian engineering approach to address crosstalk in superconducting qubit devices for quantum computation.

Experimental realization in NV ensembles
A successful demonstration of NV − -NV − decoupling using HoRD-qutrit-8 in NV ensembles with sufficiently high nitrogen density must address the following challenges; similar considerations apply to other solid state color centers such as the silicon vacancy center in diamond. Prolonging the inhomogeneous dephasing time T 2 * in experiments requires a 'shoot the alligator closest to the boat' approach, which entails mitigating the most dominant mechanism first, then systematically addressing each next-dominant factor. Consider the contributions of independent dephasing mechanisms to T 2 *, where T 2 *{·} denotes the dephasing time limit due to a specific mechanism [35]. The P1 spin bath, strain inhomogeneity, and temperature fluctuations likely dominate NV − dephasing in high nitrogen density ensembles. Successively weaker interactions likely include NV − -NV − , NV − -NV 0 , and others which are less well understood. P1 driving [32] to eliminate dephasing due to the substitutional nitrogen spin bath, and double quantum magnetometry [33,34] to eliminate broadening due to strain inhomogeneity and temperature fluctuations, should be employed together [35] with the HoRD-qutrit-8 dipolar decoupling protocol for magnetometry. State preparation and manipulation using all three spin-1 basis states has been achieved for Ramsey magnetometry using multi-frequency pulses [35], and for nanoscale temperature sensing using composite pulses between the 0 and±1 states to effect a pulse between the −1 and +1 states [64]. While the HoRD-qutrit-8 sequence requires a broader library of pulses compared to these state of the art results, leveraging the same fundamental methods that employ amplitude, phase, and frequency control provides a feasible pathway to implementation. Control of the NV − spin-1 ground state manifold can be accomplished with microwave driving. The optimal time for a Ramsey magnetometry measurement is on the order of the inhomogeneous spin dephasing time T 2 *. P1 spin bath driving combined with DQ magnetometry has resulted in improving native T 2 2 * » μs to T 30 2,DQ Drive » + μs [35], which we use to estimate the power required for fast pulses. A Rabi frequency of 2 7.7 MHz R p W =´has been obtained using a loop gap resonator with incident microwave power P 16 W » [68], with Ω R 10 MHz commonly achieved using a wire loop antenna applying a homogeneous field over the sample detection volume [59]. Consider the HoRD-qutrit-8 sequence applied once during a total measurement time τ R =32 μs, with composite pulses consisting of six 50 ns long pulses punctuating free evolution intervals τ≈4 μs. This makes the pulses short compared to the free evolution intervals, and the timescale for one cycle t c R t  short compared to the characteristic dipolar interaction strength J 0 =2π×16 kHz/ppm, even at nitrogen concentrations exceeding 10 ppm. Using multi-tone driving [35] and compressing the pulses, a subject of ongoing investigation, could shorten each composite pulse to 50 ns, making it feasible to repeat the pulse sequence multiple times during the same total measurement time while limiting the accumulated pulse duration to 10% of the free evolution time. For comparison, dynamical decoupling sequences with eight pulses have been performed in less than 2 μs, and sequences of up to 256 pulses have been demonstrated to increase the coherence time T 2 to millisecond time scales, approaching the spin-lattice relaxation T 1 limit [23].

Future work
Avenues of further inquiry for spin ensemble sensing include spin bath decoupling, the effect of multiple sensing spin orientations, and interaction-enhanced metrology. Simultaneously driving multiple bath spin and sensing spin resonances does begin to crowd the frequency space and require additional control electronics. An additional direction of research will focus on spin bath decoupling for broadband magnetometry using only NV and bias field control, without directly driving the spin bath. While the model considered here describes a single NV orientation or a preferentially oriented NV ensemble with a co-oriented bias field, many ensembles contain NVs equally distributed among the four crystal axes, a critical characteristic for vector magnetometry. The current constrained optimization approach will readily incorporate a model with all four NV orientations, yet the existence and, if found, practical utility of such a sequence both remain open questions. Interactions play a key role in quantum-enhanced measurements beyond the standard quantum limit, which generally require highly non-classical states with entanglement or squeezing [69,70]. Engineering rather than eliminating interactions represents a straightforward extension of the current formalism, although high fidelity entangled state preparation in large systems remains an experimental challenge. Ancilla-assisted frequency upconversion presents an alternative scheme for DC magnetometry [71].
To illustrate the generality of this approach, we briefly discuss a second example that arises in superconducting qubit circuits with multiple qubits coupled to the same resonator [72]. Optimizing the fidelity of resonator-mediated two qubit gates between a given pair of qubits requires maximal effective coupling between this pair while eliminating, or reducing to a negligible level, crosstalk with other qubits [73]. Definitions of success for the qubit case such as 'exclusive coupling' and 'coupling strength' yield appropriate optimization criteria, with fidelity under realistic conditions providing an ultimate performance metric.

Conclusion
We cast Hamiltonian engineering as a quantum pulse sequence search problem, identifying desired and undesired terms within the system Hamiltonian, and using average Hamiltonian theory to define achievability and optimality conditions. Defining orthonormal projection coefficients of transformed Hamiltonian terms, we provide an explicit matrix-vector formulation of the search problem, which translates directly to linear equality and inequality constraint formulation. Because minimizing both the sequence length and cardinality is desirable in practice, the optimization objective function includes both of these quantities. We explicitly show how to formulate the optimization for LP or IP. We apply this formalism to the problem of dipolar decoupling for broadband magnetometry in spin ensembles. The goal is to decouple inhomogeneous interactions while retaining susceptibility to an external physical quantity of interest in order to prolong coherence in a way that directly translates to improved sensitivity. This motivates the definition of success criteria for achievabilityclean Zeeman, and optimality-Zeeman strength. The HoRD-qubit-5 pulse sequence for spin-1/2 qubits provides an intuitive example and comparison with the well-known WHH-4 sequence shows its utility. We apply the constrained optimization formalism using LP in IBM CPLEX to discover the HoRD-qutrit-8 pulse sequence for qutrit ensembles, which achieves a clean Zeeman average Hamiltonian with Zeeman strength β=1/3, the same Zeeman strength as HoRD-qubit-5. This sequence can be used in both the single quantum and double quantum bases, making it compatible with existing common-mode noise rejection techniques and spin-bath decoupling protocols. This work represents the first time, to our knowledge, decoupling techniques compatible with Ramsey magnetometry have been presented and analyzed, extending multiple pulse sequences from AC to DC sensing. Additionally, the NV − -NV − decoupling approach proposed here overcomes a previously-viewed fundamental limit posed by sensing spin interactions.

Appendix. B
The unitary operators comprising the HoZD-qutrit sequence, discussed insection 4.5, are with equal weights w k =1 for all k.