Programmable quantum simulation by dynamic Hamiltonian engineering

Quantum simulation is a promising near term application for mesoscale quantum information processors, with the potential to solve computationally intractable problems at the scale of just a few dozen interacting quantum systems. Recent experiments in a range of technical platforms have demonstrated the basic functionality of quantum simulation applied to quantum magnetism, quantum phase transitions, and relativistic quantum mechanics. In all cases, the underlying hardware platforms restrict the achievable inter-particle interaction, forming a serious constraint on the ability to realize a versatile, programmable quantum simulator. In this work, we address this problem by developing novel sequences of unitary operations that engineer desired effective Hamiltonians in the time-domain. The result is a hybrid programmable analog simulator permitting a broad class of interacting spin-lattice models to be generated starting only with an arbitrary long-range native inter-particle interaction and single-qubit addressing. Specifically, our approach permits the generation of all symmetrically coupled translation-invariant two-body Hamiltonians with homogeneous on-site terms, a class which includes all spin-1/2 XYZ chains, but generalized to include long-range couplings. Building on previous work proving that universal simulation is possible using both entangling gates and single-qubit unitaries, we show that determining the"program"of unitary pulses to implement an arbitrary spin Hamiltonian can be formulated as a linear program that runs in polynomial time and scales efficiently in hardware resources. Our analysis extends from circuit model quantum information to adiabatic quantum evolutions, where our approach allows for the creation of non-native ground state solutions to a computation.


Introduction
Universal computation is derived from the ability to reduce complex and arbitrary calculations to a finite set of physical operations. While generic digital logic is well suited for many applications, analog simulation can provide much more efficient methods for performing certain types of computations. In this approach, computational operations of interest can be efficiently built into the underlying hardware, for instance, as analog circuits. Such approaches can greatly enhance performance at the cost of guaranteed error tolerance [17].
The same capabilities exist in the quantum regime, with the promise to realize extraordinary computational capabilities for a variety of challenging problems of interest. Analog quantum simulation may provide a near-term path to providing solutions to currently intractable computations at the scale of just a few dozen interacting quantum systems. Quantum simulation leverages well-controlled quantum coherent devices in order to study the properties of poorly understood many-body systems: e.g. interacting spins and quantum magnetism. In this vein, a wide variety of experimental demonstrations have been realized to date [5,6,3,10,9] using precisely controlled atomic systems.
The efficiency benefits associated with analog simulators are partially offset by the challenging open question as to how one may efficiently implement a "program" therein using a constrained set of basic resources. This is apparent in the case of experimental simulation of quantum magnetism; all experiments to date have been severely constrained by the native interactions achievable between the underlying quantum bits. Attempts to expand functionality have typically focused on exploiting the specific properties of a given hardware platform [18,19] (e.g. varying the coupling of spins to different normal modes of motion in a trapped-ion simulator), rather than providing a generic, technology-independent approach. Producing a general, efficient, and extensible framework for programming complex interactions and broadening the range of accessible simulations with limited hardware capabilities is thus a key requirement to expand the utility of mesoscale quantum simulation.
Here we describe a technology-independent technique to develop a hybrid quantum spin simulator with constrained resources, capable of efficiently realizing a broad class of interacting spin Hamiltonians. Beginning with a native long-range pairwise interaction between spins on a one-dimensional chain, the application of short sequences of singlespin Pauli operators can effectively generate any other Hamiltonian in the family by filtering the time-averaged relative weight of different pairwise interactions. The approach we describe frames the problem as one of functional synthesis where basis functions are generated by the pulse sequences that we define. This general methodology has met success in a variety of contexts for quantum information, including the design of algorithms for quantum memory [20,21], robust-control techniques [22,23], circuitmodel algorithm synthesis [24], and even applications in quantum-enabled sensing [25].
Through sequential and concatenated application of the pulse sequences we define, the basis functions can be added and multiplied together and are shown to generate a complete, orthogonal set of functions that span the full space of Hamiltonians in the family. The challenge of determining the appropriate operations required to map the native interaction to any other is reduced to an efficiently soluble linear program, and moreover this solution is (in a certain precise sense) optimal. In addition to providing specific examples of useful filtering tasks, we show how this technique is extensible to adiabatic quantum simulation and permits the preparation of states that are otherwise inaccessible. Specifically, we consider the family of one-dimensional Hamiltonians with translation-invariant symmetrically coupled two-body terms (i.e. XX-, Y Y -, and ZZtype couplings) extended to include long-range couplings, and homogeneous transverse on-site terms. More precisely, we model a finite N -qubit chain on a one-dimensional lattice ( Fig. 1a) with Hamiltonian H = H I +H T . Here H I is a fully connected interaction Hamiltonian where represent the interaction strength and coupling in each Cartesian direction of two qubits separated by d lattice spacings. The dot product between these implements a generally anisotropic coupling. We also incorporate a homogeneous transverse field H T = B j Z j , where the magnitude of B is controllable.
Our first main result is that, starting with any such family with polynomially decaying long-range couplings Ω d , there exists a short sequence of single-qubit Z pulses such that the time-averaged Hamiltonian is equal to any other Hamiltonian in the class, up to rescaling factor. We can quantify this statement as follows. The sequence of Z pulses contains at most O(log N ) concatenation steps, requiring a pulse number of at most O(N 2 ) with at most O(N 3 ) individual spin flips, and the rescaling factor for the filtered interaction strength is at most N −2.585 , the exponent given by − log 2 6. Thus the scheme can be fully implemented in polynomial time. Moreover, the optimal rescaling factor for our scheme can be obtained by efficiently solving a linear program. We also provide explicit demonstrations of our method for implementing power-law decay modulation in Ising spin chains and adiabatic state preparation, as well as explicit bounds on simulation error.
The remainder of this paper is organized as follows. We describe our general approach to Hamiltonian engineering in Sec. 2, introducing an extensible family of pulse sequences that allow us to generate any Ising-type coupling (in the above sense). Through combination of filters we demonstrate how these pulse sequences are in fact universal, allowing arbitrary mappings within the Hamiltonian class with efficient resource scaling in Sec. 3, with additional proof details in Appendix A. We then show how filters may be efficiently compiled to achieve a desired "program" for a quantum simulator in Sec. 4, followed by examples of what can be achieved with the method in Sec. 5. We then show how this can all be generalized to adiabatic evolutions (Sec. 6) and consider the fully general Hamiltonians of Eq. (1) in Sec. 7, including rigorous bounds on the simulation error. Our work concludes with a discussion of open questions in the literature and a summary of our results.

