Hamiltonian quantum simulation with bounded-strength controls

We propose dynamical control schemes for Hamiltonian simulation in many-body quantum systems that avoid instantaneous control operations and rely solely on realistic bounded-strength control Hamiltonians. Each simulation protocol consists of periodic repetitions of a basic control block, constructed as a suitable modification of an"Eulerian decoupling cycle,"that would otherwise implement a trivial (zero) target Hamiltonian. For an open quantum system coupled to an uncontrollable environment, our approach may be employed to engineer an effective evolution that simulates a target Hamiltonian on the system, while suppressing unwanted decoherence to the leading order. We present illustrative applications to both closed- and open-system simulation settings, with emphasis on simulation of non-local (two-body) Hamiltonians using only local (one-body) controls. In particular, we provide simulation schemes applicable to Heisenberg-coupled spin chains exposed to general linear decoherence, and show how to simulate Kitaev's honeycomb lattice Hamiltonian starting from Ising-coupled qubits, as potentially relevant to the dynamical generation of a topologically protected quantum memory. Additional implications for quantum information processing are discussed.


Introduction
The ability to accurately engineer the Hamiltonian of complex quantum systems is both a fundamental control task and a prerequisite for quantum simulation, as originally envisioned by Feynman [1,2,3].The basic idea underlying Hamiltonian simulation is to use an available quantum system, together with available (classical or quantum) control resources, to emulate the dynamical evolution that would have occurred under a different, desired Hamiltonian not directly accessible to implementation [4].From a control-theory standpoint, the simplest setting is provided by open-loop Hamiltonian engineering in the time domain [5,6], whereby coherent control over the system of interest is achieved solely based on suitably designed time-dependent modulation (most commonly sequences of control pulses), without access to ancillary quantum resources and/or measurement and feedback.While open-loop Hamiltonian engineering techniques have their origin and a long tradition in nuclear magnetic resonance (NMR) [8,7], the underlying physical principles of "coherent averaging" have recently found widespread use in the context of quantum information processing (QIP), leading in particular to dynamical symmetrization and dynamical decoupling (DD) schemes for control and decoherence suppression in open quantum systems [9,10,11,12,13,14].
As applications for quantum simulators continue to emerge across a vast array of problems in physics and chemistry, and implementations become closer to experimental reality [3,15,16], it becomes imperative to expand the repertoire of available Hamiltonian simulation procedures, while scrutinizing the validity of the relevant control assumptions.With a few exceptions (notably, the use of so-called "perturbation theory gadgets" [17]), open-loop Hamiltonian simulation schemes have largely relied thus far on the ability to implement sequences of effectively instantaneous, "bangbang" (BB) control pulses [18,19,20,21,22,23,24,25]. While this is a convenient and often reasonable first approximation, instantaneous pulses necessarily involve unbounded control amplitude and/or power, something which is out of reach for many control devices of interest and is fundamentally unphysical.In the context of DD, a general approach for achieving (to at least the leading order) the same dynamical symmetrization as in the BB limit was proposed in [26], based on the idea of continuously applying bounded-strength control Hamiltonians according to an Eulerian cycle, so-called Eulerian DD (EDD).From a Hamiltonian engineering perspective, EDD protocols translate directly into bounded-strength simulation schemes for specific effective Hamiltonians -most commonly, the trivial (zero) Hamiltonian in the case of "non-selective averaging" for quantum memory (or "time-suspension" in NMR terminology).More recently, EDD has also served as the starting point for boundedstrength gate simulation schemes in the presence of decoherence, so-called dynamically corrected gates (DCGs) for universal quantum computation [27,28,29,30].
In this work, we show that the approach of Eulerian control can be further systematically exploited to construct bounded-strength Hamiltonian simulation schemes for a broad class of target evolutions on both closed and open (finite-dimensional) quantum systems.Our techniques are device-independent and broadly applicable, thus substantially expanding the control toolbox for programming complex Hamiltonians into existing or near-term quantum simulators subject to realistic control assumptions.
The content is organized as follows.We begin in Sect.II by introducing the appropriate control-theoretic framework and by reviewing the basic principles underlying open-loop simulation via average Hamiltonian theory, along with its application to Hamiltonian simulation in the BB setting.Sect.III is devoted to constructing and analyzing simulation schemes that employ bounded-strength controls: while Sec.III.A reviews required background material on EDD, Sec.III.B introduces our new Eulerian simulation protocols for a generic closed quantum system.In Sec.
III.C we separately address the important problem of Hamiltonian simulation in the presence of slowly-correlated (non-Markovian) decoherence, identifying conditions under which a desired Hamiltonian may be enacted on the target system while simultaneously decoupling the latter from its environment, and making further contact with DCG protocols.Sect.IV presents a number of illustrative applications of our general simulation schemes in interacting multi-qubit networks.In particular, we provide explicit protocols to simulate a large family of two-body Hamiltonians in Heisenbergcoupled spin systems additionally exposed to depolarization or dephasing, as well as to achieve Kitaev's honeycomb lattice Hamiltonian starting from Ising-coupled qubits.In all cases, only local (single-qubit, possibly collective) control Hamiltonians with bounded strength are employed.A brief summary and outlook conclude in Sec.V.