General Concepts
Owing its foundations to research in NMR [26,4], the idea of altering the dynamics of interacting quantum systems has proven invaluable in the field of quantum control. For instance, in the context of dynamic error suppression [27,28,29,30] a quantum system may be effectively isolated from its environment, not by eliminating the physical interaction, but by inducing a dynamical response which acts to average out the coupling to the environment. This picture extends to interactions between quantum coherent systems as in spin decoupling and recoupling in NMR [4,12,13,31] and have been used to provide universality proofs of quantum simulators using a set of operations similar to universal quantum computation [14,15,16].
The simplest example of this type of Hamiltonian engineering is the dynamic decoupling of two spins interacting via the Hamiltonian H = ΩX 1 X 2 , where X j is the Pauli x operator for the jth particle. The application of phase-flip pulses (Z j ) at time intervals ∆t produces time-dependent operators X j (t) = w j (t)X j in the Heisenberg picture. The control propagator w j (t) = ±1 captures the discrete-time modulation with the sign changes occurring at the times of the applied pulses and time evolution operator is, with Ω eff = Ωw 1 · w 2 and [w j ] l = w j (l∆t). We thus see that the inner product of the discretized control propagators w j determines the effective coupling. It has been demonstrated that the physics underlying this simple example can be generalized for application to a collection of qubits in order to realize arbitrary effective two-qubit interactions [12]. Unfortunately existing approaches entail tremendous complexity when employed for Hamiltonian engineering in quantum simulations. Our goal is to produce a simple extensible set of pulse-sequence constructions suitable for the task of quantum simulation. Our method takes advantage of the symmetry of a ubiquitous architecture in physical quantum simulators, the translational invariant onedimensional spin chain with long-range coupling.

Filter construction
We now show how to build on this very simple example and define a suite of pulse sequences that may be used to efficiently program an effective long-range pairwise spincoupling Hamiltonian in a multi-qubit system. Consider a Hamiltonian of the form in  Eq. (1); this form of long-range spin interaction is extremely important for quantummany-body physics studies of exotic forms of quantum magnetism and is commonly realized e.g. in trapped-ion spin simulators and optical lattice simulators [32,33,34,5,35].
We begin by treating the special case of Ising-type interactions, and determine the relevant pulse sequences that permit our desired transformations on the interqubit couplings to be realised exactly. Later we will revisit the generalised anisotropic Heisenberg-type coupling with homogeneous on-site terms and discuss what modifications must be made in order to effectively treat the more general model. We look to craft a set of pulse sequences that permit the realization of an arbitrary effective Hamiltonian within this class by modifying the native interaction so that In this approach both the original and effective Hamiltonians are translationally invariant.
To simulate the time-evolution under a Hamiltonian for a time T , we break T into k segments of equal duration ∆t k = T /k, with pulses applied as required at the transitions between time segments [15]. This results in a family of control propagators {w (k) j } labelled by k with j indexing the qubits. We define operatorsP k,l that are applied both at the beginning and end of the lth time segment ‡. Each is a product of Z j operators on some subset of the N different qubits in the chain. The result of applying these pulse sequences is a modification of the time evolution during the lth time segment to U k,l (∆t k ) =P k,l exp −iH∆t k P k,l . TheP k,l operators generate sign changes in the control propagators of individual qubits and the sign ofP k,l X jPk,l is encoded by the lth element of w The total evolution operator is given by the product over all time periods as The dynamically modified coupling is represented through the function f k (d), which we define to be the vector dot product of the control propagators w j+d over the discrete time segments.
With this framework we introduce the primary class of pulse sequence we construct, Λ k , defined by the pulse operators Here, square brackets around the operators P k,j indicate that we are defining a set of operators acting simultaneously on all qubits in the chain. We also denote rounding x down (respectively, up) to the nearest integer by x (respectively, x ). In this construction the exponent on the ith Pauli Z operator is always an integer so that, depending on the parity, either Z i or the identity is applied to the ith qubit and j denotes the relevant time bin. These sequences modify Ω d by the factor where {x} denotes the fractional part of x. This filter varies as a triangle wave as a function of qubit distance (as shown in Fig. 2a), with period 2k and range f Λ k (d) ∈ [−1, 1] which allows for considerable flexibility in the effective coupling. For example, it is possible to switch the sign of the interaction from ferromagnetic to antiferromagnetic, or visa versa, as a function of distance, d.
The second class of filters we introduce is defined as Γ k : The filter Γ k is capable of providing a relative enhancement to the interaction strength of qubits separated by integer multiples of k lattice spacings. The Γ k map modifies the system's evolution according to Because of the selectivity of the filter's action in providing a periodic (in d) relative enhancement in coupling, we refer to this as a boost filter.

A complete set of interactions
The simulation of arbitrary Hamiltonians mandates much more complex functional dependences than simple triangle waves or periodic boosts. Previous work [36] has shown that triangle waves coupled with phase shifts form a complete basis for functional synthesis similar to sines and cosines. Despite being unable to generate phase shifted triangle waves, we show that it is possible to achieve arbitrary multiqubit interactions when we augment positive linear combination of filters with products of filters. Utilizing both operations overcomes the limitations arising from accessing only positive linear combinations and "cosine"-like triangle waves.

Universal filter space
A positive linear combination such as, , is generated through sequential application of the pulse sequences with the relative weights being determined by the time taken to implement the pulse sequence; cf. Eq. 4. Similarly, products of filters such as are achieved via sequence concatenation, in which sequences are "nested" within the free evolution periods of one another. These two operations are simple but extremely powerful; with these it is possible to engineer arbitrary translationally invariant qubit couplings within the family we consider. Achieving this among N qubits implies the ability to perform functional synthesis in an (N − 1)-dimensional space, where each dimension represents the coupling between qubits separated by a particular distance, d. A point in this space therefore represents an arbitrary combination of interaction strengths over the relevant values of d. In this space, we define C N to be the convex set of filters f k (d) generated by both positive linear combination and concatenation of our basic filters. While this set takes a complex geometric form, we will show that it contains an (N − 1)-dimensional hypercube around the origin, implying that these filters achieve an arbitrary f (d) and hence an arbitrary Ω eff (d). We call this cube "universal", in analogy with universal computation, but here Traces plotted here have k = 1 (black, highest frequency) to k = 8 (violet, lowest frequency). Lower Panel: Basis of boost filters arising from f Γ k (d), for k = 3 to k = 6 using the same color encoding as in the upper panel. b) A representation of the set C 3 and the universal filter space contained within it, represented by the square confined in size to fit within the convex set defined along axes Ω d=1 and Ω d=2 , the only possible values for three spins. Black points represent the extremal values of accessible filters, (1, 1) = Λ 0 , (−1, 1) = Λ 1 and (0, −1) = Λ 2 in the two dimensions, bounding a convex set denoted by the blue line. Only Λ 0 and Λ 2 are required to construct any other triangle wave where k ≥ 2, since all higher order waves are a weighted average of these two maps when d ≤ 2, and concatenation of sequences does not yield new filters: restricted to translation-invariant two-body spin Hamiltonians with on-site transverse fields. In Figure 2, we show C 3 and C 4 , including the universal set contained therein.