Control-theoretic framework
We consider a quantum system S, with associated Hilbert space H, whose evolution is described by a time-independent Hamiltonian H.As mentioned, Hamiltonian simulation is the task of making S evolve under some other time-independent target Hamiltonian, say, H. Without loss of generality, both the input and the target Hamiltonians may be taken to be traceless.Two related scenarios are worth distinguishing for QIP purposes: • Closed-system simulation, in which case S coincides with the quantum system of interest, S (also referred to as the "target" henceforth), which undergoes purely unitary (coherent) dynamics; • Open-system simulation, in which case S is a bipartite system on H ≡ H S ⊗ H B , where B represents an uncontrollable environment (also referred to as bath henceforth), and the reduced dynamics of the target system S is non-unitary in general.
In both cases, we shall assume the target system S to be a network of interacting qudits, hence H S (C d ) ⊗n , for finite d and n.In the general open-system scenario, the joint Hamiltonian on H may be expressed in the following form, where the operators H S (H B ) and S α (B α ) act on H S (H B ) respectively, and all the bath operators are assumed to be norm-bounded, but otherwise unspecified (potentially unknown).The closed-system setting is recovered from Eq. (1) in the limit S α ≡ 0. Likewise, we may express the target Hamiltonian H in a similar form, with two simulation tasks being of special relevance: Sα ≡ 0, in which case the objective is to realize a desired system Hamiltonian HS while dynamically decoupling S from its bath B, thereby suppressing unwanted decoherence [11]; or, more generally, H S → HS and S α → Sα , where the simulated, dynamically symmetrized error generators Sα may allow for decoherence-free subspaces or subsystems to exist [13,31].
Formally, the dynamics is modified by an open-loop controller acting on the target system according to where the operators {X u = X † u } and the (real) functions {f u (t)} represent the available control Hamiltonians and the corresponding, generally time-dependent, control inputs respectively.Clearly, if the Hamiltonian ( H − H) is contained in the admissible control set, the corresponding control problem is trivial and the desired time-evolution, Ũ (t) = e −i Ht , t ≥ 0, can be exactly simulated continuously in time.However, this level of control need not be available in settings of interest, including open quantum systems where control actions are necessarily restricted to the target system S alone, H c (t) ≡ H c (t) ⊗ I B in Eq. ( 2).Following the general idea of "analog" quantum simulation [3], we shall assume in what follows a restricted set of control Hamiltonians (in a sense to be made more precise later) and focus on the task of approximately simulating the desired time evolution Ũ (t) at a final time t = Tf , or more generally, stroboscopically in time, that is, at instants t = tM , where and T is a fixed minimum time interval.Choosing T sufficiently small allows in principle any desired accuracy in the approximation to be met, with the limit T → 0 formally recovering the continuous limit.Specifically, let U (t) and U c (t) denote the unitary propagators associated to the total and the control Hamiltonians in Eq. ( 2), respectively: where we have set = 1 and T indicates time-ordering, as usual.Then, for a given pair (H, H), we shall provide sufficient conditions for H to be "reachable" from H and, if so, devise a cyclic control procedure such that the resulting controlled propagator where T c is the cycle time of the controller, that is, U c (t + T c ) = U c (t).In general, we shall allow for T c to differ from T , corresponding to an overall scale factor in the simulated time, as it will become apparent later.If, for a fixed input Hamiltonian H, arbitrary target Hamiltonians are reachable for given control resources, the simulation scheme is referred to as universal.In this case, complete controllability must be ensured by the tunable Hamiltonians X u in conjunction with the "drift" H S [6].In contrast, we shall be especially interested in situations where control over S is more limited.
Similar to DD protocols, Hamiltonian simulation protocols are most easily constructed and analyzed by effecting a transformation to the "toggling" frame associated to U c (t) in Eq. (4) [7,11,14].That is, evolution in the toggling frame is generated by the time-dependent, control-modulated Hamiltonian with the corresponding toggling-frame propagator U (t) being related to the physical propagator in Eq. ( 3) by U (t) = U c (t)U (t).Since the control propagator is cyclic and H is time-independent, it follows that U (t M ) = U (t M ) and, furthermore, H (t) acquires the periodicity of the controller, U (t M ) = [U (T c )] M .Thus, the stroboscopic controlled dynamics of the system is determined by Average Hamiltonian theory [7,35] may then be invoked to associate an effective timeindependent Hamiltonian H to the evolution in the toggling-frame: where H is determined by the Magnus expansion [32], H = H(0) + H(1) + H(2) + . . .Explicitly, the leading-order term, determining evolution over a cycle up to the first order in time, is given by with (absolute) convergence being ensured as long as t H < π [34].Subject to convergence condition, higher-order corrections for evolution over time t can also be upper-bounded by (see Lemma 4 in [33]) Ideally, one would like to achieve HT c = H T , so that equality would hold in Eq. ( 5) for all M ∈ N. In what follow, we shall primarily focus on achieving first-order simulation instead, by engineering the control propagator U c (t) in such a way that whereby, using Eq. ( 10) with κ = 1, In general, the accuracy of the approximation in Eq. ( 11) improves as the "fast control limit", T c → 0, is approached.Physically, this corresponds to requiring that the shortest control time scale (pulse separation) involved in the control sequence be sufficiently small relative to the shortest correlation time of the dynamics induced by H [35,36].
While the problem of constructing general high-order Hamiltonian simulation schemes is of separate interest, we stress that second-order simulation can be readily achieved, in principle, by ensuring that Since all odd-order Magnus corrections vanish in this case [35], it follows (by using again Eq. (10), with κ = 2), that HT c = H T + O[( H T c ) 3 ], correspondingly boosting the accuracy of the simulation.

Hamiltonian simulation with bang-bang controls
BB Hamiltonian simulation provides the simplest control setting for achieving the intended objective, given in Eq. (5).Two main assumptions are involved: (i) First, we must be able to express the target Hamiltonian H as where {U j } are unitary operators on S and the {w j } non-negative real numbers (not all zero).(ii) Second, the available control resources include a discrete set of instantaneous pulses {P j } on S, whose application results in a piecewise-constant control propagator U c (t) over [0, T c ], with corresponding toggling-frame propagators {U j }, U j ≡ j k=1 P k , U 1 = I S [9,14].Assumptions (i)-(ii) together allow for the time-average in Eq. ( 9) to be mapped to a convex (positive-weighted) sum.Eq. ( 13) may be interpreted as a sufficient condition for the target Hamiltonian H to be reachable from H given open-loop unitary control on S alone.Reachable Hamiltonians must thus be at least as "disordered" as the input one in the sense of majorization [21,14].Specifically, Eq. ( 13) leads naturally to the following BB simulation scheme.Given simulation weights {w j }, define the following simulation intervals and timing pattern: A piecewise-constant control propagator for the basic simulation block to be repeated may then be constructed as follows: By using Eq. ( 9), it is immediate to verify that resulting in the desired controlled evolution, Eqs. ( 11)- (12), provided that the convergence conditions for first-order simulation under H are obeyed.Since, in practice, technological limitations always constrain the cycle duration to a finite minimum value T c > 0, such conditions ultimately determine the maximum simulated time tM up to which evolution under H may be reliably simulated using the physical Hamiltonian H.
In analogy with BB DD schemes, realizing the prescription of Eq. ( 15) requires to discontinuously change the control propagator from U j to U j+1 = (U j+1 U † j )U j , via an instantaneous BB pulse U j+1 U † j at the jth endpoint t j .As a result, despite its conceptual simplicity, BB simulation is unrealistic whenever large control amplitudes are not an option, and the evolution induced by H during the application of a control pulse must be considered from the outset.This demands redesigning the basic control block in such a way that the actions of H and H c (t) are simultaneously accounted for.