Efficient Resource Scaling of the Method
Proving access to a universal filter space still leaves open the question of how the actual control operations affect the scaling of the simulator's performance. We prove that it is possible to efficiently generate an arbitrary desired Ω eff (d) using linear combination and concatenation of the primary sets of filters we define. A key result of our method, described in detail in Appendix A, is that the realization of an arbitrary interaction Hamiltonian within this universal space may be accomplished using at most O(log N ) concatenation steps. These are the most resource-intensive operations we employ, requiring a pulse number of at most O(N 2 ) with at most O(N 3 ) individual spin flips. Arriving at this result is a complex exercise, and finding such an efficient resource scaling was not expected at the outset of this work.
Our ability to efficiently construct the universal filter space has significant impact on the runtime of filter implementation. The cost of universality is a rescaling of the interaction strength, reflected in the reduced volume of the hypercube relative to the extremal points of the convex set (see Fig. 2b-c). This reduced effective interaction strength within the universal filter space physically corresponds to an increase in the total evolution time required to achieve a target interaction form.
Naively, requiring O(log N ) layers of concatenation in realizing the universal filter space would result in a bound of N O(log N ) , which is only quasi-polynomial. To test this we calculate the quantity s(Q), the "worst-case" strength of the filtered interaction for qubits separated a distance Q > N/2, whose inverse sets the upper bound on runtime (see Appendix A for details of the construction of this equation). This expression is derived from the concatenated sequence described by Eq. A.5, and it is simply the product of the largest reduction in strength for each of the m concatenations. This is equivalent to assuming the smallest value of the interaction strength reduction for each of the layers of concatenation, thus providing a lower bound on the final interaction strength. We note that we can focus only on the last half of the chain without loss of generality because the filter strength typically decreases as a function of d.
In Fig. 2d we numerically evaluate this function for up to N = 2 20 , revealing exotic fractal-like features due to a recursive structure in the definition in terms of the binary digits of Q. More importantly, we find that the special structure of our pulse sequences gives a much stronger runtime upper bound with polynomial scaling, than that suggested by the naive expectation. Examining the results, we conjecture that a lower bound for this function is given by Q −2.585 (red line, Fig. 2d), the exponent given by − log 2 6. This is an exceptional observation, highlighting the practical relevance of this approach.
We may also consider the relevance of control errors in instances where complex concatenated sequences are employed. A chief source among these is the role of noninstantaneous control may be simply bounded by minimizing the ratio of τ π /T for a given sequence. Since we have demonstrated that the worst-case pulse number required to enact an arbitrary transformation scales polynomially with qubit number, we may therefore also bound the growth of this error. Error-resistant compensating pulses may also be employed, as has been discussed in the literature [37,38], or recent results assuming bounded-strength controls may be considered [39].

Efficient Linear Program for Finding Optimal Pulse Sequences
The existence of the universal filtering space guarantees that some pulse sequence exists in order to provide a desired transformation. The process of mapping a desired Hamiltonian for simulation to an appropriate and efficient "program" of filters is accomplished as a problem in linear programming. The algorithm presented here actually finds a valid pulse sequence (if one exists given the input filters) and minimizes the total amount of time spent applying pulses. Our results show that the problem of compiling pulse sequences can be solved efficiently in N , the number of spins, subject to some mild caveats. Additionally, many fast, practical implementations of linear programming algorithms exist to solve this problem at scales of interest for the implementation of useful quantum simulations in polynomial time.
Consider a set of elementary filter functions f 1 , . . . , f m , which we think of as vectors of dimension N − 1. These could come from e.g. the basic Λ and Γ filters defined above, together with k levels of concatenation for some fixed k. We can also add any other filters so long as we restrict to a polynomial upper bound on m, i.e. m = O(N c ) for some constant c independent of N , in order to ensure that the algorithm below runs in polynomial time. In particular, we can add the O(N ) specific "basis" we constructed in the Appendix at k = O(log N ) levels of concatenation which we use there to prove the universality of or scheme.
Given this set of dynamic filters, we wish to compile a positive linear combination of pulse sequences that will generate a given Hamiltonian. As mentioned above, this problem can be cast as a linear program, which is an efficient method for finding optima of linear objective functions subject to linear equality and inequality constraints [40]. Using standard linear algebra routines, this procedure may be implemented numerically and provides provably optimal solutions in terms of the total evolution time; the runtime of a linear program with n input variables and poly(n) constraints is polynomial in n.
To see that our problem is a linear program, we make the following observations. Given a desired vector of couplings Ω eff , we wish to find time steps t j ≥ 0 such that Here we actually have a vector equation, and the division on the righthand side is done elementwise. Suppose that each filter f j has a cost c j associated to it (for simplicity we can imagine that all of these costs are equal.) Then a linear program which will compile a given pulse sequence to generate the effective coupling Ω eff is The form of the program presented above may be converted to the "standard" form for linear programming by adding slack variables with additional equality constraints for each of the components of the vector equation. Finding a solution requires that the selected filter set is chosen to be universal. Regardless, the program will terminate in time polynomial in m and N and will minimize the total amount of evolution time required to implement the desired coupling.
Further improvements might be desired, such as having only a few nonzero values for t j (that is, requiring implementation of only a few different filters). By Carathéodory's theorem [41], only N of the filters need to be nonzero for any fixed value of Ω eff /Ω in order to guarantee a solution. However, finding such sparse solutions is in general NP-hard, so one would likely have to resort to heuristic methods to make improvements along these lines.

Example programs
Some examples of Hamiltonian engineering will help to reveal the power, generality, and utility of our approach. A simple example is the application of the filter f Λ k (d) + I , where I indicates free evolution of equal duration. This particular filter ensures that all qubits separated by distance k will be fully decoupled, combining a period of pulsed modulation with free-evolution, and exploiting the fact that f Λ k (k) = −1. Using concatenation, we can eliminate all interactions in the spin chain except for, e.g., nearestneighbor interactions. This Hamiltonian is vital not only for simulation of quantum magnetism, but also many other quantum information protocols [42,43].
A particularly useful example of a relevant Hamiltonian mapping relates to problems in quantum magnetism [44,45,46] where long-range qubit interactions can be engineered to scale as Ω d ∝ d −α , α ∈ [0, 3]: a form of interaction that arises in phononmediated spin simulators using trapped ions [47,48,35]. In practice, many simulators cannot reach the achievable limits of this range, or there may be a desire to induce a scaling outside of the range of this native interaction.
Here we provide details of how, using the pulse sequences we've introduced, one may adjust power-law interactions beyond the accessible bounds native to the underlying hardware. As an example, consider a system whose qubits interact through a Coulomb interaction such that Ω d ∝ Ω 1 /d. In order to use this system to simulate a many-body system whose particles interact through a 1/d 2 potential, we must construct a filter with functional dependence f (d) ∝ 1/d so that Ω d → Ω d /d. Such a filter may obviously be applied repeatedly to modify arbitrary power-law interactions.
For a system of N qubits, this transformation may be achieved through the sequential application of increasingly complex concatenations of Λ N +1 . Any Λ k , k > N can be used, but k = N + 1 minimizes the number of pulses in the primitive filter. In a system of N qubits, the map generated by Λ k is linear over the entire chain when k > N so that it can be described as highlights that this is equal to the linear term in the Taylor expansion of k/2 d around d = k/2, with a radius of convergence of k/2; We define the evolution operator U 1/d n which is generated by an interaction that obeys a power law of Ω d = Ω 1 /d n . Consider the following filtered evolution operator V , Here the exponent on the symbol Λ x k denotes concatenation to x levels, and the last line follows by summing the geometric series. The effective interaction strength now obeys a new inverse power law where the exponent has been incremented by a single power. More generally, one may choose different timings for the various filtering operations in order to construct a more general dynamical mapping function, The function g(d, {α i }) is restricted by the fact that the coefficients {α i } are necessarily positive. This implies that the function f must have Taylor expansion coefficients that alternate in sign, meaning that the technique is capable of producing any inverse power law filter but is incapable of producing polynomials with positive exponents. Another possibility to realize these power laws is to use the linear program of the previous section together with a family of pulse sequences to try to find the most time-efficient power law across all N of the spins.