Eulerian simulation of the trivial Hamiltonian
The key to overcome the disadvantages of BB Hamiltonian simulation is to ensure that the control propagator varies smoothly (continuously) in time during each control cycle.We achieve this goal by relying on Eulerian control design [26].To introduce the necessary group-theoretical background, we begin by revisiting how, for the special case of a target identity evolution (that is, H ≡ 0, also corresponding to a "noop" gate, in terms of the end-time simulated propagator), EDD can be naturally interpreted as a bounded-strength simulation scheme.
In the Eulerian approach, the available control resources include a discrete set of unitary operations on S, say, {U γ }, γ = 1, . . ., L, which are realized over a finite time interval ∆ through application of bounded-strength control Hamiltonians Note that the choice of the control Hamiltonians h γ (t) is not unique, which allows for implementation flexibility.The unitaries {U γ } are identified with the image of a generating set of a finite group under a faithful, unitary, projective representation ρ [26].That is, let G ≡ {g} be a finite group of order |G|, such that each element may be written as an ordered product of elements in a generating set Γ ≡ {γ} of order |Γ| = L, g → ρ(g) ≡ U g be the representation map [37], and G ≡ {U g }.The Cayley graph C(G, Γ) of G relative to Γ can be thought of as pictorially representing all elements of G as strings of generators in Γ.Each vertex represents a group element and a vertex g is connected to another vertex g by a directed edge "colored" (labeled) with generator γ if and only if g = γg.The number of edges in C(G, Γ) is thus equal to N ≡ |Γ||G|.
Because a Cayley graph is regular, it always has an Eulerian cycle that visits each edge exactly once and starts (and ends) on the same vertex [38,39].Let us denote with C ≡ (γ 1 , . . ., γ N ) the ordered list of generators defining an Eulerian cycle on C(G, Γ) which, without loss of generality, starts (and ends) at the identity element of G.
Once a control Hamiltonian for implementing each generator as in Eq. ( 17) is chosen, an EDD protocol is constructed by assigning a cycle time as T c ≡ N ∆ and by applying the control Hamiltonians h γ (t) sequentially in time, following the order determined by the Eulerian cycle C.Thus, where U γ j is the image of the generator labeling the jth edge in C. As established in [26], the lowest-order average Hamiltonian associated to the above EDD cycle has the form H(0 , where for any operator A acting on H S , the map projects onto the centralizer of G (i.e., Π G (A) commutes with all U g ∈ G), and implements an average of H over both the control interval and the group generators.Accordingly, bounded-strength simulation of H = 0 is achieved provided that the following DD condition is obeyed: By Schur's lemma, this is automatically ensured if the group representation acts irreducibly on H S .Formally, the BB limit may be recovered by letting F Γ (A) ≡ A for all A [26], reflecting the ability to directly implement all the group elements (with no overhead, as if |Γ| = 1) and corresponding to uniform simulation weights w j = 1/|G|.

Eulerian simulation protocols beyond noop: Construction
We show how the Eulerian cycle method can be extended to bounded-strength simulation of a non-trivial class of target Hamiltonians.We assume that H may be expressed as a convex unitary mixture of the group representatives U g , We construct the desired control protocol starting from an Eulerian cycle C = (γ 1 , . . ., γ N ) on C(G, Γ).Specifically, the idea is to append to each of the N control slots that define an EDD scheme a free-evolution (or "coasting") period of suitable duration, in such a way that the net simulated Hamiltonian is modified from H = 0 to H = 0 as given in Eq. (22).A pictorial representation of the basic control block is given in Fig. 1.As in Eq. ( 17), let ∆ denote the minimum time duration required to implement each generator, hence, to smoothly change the control propagator from a value U g to U g along the cycle.While such "ramping up" control intervals have all the same length, each "coasting" interval is designed to keep the control propagator constant at U g for a duration determined by the corresponding weight w g .Since the control is switched off during coasting, continuity of the overall control Hamiltonian H c (t) may be ensured, if desired, by requiring that in addition to the bounded-strength constraint.An Eulerian simulation protocol may be formally specified as follows.As before, let the jth time interval be denoted as [t j−1 , t j ], j = 1, . . ., N , with t 0 = 0 and t N defining Schematics of an Eulerian simulation protocol.The basic control block consists of N time intervals, each involving a "ramping-up" sub-interval of fixed duration ∆, during which H c (t) = 0, followed by a "coasting" (free evolution) period of variable duration Θ k , Eq. ( 24), during which no control is applied.During the jth ramping-up sub-interval we apply h γj , i.e., the control Hamiltonian that realizes the generator γ j , smoothly changing the control propagator from U gj−1 to U gj .In this way, the control protocol corresponding to Eqs. ( 26)-( 27) is implemented.By construction, a standard EDD cycle with H = 0 is recovered by letting Θ k → 0 for all k, while in the limit ∆ → 0 standard BB simulation of H is implemented.
the cycle time T c .For each j, let τ g j ≡ w g j T as in the BB case.The duration of the jth coasting period is then assigned as resulting in the following timing pattern {t j } [compare to Eq. ( 14)]: As the expression for the cycle times makes it clear, the resulting protocol may be equivalently interpreted in two ways: starting from an EDD cycle, corresponding to N ∆ and H = 0, we introduce the coasting periods to allow for non-trivial simulated dynamics to emerge; or, starting from a BB simulation scheme for H, corresponding to W T , we introduce the ramping-up periods to allow for control Hamiltonians to be smoothly switched over ∆.Either way, bounded-strength protocols imply a time-overhead N ∆ relative to the BB case, recovering the BB limit as ∆ → 0 as expected.Explicitly, the control propagator for Eulerian simulation has the form: The resulting first-order Hamiltonian H(0) under Eulerian simulation is derived by evaluating the time-average in Eq. ( 9) with the control propagator given by Eqs. ( 26)- (27).Direct calculation along the lines of [26] yields: where the last equality follows from two basic properties of Eulerian cycles: firstly, the list {g 0 , g 1 , . . ., g N −1 } (and also {g 1 , g 2 , . . ., g N }) of the vertices that are being visited contains each element g ∈ G precisely |Γ| times; secondly, in traversing the Cayley graph, each group element g is left exactly once by a γ-labeled edge for each generator γ ∈ Γ.Thus, by recalling the definitions given in Eqs. ( 19) and ( 20), we finally obtain which indeed achieves the intended first-order simulation goal, Eqs. ( 11)-( 12), as long as convergence holds and the DD condition of Eq. ( 21) is obeyed.
The simulation accuracy may be improved by symmetrizing U EUS c (t) in time.In analogy to symmetrized EDD protocols [9], this can be easily accomplished by running the protocol and then suitably running it again in reverse.Specifically, let the duration of the coasting interval be changed as Θ j → Θ j /2.Run the protocol as described above until time t = N ∆ + 1 2 W T .Then, from time t = N ∆ + 1 2 W T until time t = T c = 2N ∆ + W T , modify Eqs. ( 26)- (27) as follows: for j = N, . . ., 1. Provided that one is able to implement u γ j (∆ − δ), we again obtain hence ensuring that H(1) = 0.

Eulerian simulation while decoupling from an environment
The ability to implement a desired Hamiltonian on the target system S, while switching off (at least to the leading order) the coupling to an uncontrollable environment B, is highly relevant to realistic applications.That is, with reference to Eq. ( 1), the objective is now to simultaneously achieve HS ≡ H target , Sα ≡ 0, by unitary control operations acting on S alone.Because the first-order Magnus term H(0) is additive [recall Eq. ( 9)], it is appropriate to treat each summand of H individually, leading to a relevant average Hamiltonian of the form where for a generic operator on H S we let We can then apply the analysis of Sec.3.2 to the internal system Hamiltonian ( HS ) and each error generator ( Sα ) separately, to obtain in both cases a simulated operator of the form given in Eq. ( 28): Since the task is to decouple S from B while maintaining the non-trivial evolution due to HS = H target , the reachability condition of Eq. ( 22) must now ensure that Accordingly, it is necessary to extend the DD assumption of Eq. ( 21) to become such that Ā = ( T /T c ) Ã holds for each of the summands in H. Altogether we recover It is interesting in this context to highlight some similarities and differences with DCGs [27], which also use Eulerian control as their starting point and are specifically designed to achieve a desired unitary evolution on the target system while simultaneously removing decoherence to the leading [27,28,30] or, in principle, arbitrarily high order [29].By construction, the open-system simulation procedure just described does provide a first-order DCG implementation for the target gate Q ≡ exp(−i HS Tf ): in particular, the requirement that Eqs. ( 29)- (30) be obeyed together (for the same weights w g ) is effectively equivalent to evading the "no-go theorem" for black-box DCG constructions established in [28], with the coasting intervals and the resulting "augmented" Cayley graph playing a role similar in spirit to a (first-order) "balancepair" implementation.Despite these formal similarities, a number of differences exist between the two approaches: first, an obvious yet important difference is that DCG constructions focus directly on synthesizing a desired unitary propagator, as opposed to a desired Hamiltonian generator.Second, while the internal system Hamiltonian, H S , is a crucial input in a Hamiltonian simulation problem, it is effectively treated as an unwanted error contribution in analytical DCG constructions, in which case complete controllability over the target system must be supplied by the controls alone.Although in more general (optimal-control inspired) DCG constructions [30], limited external control is assumed and H S may become essential for universality to be maintained, emphasis remains, as noted above, on end-time synthesis of a target propagator.Finally, a main intended application of DCGs is realizing low-error single-and two-qubit gates for use within fault-tolerant quantum computing architectures, as opposed to robust Hamiltonian engineering for many-body quantum simulators which is our focus here.

Eulerian simulation protocols: Requirements
Before presenting explicit illustrative applications, we summarize and critically assess the various requirements that should be obeyed for Eulerian simulation to achieve the intended control objective of Eq. ( 5) in a closed or, respectively, open-system setting: (i) Time independence.Both the internal Hamiltonian H and the target Hamiltonian H are taken to be time-independent (and, without loss of generality, traceless).
(ii) Reachability.The target Hamiltonian H must be reachable from H, that is, there must be a control group G, with a faithful, unitary projective representation mapping g → ρ(g) = U g , such that Eq. ( 22) holds.For dynamically-corrected Eulerian simulation in the presence of an environment, this requires, as noted, that for the same weights {w g }, the desired system Hamiltonian is reachable from H S while the trivial (zero) Hamiltonian is reachable from each error generator S α separately, such that both Eqs. ( 29)-( 30) hold together.
(iii) Bounded control.For each generator γ of the chosen control group G, we need access to bounded control Hamiltonians h γ (t), such that application of h γ (t) over a time interval of duration ∆ realizes the group representative U γ = ρ(γ) = u γ (∆), additionally subject (if desired) to the continuity condition of Eq. ( 23).
(iv) Decoupling conditions.Suitable DD conditions, Eq. ( 21) in a closed system or Eqs. ( 31)- (32) in the open-system error-corrected case, must be fulfilled, in order for undesired contributions to the simulated Hamiltonians to be averaged out by symmetry to the leading order.
(v) Time-efficiency.If the choice of G is not unique for given (H, H), the smallest group should be chosen, in order to keep the number of intervals per cycle, N = |G||Γ|, to a minimum.In particular, efficient Hamiltonian simulation requires that |G| (hence also |Γ|) scales (at most) polynomially with the number of subsystems n.
The key simplification that the time-independence Assumption (i) introduces into the problem is that the periodicity of the control action is directly transferred to the toggling-frame Hamiltonian of Eq. ( 6), allowing one to simply focus on single-cycle evolution.Although this assumption is strictly not fundamental, general time-dependent Hamiltonians may need to be dealt with on a case-by-case basis (see also [40,41,42]).A situation of special practical relevance arises in this context for open systems exposed to classical noise, in which case H B C and the system-bath interaction in Eq. ( 1) is effectively replaced by a classical, time-dependent stochastic field.Similar to DD and DCG schemes, Eulerian simulation protocols remain applicable as long as the noise process is stationary and exhibiting correlations over sufficiently long time scales [9,43].
The reachability Assumption (ii) is a prerequisite for Eulerian Hamiltonian simulation schemes.Although BB Hamiltonian simulation need not be group-based, most BB schemes follow this design principle alike.Assumption (iii), restricting the admissible control resources to physical Hamiltonians with bounded amplitude (thus finite control durations, as opposed to instantaneous implementation of arbitrary group unitaries as in the BB case) is a basic assumption of the Eulerian control approach.As remarked, our premise is that the available Hamiltonian control is limited, restricted to only the target system if the latter is coupled to an environment, and typically non-universal on H S ; in particular, we cannot directly express H = H + H c and apply H c = H − H, or else the problem would be trivial.In addition to errorcorrected Hamiltonian simulation in open quantum systems, scenarios of great practical interest may arise when the control Hamiltonians are subject to more restrictive locality constraints than the system and target Hamiltonians are (e.g., two-body simulation with only local controls, see also Sec. 4.1).
The required decoupling conditions in Assumption (iv) are automatically obeyed if the representation ρ acts irreducibly on H S .This follows directly from Schur's lemma, together with the fact that the map F Γ defined in Eq. ( 20) is trace-preserving, and both H S and S α can be taken to be traceless.While convenient, irreducibility is not, however, a requirement.When the representation ρ is reducible, care must be taken in order to ensure that Assumption (iv) is nevertheless obeyed.It should be stressed that this is possible independently of the target Hamiltonian H. Therefore, if the choice (G, ρ) works for one Eulerian simulation scheme (whether ρ is irreducible or not), then it can be used for Eulerian simulation with any target H that belongs to the reachable set from H, that is, that can satisfy Eq. (22).
We close this discussion by recalling that it is always possible for a finite-dimensional target system S to find a control group G for which both Assumptions (ii) and (iv) are satisfied, by resorting to the concept of a transformer [22,14].A transformer is a pair (G, ρ), where G is a finite group and ρ : G → U(H S ), g → ρ(g) = U g is a faithful, unitary, projective representation such that, for any traceless Hermitian operators A and B on H S with A = 0, one may express We illustrate this general idea in the simplest case of a single qubit, H = H S = C 2 .Let X, Y, Z denote the Pauli matrices and R the unitary matrix defined by which corresponds to a rotation by an angle 4π/3 about an axis n ≡ (1, 1, 1)/ √ 3. Direct calculation shows that R 3 = I and that conjugation by R cyclically shifts the Paulimatrices, i.e., R † XR = Y, R † Y R = Z, and R † ZR = X.Consider now the group G given by the presentation Using the defining relations of this group, its elements can always be written as x a z b r c , where a, b ∈ {0, 1} and c ∈ {0, 1, 2}.Clearly, the assignment ρ given by x → X, y → Y, z → Z, r → R yields a faithful, unitary, irreducible representation since the Pauli matrices commute up to phase.It is shown in [22] that the pair (G, ρ) defines a transformer in the sense given above, namely, any 2 × 2 traceless matrix B may be reached from any fixed 2 × 2 traceless, nonzero matrix A, for suitable nonnegative weights w g .The irreducibility property for any transformer pair can be easily established by contradiction [44].
A drawback of the transformer formalism is that general transformer groups tend to be large, making purely transformer-based simulation schemes inefficient.In practice, given the native system Hamiltonian H S , the challenge is to find a group G that grants a reasonably efficient scheme while satisfying Assumptions (ii) and (iv), and subject to the ability to implement the required control operations.As we shall see next, transformerinspired ideas may still prove useful in devising simulation schemes in the presence, for instance, of additional symmetry conditions.

Illustrative applications
In this section, we explicitly analyze simple yet paradigmatic Hamiltonian simulation tasks motivated by QIP applications.While a number of other interesting examples and generalizations may be envisioned (as also further discussed in the Conclusions), our goal here is to give a concrete sense of the usefulness and versatility of our Eulerian simulation approach in physically realistic control settings.In particular, we focus on achieving non-local Hamiltonian simulation using only bounded-strength local (singlequbit) control, in both closed and open multi-qubit systems.

Eulerian simulation in closed Heisenberg-coupled qubit networks
Let us start from the simplest case of a system consisting of n = 2 qubits, interacting via an isotropic Heisenberg Hamiltonian of the form where J has units of energy and the second equality defines an equivalent compact notation.We are interested in a class of target XYZ Hamiltonians of the form For instance, J x = J y = ±J, J z = 0 corresponds to an isotropic XX model, whereas if J x = J y with J z = 0, an XXZ interaction is obtained, the special value J z = ∓2J corresponding to the important case of a dipolar Hamiltonian.The construction of a simulation protocol starts from observing that Hamiltonians as in Eq. ( 34) are reachable from H, in the sense of Eq. ( 22), based on single-qubit control only.Specifically, let G ≡ Z 2 × Z 2 ≡ Z 2 2 , and let the representation ρ map (n, m) ∈ G to X n Z m ⊗ I.That is, G is mapped to the following set of unitaries: Choosing the generators of G to be (1, 0) → γ x,1 = X 1 and (0, 1) → γ z,1 = Z 1 , we assume that we have access to the control Hamiltonians where the control inputs f x (t) and f z (t) satisfy f u (0) = 0 = f u (∆) and ∆ 0 f u (τ )dτ = π/2, for u = x, z.Recalling Eq. ( 17), this yields the control propagators with u x (∆) = X 1 and u z (∆) = Z 1 (up to phase), as desired.
Note that for any single-qubit Hamiltonians A and B, averaging over the unitary group in Eq. ( 35) results in the following projection super-operator: In general, the map F Γ is trace-preserving and, in this case, it acts non-trivially only on the first qubit.Thus, F Γ is trace-preserving on the first qubit.Since each term in H is traceless in the first qubit, the decoupling condition Π G [F Γ (H)] = 0 follows directly from Eq. ( 36), even though the relevant representation ρ is, manifestly, reducible.Having satisfied our main requirements for Eulerian simulation, reachability of XYZ Hamiltonians as in Eq. ( 34) is equivalent to the existence of a solution to the following set of conditions: for non-negative weights w g .While infinitely many choices exist in general, minimizing the total weight W = g w g keeps the simulation time overhead to a minimum.For instance, it is easy to verify that a dipolar Hamiltonian of the form may be simulated with minimum time overhead by choosing weights The Cayley graph associated with the resulting Eulerian simulation protocol is depicted in Fig. 2, with the explicit timing structure of the control block as in Fig. 1 and  34), which is proportional to the time τ g = w g T spent at vertex g during the coasting subinterval; see also Fig. 1.
N = 2 × 4 = 8 control segments per block.It is worth observing that although the weights w X 1 and w Y 1 are zero in the particular case at hand, all group members of G are nonetheless required, and the unitaries X 1 and Y 1 still show up in the simulation scheme (during the ramping-up sub-intervals, as evident from Eq. ( 26)).This is crucial to guarantee that the unwanted F Γ term is projected out.The above analysis and simulation protocols can be easily generalized to a chain of n qubits (or spins), subject to nearest-neighbor (NN) homogeneous Heisenberg couplings, that is, a Hamiltonian of the form where for later reference we have introduced the standard compact notation σ i ≡ (X i , Y i , Z i ) and we assume for concreteness that n is even.In this case, we need only change the unitary representation ρ of Z 2 × Z 2 to be defined by the two generators Physically, the required generators γ x,odd and γ z,odd correspond to control Hamiltonians that are still just sums of 1-local terms, and that act non-trivially on odd qubits only: We expect that the design of Eulerian simulation schemes for more general scenarios where both the input and the target (H, H) are arbitrary two-body Hamiltonians (including, for instance, long-range couplings) will greatly benefit from the existence of combinatorial approaches for constructing efficient DD groups [45,41].A more indepth analysis of this topic is, however, beyond our current scope.