Application to Adiabatic Protocols
Next, we describe how this approach may be applied to adiabatic simulators and preparation of entangled ground-states of designer Hamiltonians [49]. Say we wish to adiabatically evolve to the ground state of a target Hamiltonian, H t , but can only turn on the available Hamiltonian H a in a particular experimental apparatus. Returning again to our 1D qubit chain, we initialize in |g s , which is the ground state of a simple Hamiltonian such as H s = − ω 2 j Z j , and allow the system to evolve under the timedependent Hamiltonian, H(t) = 1 − t τ H s + t τ H a . If the evolution is adiabatic, then | g a |Û a (τ )|g s | 2 ≈ 1 at the end of the interaction. Now we apply stroboscopic filtering to drive the system to the new ground state of H t , |g t . The filtering must be done repeatedly throughout the adiabatic evolution since the Hamiltonian generally does not commute with itself at different times. The rate of application is set such that, to first order, the Hamiltonian is constant over the pulse period, ∆t, giving a total evolution defined by R l=1 F Û a (l∆t, (l − 1)∆t) with R = τ /∆t the number of filtering operations performed during the adiabatic ramp, F the dynamic filter, andÛ a (t 1 , t 2 ) the time-ordered evolution from t 1 to t 2 .
Using the filtering approach in adiabatic simulators to reach a target ground state |g t , requires repeatedly applying the sequence at a rate fast enough that the evolving Hamiltonian is approximately constant over the time it takes to apply the pulse sequence. Defining a filter F that maps the evolution operator generated by H a to one generated by H t gives a piecewise-constant filtered adiabatic evolution operator where R = τ /∆t is the number of filtering operations performed during the adiabatic ramp. The state then evolves during these discrete timesteps as Here H (j) f is the filtered Hamiltonian at t = j∆t. The error accrued due to the Hamiltonian changing during the filtering operations can be calculated as the overlap of |ψ with a state |φ whose evolution includes the firstorder time-dependence of the Hamiltonian during filtering operations. In the adiabatic approximation, the time evolution operator is given by exp[−i t+∆t t H(t )dt ] which can be approximated by exp[−i(H(t)∆t + 1 2Ḣ (t)∆t 2 )] when ∆t is small. In cases like ours where,P j d dt H(t)P j = d dtP j H(t)P j we can write which is meant to capture the detrimental effect on the filtering operation due to the changing Hamiltonian to lowest order in ∆t. The fidelity of the filtered adiabatic protocol can be defined as the overlap of |ψ and |φ as defined in Eqs. (20) and (21). Thus, we have by the Baker-Campbell-Hausdorff formula. The sum R j=1Ḣ (j) ∆t can be approximated by an integral when ∆t is small enough so that the overlap is, From this last expression, we see that the error can be thought of as the probability of the Hamiltonian 1 2 (H(t f ) − H(t 1 )) evolving the system out of the state |g s in a time ∆t. To lowest order in ∆t, this error is Var(H(t f )−H(t 1 ))∆t 2 /4 2 where the variance is calculated with respect to |g s . In the adiabatic evolution described in this article, |g s is the ground state of H(t 1 ), meaning that the expression for the error can be further simplified to Var(H(t f ))∆t 2 /4 2 .
To test this, we numerically integrate the Schrödinger equation for a system of 4 qubits undergoing adiabatic evolution, starting from a distance-independent spin coupling, Ω d ≡ Ω, and filtered to produce only nearest-neighbor coupling. We calculate the error accrued due to deviation from the assumption of a piecewise-constant Hamiltonian during the filtering operations as the overlap of the ideal state with a state whose evolution includes a first-order time-dependence of the Hamiltonian during filtering and find that the infidelity decreases approximately as (∆t) 2 . In this case, the overlap between the ground state of the target Hamiltonian and that of the unfiltered Hamiltonian is | g t |g a | 2 = 0.33, but with the application of dynamic filters, surpasses 10 −2 after just six filtering steps.