Error-corrected Eulerian simulation in open Heisenberg-coupled qubit networks
Imagine now that the Heisenberg-coupled system S considered in the previous section is coupled to an environment B, and the task is to achieve the desired XYZ Hamiltonian simulation while also removing arbitrary linear decoherence to the leading order.The total input Hamiltonian has the form (38) where H B and B u,i , for each i and u = x, y, z, are operators acting on H B , whose norm is sufficiently small to ensure convergence of the relevant Magnus series, similar to first-order DCG constructions [27,28].The target Hamiltonian then reads in terms of suitable coupling-strength parameters J u as in Eq. (34).As before, we start by analyzing the case of n = 2 qubits in full detail.Our strategy to synthesize a dynamically corrected simulation scheme involves two stages: (i) We will first decouple S from B, while leaving the system Hamiltonian H S = H iso unaffected; (ii) We will then apply the closed-system protocol of Sec.4.1 to convert H iso into the target system Hamiltonian HS = H XYZ .Once a suitable group and weights are identified in this way, both stages are carried out simultaneously in application.
A suitable DD group able to suppress general linear decoherence is provided by G DD = Z 2 × Z 2 , under the n-fold tensor power representation yielding (see also [28]): generated, for instance, by the two collective generators γ x,all = X (all) = X 1 X 2 and γ z,all = Z (all) = Z 1 Z 2 .In addition to the order of G GL being minimal, with |G GL | = 4 independently of n, step (i) above is automatically satisfied for the input Hamiltonian at hand, since Given a generic operator A on H = H S ⊗ H B , we may define the superoperator Φ DD as A + X (all) AX (all) + Y (all) AY (all) + Z (all) AZ (all) , corresponding to weights {w h } given by w In step (ii), we still rely on the group Z 2 × Z 2 , but now under a different representation.We choose the representation yielding the set G 1 of Eq. (35), with the same single-qubit generators γ x,1 = X 1 , γ z,1 = Z 1 , and the corresponding weights {w g 1 } determined by the solution of Eqs.(37).Define the superoperator Φ 1 to act as Then the combined action of the two superoperators Φ DD and Φ 1 yields where , with unitary representation elements corresponding to the full Pauli group on two qubits: The above representation is irreducible, with Π G implementing the complete depolarizing channel on two qubits: for every input A. Together with the fact that all of the system terms in H are traceless and F Γ is trace-preserving, this ensures that the DD conditions of Eqs. ( 31)-( 32 A practically important case, where simpler simulation schemes are possible, occurs if qubits couple to their environment along a fixed axis, effectively corresponding to pure dephasing -say, for concreteness, that B y,i = 0 = B z,i for i = 1, 2 in Eq. (38).A smaller DD group suffices in this case [28], namely G DD = Z 2 , represented again in terms of collective qubit rotations, and generated by the single element γ z,all .Clearly, the commutation relationship in Eq. ( 39) is maintained, still allowing our two-step procedure to be followed.In this case, the combined group for simulation is 2 , with |G| = 8, |Γ| = 3, reducibly represented as follows on the two-qubit space: Suppose, for instance, that the task is to simulate a dipolar Hamiltonian H dip as in Sec.4.1.By following the above general procedure, with weights {w h } given by w I = w Z 1 Z 2 = 1/2 for G D alone, it is easy to see that Eq. ( 40) simplifies, leading to simulation weights w I = 1/4, w Z 1 = 3/4 = w Z 2 , w Z 1 Z 2 = 1/4, with the remaining 4 weights equal to 0. While this implies that the simulation can now be achieved with only N = 8 × 3 = 24 segments per cycle and minimum weight W = 2, care is needed in ensuring that the DD conditions in Eqs. ( 31)- (32), are still obeyed.This may be checked by inspection.In particular, the fact that Π G [F Γ (X i )] = 0 for i = 1, 2 follows by analyzing the structure of each toggling-frame "error Hamiltonian", u † γ j (t)X i u γ j (t), for γ j ∈ Γ = {X 1 , X 2 , Z 1 + Z 2 }, and verifying that no term proportional to Z 2 is generated, that would be left uncorrected by averaging over the representation in Eq. (41).Likewise, the fact that Π G [F Γ (H S )] = 0 for H S = H iso may be directly established by a similar calculation, or by using the trace argument in Sec.4.1 for the two group generators γ x,1 = X 1 and γ z,1 = Z 1 , while also noting that for the third generator γ z,all = Z 1 Z 2 , we have F Z 1 Z 2 (H iso ) = H iso and the latter is decoupled by the representation in Eq. ( 41), Π G (H iso ) = 0. Thus, Eulerian Hamiltonian simulation in the presence of single-axis errors can be efficiently achieved.
Again, the schemes we have just presented for n = 2 can be generalized to a chain consisting of n spins, which interact according to a NN Heisenberg interaction and are each linearly coupled to the environment, according to Eq. (38).In this case, exploiting the results of Sec.4.1, a useful group for simulation is provided by G Z 4 2 , under the unitary representation {U g } ≡ G GL × G odd , corresponding to generators γ x,all , γ z,all , γ x,odd , γ z,odd , all of which can be implemented using only 1-local (single-qubit) Hamiltonians.As before, each simulation cycle will consist in the general case of arbitrary linear decoherence of N = 16 × 4 = 64 time segments.Despite the reducibility of the above representation (with the full Pauli group on n qubits consisting of 4 n elements), the DD conditions given by Eqs. ( 31)- (32) remain valid for reasons similar to those outlined for n = 2 under pure dephasing.

Eulerian simulation of Kitaev's honeycomb lattice Hamiltonian
We return to Eulerian simulation in closed quantum systems, but tackle a more complicated Hamiltonian of paradigmatic relevance to topological quantum memories, namely, Kitaev's honeycomb lattice model [46].Suppose that the target system consists of a network of qubits arranged on a honeycomb lattice and interacting via NN Ising couplings.The relevant Hamiltonian H is graphically displayed in Fig. 3  The basic idea to accomplish this simulation is to exploit the matrix R given in Eq. ( 33), in conjunction with the symmetry of our problem: since all Hamiltonian terms are precisely two-local and of the homogeneous form σ ⊗ σ, it will be possible to avoid using the full machinery of a transformer.Consider the group G generated by the three unitaries, ρ X , τ X , and R global , where ρ X , shown in Fig. 4(left) with σ = X, has X's on every second forward-slash, τ X , shown in Fig. 4(center) with σ = X, has X's on every second back-slash, and R global , shown in Fig. 4(right), has R applied to every vertex.These unitaries can be generated by one-local Hamiltonians.By repeatedly conjugating ρ X and τ X with R global , we immediately see that we can also perform ρ σ and τ σ , shown in Fig. 4, for any σ = X, Y, Z.Note that up to phase, all such ρ and τ commute.Because conjugation by R maps Pauli matrices to Pauli matrices, for any Pauli σ we have Rσ = (RσR −1 )R = σ R, where σ is another Pauli matrix.Thus, up to phase, we can write any element of G in the canonical form where and R a global only appears on the right.To construct an Eulerian simulation protocol we must be able to choose w g so that H is reachable from H, i.e., obeys Eq. ( 22), while ensuring that the DD condition of Eq. ( 21) is also fulfilled.We start from the fact that Observe that when U g = ρ X , all forward-slash edges connect vertices that are acted upon by either I ⊗ I or X ⊗ X, while all other edges connect vertices that are operated by X ⊗ I. Consequently, 1 2 I † HI + 1 2 ρ † X Hρ X removes all Hamiltonian terms except for those along the forward-slashes; upon conjugating by R global , we may then convert these surviving ZZ terms to XX terms, as desired.To summarize, gives the Hamiltonian shown in Fig. 5(left).Similarly, the effect of 1 2 I † HI + 1 2 τ † X Hτ X is to leave precisely the back-slash edges, which can be converted from ZZ to Y Y by conjugation by R 2 global .Thus, ) gives the Hamiltonian shown in Fig. 5(center).Lastly, it is not hard to see that the product ρ X τ X has X's on every second row of verticals; accordingly, Φ ZZ (H) ≡ 1 2 isolates precisely the verticals, giving the Hamiltonian shown in Fig. 5(right).In this case, no R-conjugation is necessary since we wish to maintain ZZ edges along the verticals.Putting all these steps together, we conclude that thus providing the desired weights for the Eulerian protocol.Since there are |Γ| = 3 generators and, from Eq. ( 42), |G| = 4 × 4 × 3 = 48 group elements, each control block consists of N = 144 time intervals.Lastly, we must verify that Eq. ( 21) holds.Note that F Γ (H) acts via conjugating each vertex by unitaries (since the generating pulses are one-local), and since such an operation is trace-preserving at each vertex, this necessarily takes the precisely twolocal terms in H to precisely two-local terms in F Γ (H).Since no one-local terms can arise, all terms are of the form σ Due to the canonical form of our group elements, Eq. ( 42), the action of Π G reads where τ ∈ {I, τ X , τ Y , τ Z } and ρ ∈ {I, ρ X , ρ Y , ρ Z }, respectively.Just as the map 1 2 IHI + 1 2 ρ X Hρ X removes all non-forward-slash ZZ terms, the map ρ ρF Γ (H)ρ depolarizes precisely one vertex of each pair of non-forward-slash vertices, and therefore suppresses all non-forward-slash terms.With only forward-slash terms remaining, τ τ [ ρ ρF Γ (H)ρ]τ = 0, since the τ -sum removes all non-back-slash terms.Thus, we conclude that Π G [F Γ (H)] = 0, as desired.

Conclusion and outlook
We have shown that the Eulerian cycle technique successfully employed in both dynamical decoupling schemes and dynamically corrected gates can be extended to also enable Hamiltonian quantum simulation with realistic bounded-strength controls.For given internal dynamics and control resources, we have characterized the family of reachable target Hamiltonians and provided constructive open-loop control protocols for stroboscopically implementing a desired evolution in the family with accuracy (at least) up to the second order in the sense of average Hamiltonian theory.We have additionally shown how Hamiltonian simulation may be accomplished in an open quantum system while simultaneously suppressing unwanted decoherence, provided that appropriate time-scale requirements and decoupling conditions are fulfilled.The usefulness and flexibility of our Eulerian simulation techniques have been explicitly illustrated through several QIP-motivated examples involving both unitary and open-system dynamics on interacting qubit networks.In all cases, access to purely local (single-qubit) control Hamiltonian is assumed, subject to finite-amplitude constraint.
It is our hope that our results may be of immediate relevance to ongoing efforts for developing and programming quantum simulators in the laboratory.A a number of possible generalizations and further theoretical questions may be worth considering.As an additional simulation problem dual to the one we analyzed for Heisenberg-coupled spin chains, exploring schemes where a target Heisenberg Hamiltonian is generated out of Ising couplings only would be of special interest, given the experimental availability of the latter in existing large-scale trapped-ion simulators [16].Likewise, it could be useful to explore whether bounded-strength simulation as proposed here may be made compatible with open-loop filtering techniques for modulating coupling strengths, such as proposed in [47], as well as in [48] in conjunction with non-unitary control via field gradients.Building on existing results for dynamical decoupling schemes [42], the use and possible advantages of randomized simulation schemes in terms of efficiency and/or robustness may be yet another venue of investigation, especially in connection with large control groups.Lastly, it remains an important open question to determine whether simulation schemes able to guarantee a minimum fidelity over long evolution times may be devised, in the spirit of [49] for the particular case of the zero Hamiltonian.

Figure 1 .
Figure 1.Schematics of an Eulerian simulation protocol.The basic control block consists of N time intervals, each involving a "ramping-up" sub-interval of fixed duration ∆, during which H c (t) = 0, followed by a "coasting" (free evolution) period of variable duration Θ k , Eq. (24), during which no control is applied.During the jth ramping-up sub-interval we apply h γj , i.e., the control Hamiltonian that realizes the generator γ j , smoothly changing the control propagator from U gj−1 to U gj .In this way, the control protocol corresponding to Eqs. (26)-(27) is implemented.By construction, a standard EDD cycle with H = 0 is recovered by letting Θ k → 0 for all k, while in the limit ∆ → 0 standard BB simulation of H is implemented.

Figure 2 .
Figure 2. Cayley graph for the Eulerian simulation of the dipolar Hamiltonian in Heisenberg-coupled qubits.Vertices are labeled by group elements; edges are labeled by group generators.Numbers in parentheses next to vertices indicate the weights w g of the corresponding group elements g in Eq. (34), which is proportional to the time τ g = w g T spent at vertex g during the coasting subinterval; see also Fig.1.
) are satisfied.Since |G| = 16 and |Γ| = 4, the resulting Eulerian simulation cycle will involve in general N = 64 time segments, with the number of non-zero weights (hence the total weight W and the time-overhead of the simulation) being determined by the details of the error model and/or the target Hamiltonian.