Moving Beyond Ising Interactions
We now relax the original constraints we placed on our interaction Hamiltonian in order to simplify the analytic treatments presented above in demonstrating the basic functionality of our approach. All of these results continue to hold when our Ising-type model is generalised to include both the homogeneous on-site terms and the presence of long-range Heisenberg-type interactions.

Heisenberg Interactions
Similar to the work of Hodges et al. [50], an Ising interaction (e.g. XX) can be modified into a Heisenberg-type (S · S) by including global basis-changing π/2 pulses in the sequence of applied filters. Similarly, any native Heisenberg interaction may be modified into any other using pulse sequences applied in the appropriate basis and concatenated We may rigorously bound the error and runtime as a polynomial function of T and N , finding only a constant multiplicative penalty in run-time and Trotter error (see below in Section 7.2). Suppose that our target Hamiltonian is of the form (cf. Eq. 1) For example, if Ω v j = 0 for v = x, y, z and j ≥ 2, then this is the standard spin-1/2 XYZ model. This Hamiltonian can be engineered through a simple extension of the method we describe above. Assuming the native interaction Hamiltonian consist solely of X j operators, we build up the time evolution operator using another layer of Trotterization.
To be more precise, we generate the Y j interactions by applying global π/2 pulses about the Z axis, and the Z j interactions are generated by global π/2 pulses about the Y axis. If the filtering sequence is represented by the superoperator F, and the global π/2 pulses are represented by G, the Heisenberg evolution U H is given by, While this method necessarily incurs a Trotter-type error in the implementation, we demonstrate below that this error can be bounded, and provide such a bound for the simulation error (that is, for the case of perfect gates).

Inclusion of a Transverse Field
Previously, in the case of the Ising coupling, the transverse field was set to zero so that all of the terms in the Hamiltonian would commute with one another. This resulted in an exact mapping from an initial Hamiltonian of the form (28) to another Hamiltonian of the same form, but with the distance function being modified The functional dependence f (d) can take an arbitrary form, up to an overall rescaling factor.
To incorporate a homogeneous transverse-field term of the form we use a Suzuki-Trotter-type decomposition of the evolution operator [51,52,53]. We assume that the overall strength B of the transverse field is a tunable parameter so that we can match its rescaled strength relative to the rescaled coupling strengthΩ d by the same overall scale factor. Given this assumption, there are two natural ways to analyze this situation: an always-on field and a field which can be switched on and off rapidly.
In some implementations, it is further reasonable to expect that the interactions can be switched on and off rapidly (e.g., some trapped ion experiments [54]). We will also analyze this case, since in fact it is the easiest to analyze and gives insight into the other results.