Figure 3 .
Figure 3. Input and target Hamiltonians on a 2D honeycomb lattice, where qubits are placed at each vertex.Left: The system Hamiltonian H describes a system where all adjacent vertices have ZZ Ising couplings.Right: The target Hamiltonian H realizes Kitaev's honeycomb lattice model, with XX, Y Y , and ZZ couplings depending on the type of the edge.
(left), where vertices represent qubits and edges represent two-qubit couplings of the form Z k Z , with vertices k and being adjacent in the graph and Z k indicating, as before, the Pauli Z operator acting non-trivially only on qubit k.The target Hamiltonian H is shown in Fig.3(right), where some of the edges are now of the form X k X and Y k Y .In accordance with the figure, we shall also call the XX-edges forward-slashes, the Y Y -edges backslashes, and the ZZ-edges verticals henceforth.

Figure 4 .
Figure 4. Pictorial representation of different control operations.Left: The unitary ρ σ , with σ on the vertices of every second forward-slash and I on all other vertices, where σ is a fixed X, Y, or Z operator.When σ = X, this is the generator ρ X .Center:The unitary τ σ , with σ on the vertices of every second back-slash, where σ is a fixed X, Y, or Z operator.When σ = X this is the generator τ X .Right: The generator R global , with R at every vertex.

Figure 5 .
Figure 5. Pictorial representation of different simulation superoperators (see text).Left: Action of the superoperator Φ XX , leaving XX terms at forward-slashes only.Center: Action of the superoperator Φ Y Y , leaving Y Y terms at back-slashes only.Right: Action of the superoperator Φ ZZ , leaving ZZ terms at verticals only.
k and are adjacent vertices and σ u , σ v ∈ {X, Y, Z}.Thus, we may write