Bounding the Trotter error in a simulation
We wish to bound the deviation between the exact and the approximate (Suzuki-Trotter) evolution, assuming that all of our pulse sequences are accurate and perfectly timed. We denote the ideal evolution operator by where jP j H IPj ∆t =H I t. The approximate evolution is defined in terms of the second-order integrator (the "split-step" method), given by the following formula: Here the H j are the allowed evolutions over a given time interval, and their exact definition depends on which resources we allow ourselves in the protocol, but they are related to the ideal effective Hamiltonian (here called simply H) by We will discuss the specific choices below. We have also introduced a parameter r, the Trotter number, which counts the number of steps in our Suzuki-Trotter expansion and hence controls our error.
Note that we cannot use the higher-order integrators which provide faster convergence because they all require evolution for negative values of t [53]. If rapidly flipping the sign of the interactions and the transverse field is feasible in a given system, then these higher-order integrators could be used to give sharper error bounds than the ones given below.
Given this very general framework, we can upper bound the error as follows [55] Here the error is in terms of the operator norm (largest singular value). Note that H ≤ N 2 /2 + BN by the triangle inequality, so this quantity is always polynomial in the number of spins. This bound on the norm of the difference has an operational meaning: it gives an upper bound on the trace distance between the evolved states, namely for any initial state ρ and unitaries U and V .
To make the error bound in Eq. (33) more concrete, we need to determine m for a given protocol. The form of the bound above indicates that a smaller number of different Hamiltonians in use results in improved error scaling.
The simplest case to consider is that where both the interactions and the transverse field can be switched on and off rapidly. In this case we can use the above bound with m = 2, in which we simply switch between (1) the exact evolution of the effective Hamiltonian with zero transverse field, and (2) the transverse field alone with no interactions. That is, we have Note thatH I can be implemented exactly by subdividing each Suzuki-Trotter step into smaller segments, as discussed in the main text. Next, consider the case where the interactions are always on, but the transverse field can be switched on and off rapidly. Here a good strategy is to turn on the transverse field in short, strong bursts, and turn it off while applying the filter functions. The control Hamiltonians during each of the m = 2 Suzuki-Trotter steps are given by Here we introduce an additional parameter δ 1, and we incur a small error in the final Hamiltonian, We can choose δ as small as possible to be compatible with the assumption of accurate and rapid switching of the transverse field. Thus, the error bound is nearly the same as the previous m = 2 bound, except that we add a small error from the inexact decomposition.
Next, consider the most pessimistic case where both the interactions and the transverse field are not rapidly switchable. In this case the best strategy seems to be to use a different Suzuki-Trotter step for each increment in a pulse sequence. In this case, we have no choice but to use the bound in Eq. (33) by setting m equal to the total number of pulses in the pulse sequence. For a chain of length N , we can always bound the maximum number of pulses required by m ≤ N , which follows from Carathéodory's theorem [41]. This gives a pessimistic but nonetheless polynomial bound on the required Trotter number to achieve a fixed error.
Finally, we remark that the same methods for bounding the Trotter error in this section also apply in Section 7. Since this more general case includes, at most, three times the number of pairwise coupling terms in the Hamiltonian, this implies that the error remains polynomial in the number of qubits.

Conclusion
In conclusion, we have provided a general framework for dynamic filtering of spin Hamiltonians that enables programmable quantum simulators using single-qubit Pauli operations and a native long-range spin coupling. The scheme we introduce provides a system designer with the ability to perform numerical decomposition of a large class of realizable couplings into the basis of available filters, thus providing a prescription of how to "program" the desired interactions. We have further showed how this technique can be applied to augment the adiabatic evolution of a spin Hamiltonian in a form of hybrid adiabatic quantum simulator.
We are excited by the fact that our approach to combine such filters in realizing a desired program for quantum simulation represents a contribution to a growing body of literature focused on the benefits of functional analysis in the context of protocol, algorithm, and circuit development for quantum information. In particular, following completion of this work we discovered significant similarity between our approach and the synthesis of diagonal operators for circuit-model quantum computation [24]. Given these connections and the breadth of applicability of functional analysis in quantum information, we hope that our results enable new capabilities well beyond quantum simulation.
The results presented herein demonstrate the utility of the pulse sequences we define for efficiently modifying native Hamiltonians in quantum simulators. From the perspective of control engineering these results demonstrate that via the quantum control protocols we have defined, the 1D translationally invariant quantum simulator in which we are interested is efficiently controllable, requiring resources scaling only polynomially in time and energy (using pulse number as a proxy). Assuming only Isingtype long range interactions, we may realise a much broader class of Hamiltonians including spin-1/2 XY Z Heisenberg-type models with homogeneous on-stie terms and much more exotic multi-axis spin-coupling Hamiltonians with bounded error and polynomial resources. All of our pulse sequences and methods of filter combination may be applied in order to achieve these unusual models with the small associated penalties in time and resources.
While we have provided explicit examples of how one might tune the power-law of a long-range spin coupling or completely cancel undesired spin couplings on a 1D lattice, it is an interesting open problem to find additional examples where our method can be applied efficiently. This includes, for instance, the efficient generation of an effective Heisenberg Hamiltonian using only Ising couplings and single-qubit control such as that in Ref. [50], or generalizations of our method to systems with a natural 2D geometry of spins. A fundamental open problem is also finding optimal primary filter constructions that minimize the overhead in evolution time for a desired effective Hamiltonian. In addition we wish to explore the possibility of finding a rigorous proof of polynomial runtime scaling, rather than relying only on numeric evidence supporting our claim. Finally, it is an exciting open question to explore how these techniques may be coupled with notions of fault-tolerance in quantum simulation. We note that after completion of this work we became aware of the use of similar language of Hamiltonian "filtering" for quantum simulation, treating quantum simulation and spin transport in a related framework [56].