Parameterized Hamiltonian simulation using quantum optimal control

Analog quantum simulation offers a hardware-specific approach to studying quantum dynamics, but mapping a model Hamiltonian onto the available device parameters requires matching the hardware dynamics. We introduce a paradigm for quantum Hamiltonian simulation that leverages digital decomposition techniques and optimal control to perform analog simulation. We validate this approach by constructing the optimal analog controls for a superconducting transmon device to emulate the dynamics of an extended Bose-Hubbard model. We demonstrate the role of control time, digital error, and pulse complexity, and we explore the accuracy and robustness of these controls. We conclude by discussing the opportunity for implementing this paradigm in near-term quantum devices.


I. INTRODUCTION
Estimating dynamical properties of entangled manybody quantum states is critical for scientific discovery and engineering robust quantum technologies [1,2]. Unfortunately, even the most sophisticated algorithms are limited to simulating quantum dynamics at small times due to the exponential growth in memory requirements with increasing quantum state entanglement [3][4][5]. Quantum computers, on the other hand, are capable of efficient simulation of quantum dynamics because even highly entangled quantum states can be represented efficiently in a quantum register [6].
Simulating the dynamics of a model quantum system using a controllable quantum device is referred to as quantum simulation [1]. Two leading paradigms for quantum simulation are based on digital and analog representations. In digital quantum simulation, sequences of discrete quantum logic gates are compiled to generate unitary operators that implement time evolution of the model system, while in analog quantum simulation the quantum device is controlled so that its own evolution imitates the dynamics of the model system [1].
Several approaches have combined digital quantum logic with analog quantum evolution, especially for quantum simulation [7][8][9][10][11][12][13][14][15]. These approaches augment digital simulation methods with specialized analog evolution to extend the types of quantum states that can be prepared. However, a typical hurdle for analog quantum simulation is the design or discovery of a controllable * pkairys@vols.utk.edu † This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
quantum device whose Hamiltonian is isomorphic to the model Hamiltonian [1]. In many cases, construction of the device itself limits this isomorphism to small regime of model parameters and, therefore, the class of possible simulations.
Here we address this limitation by introducing a framework to decompose analog quantum simulation into a sequence of individually controlled unitaries. By designing sequences of unitaries with optimal control, this framework demonstrates as a direct route to generating model dynamics of interest. This generalizes recent work by Holland et al., which used optimal control for analog simulation of nuclear physics [16]. In addition, we show that the resulting decomposition partitions the problem of identifying optimal analog controls into a linear set of independent optimization problems. This enhances the computational efficiency of designing the simulations and expands the class simulatable models. In justifying the feasibility of the optimal control in a simulation paradigm, we introduce a strategy to realize quantum simulation over a wide range of model parameters without additional control optimizations.
We structure the remainder of the work as follows: in Sec. II we describe our proposed simulation paradigm. In Sec. III we demonstrate the feasibility of the paradigm in three steps: Sec. III A defines the model system of interest, Sec. III B defines the controllable device, and Sec. III C discusses the numerical results of our demonstration. Finally we conclude and discuss our contributions in Sec. IV.

II. PARAMETERIZED SIMULATION WITH OPTIMAL CONTROL
Quantum simulations require specification of a model quantum system and a controllable quantum system. Typically this is accomplished by defining a parameterized model Hamiltonian H M ( λ) and a parameterized, time-dependent device Hamiltonian H D ( α, t). The parameters of the device Hamiltonian α represent controllable parameters such as tunable electric or magnetic arXiv:2105.02153v1 [quant-ph] 5 May 2021 fields. The parameters of the model Hamiltonian λ are specified over some finite range λ l,min ≤ λ l ≤ λ l,max on which interesting dynamical physics are sought.
The task of quantum simulation is that the dynamics of the quantum device over a control time T c , and the dynamics of the model system over a time T s , are equivalent. In practice, this equivalence is replaced by requiring the dynamics to be sufficiently close according to some distance measure. We consider a measure between two unitary operators given by the infidelity which is invariant to a global phase and takes a positive value between zero and one [17]. We assume the Hamiltonians are defined on a composite Hilbert space of N d-dimensional quantum systems, H = N i H i , where each Hamiltonian is a sum of k-local operators. A k-local operator acts as the identity on all but k of the d-dimensional local Hilbert spaces, H i . We will assume that the model Hamiltonian has the form where H M,l is k-local. Then, we consider a device Hamiltonian composed as a linear combination of operators with the same locality as the model Hamiltonian, i.e., both H M,l and H D,l act on the same k-local subspace l. Furthermore, we assume that the subsystems acted on by each local term are addressable via a chosen set of device parameters, i.e., ∀ l ∃ α s.t. H D,l ( α, t) = 0 ∀ l = l. This structure of the device Hamiltonian ensures that any subspace on which model evolution occurs can be addressed independently by a choice of device parameters. One approach to realizing quantum simulation is to numerically determine the parameters α for a control time T c that generates the target unitary dynamics for a specific parameter set λ and model simulation time T s as There are two consequence of this approach. First, the model and device unitaries become exponentially large with increasing quantum system size. Second, for every desired set of parameters λ, T s , a new set of control parameters and control time must be determined. This strategy is more than exponentially inefficient. We overcome these bottlenecks by assuming the model Hamiltonian parameters are bounded and discretely parameterized. Explicitly, we assume that each parameter is bounded as λ l,min ≤ λ l ≤ λ l,max ∀ l and discretized as a grid with step size ∆λ l > 0 for each parameter l such that λ l = λ l,min + n l ∆λ l . Next, we use the Trotter decomposition to convert the global model unitary operator into a product of k-local unitary operators [1,18] = lim One typically truncates the series at an order q >> 1 that gives a desirable global error defined by the conservative error bound Based on Eq. (8), a target unitary generated by a Hamiltonian of L parameters requires determining only 2L unitaries. Furthermore, to implement these unitary operators in a quantum device, one solves 2L optimization problems corresponding to λ l,min as min α g(U D,l ( α, T c ), U M,l (λ l,min , T s /q)), (10) and ∆λ l as where U D,l is the unitary operator generated by addressing the subspace on which H M,l acts as The above result demonstrates that using optimal control within the context of quantum simulation is efficient, in principle. In practice, additional structure to the model Hamiltonian permits more efficient implementations. For example, if the model Hamiltonian is sufficiently uniform, H M,l ≈ H M,l ∀ l, l , and the device Hamiltonian is sufficiently uniform, H D,l ≈ H D,l ∀ l.l , then one can develop control ansatze that are capable of determining the controls across any subset of the device, requiring only minor parameter refinements (or calibrations) for each subsystem.
The former condition implies that the model Hamiltonian has internal symmetry as is common in condensed matter or chemical physics [1]. The latter condition, that the quantum device is relatively uniform, will typically be satisfied from an engineering perspective as it decreases device complexity and cost [19,20]. We discuss extensions and other applications of this approach in Appendix A.

III. DEMONSTRATION
We next demonstrate the validity and feasibility of optimal control based quantum simulations through a demonstration: emulating the dynamics of the extended Bose-Hubbard (EBH) model in a device of interacting transmon superconducting circuit elements.

A. Application to Bose-Hubbard dynamics
The Bose-Hubbard (BH) model describes spinless bosons moving in a discrete space, like a lattice or graph. The BH model is particularly relevant to studying ensembles of cold atoms in optical lattices, light-matter interactions, multi-particle quantum walks, and dynamical quantum phase transitions [4,[21][22][23]]. The EBH model includes an additional term that has been shown to be a building block for novel phases of matter including topologically entangled quantum spin liquid phases, making it of immediate interest for near-term quantum simulation [24].
Recently, analog quantum simulations of the BH model have been performed in systems of capacitively-coupled, tunable-frequency transmons [23,25]. This is possible because capacitively coupled transmons in the dispersive regime have a Hamiltonian isomorphic to the BH Hamiltonian. However, parameter ranges for the simulation are limited due to hardware constraints and the device hardware used in these simulations do not support analog implementations of the EBH model, which is the sum of the typical BH model Hamiltonian with an additional extension term [26].
For a one-dimensional chain of lattice sites, the instance of the EBH Hamiltonian considered in this work is where V quantifies the strength of the two-site potential, J quantifies the kinetic energy of particle hopping, and the operatorsb † i ,b i ,n i are the bosonic creation, annihilation operators, and number operators acting on site i, respectively. In this representation, H EBH , is the sum of two-local operators and thus is a variant of Eq. (5).
In our demonstration, we explore dynamical properties in the subspace spanned by states of {|0 i , |1 i } ∀ i for various parameter values V /J = −1, −2, −3. This is accomplished by fixing J = −0.1/2π GHz, V min = ∆V = 0.1/2π GHz. Then, the Hamiltonian Eq. (13) is parameterized by a single, discrete parameter n v = 1, 2, 3 that corresponds the weight V /J, according to the decomposition presented in Sec. II.
Following with the approach outlined in Sec. II, we transform the global time-evolution operator into a product of two-site unitary operators using a Trotter decomposition as where the parameter n v determines the magnitude of V /J, as in Eq. (7). We note that for each pair of sites, i, i + 1, there is the BH evolution unitary and the EBH potential evolution unitary As mentioned in Sec. II these operators are locally the same for any two sites i, j and i , j . Provided that the device Hamiltonian is sufficiently uniform, we will assume that a suitable control ansatz used to find optimal controls for a single pair of sites will also be valid for any other pair of sites across the whole device. Essentially reducing the optimal control task to the solution of two optimization problems. One is for BH evolution, and the other one is for EBH potential evolution, where the superscripts labeling sites i, i + 1 have been removed to emphasize these local unitaries are invariant with respect to site labeling. The lone hyper-parameter introduced in Sec. II is the order of the Trotter decomposition, q. An accurate simulation requires that the ratio Ts q << 1. However, varying q also changes the unitaries to determine via optimal control, and we evaluate the role of q on the fidelity of the optimal controls by fixing T s = 1 ns and varying the Trotter order 1 < q < 10 in Sec. III C.

B. Device Hamiltonian
Quantum devices composed of superconducting circuits are a leading paradigm for quantum computing and systems of tunable-coupling, tunable-frequency transmons are a particularly mature instance of that technology at large scales [19,[27][28][29][30]. Tunable couplers can be used to effectively turn-off interactions between transmons when those interactions are undesired. This allows addressing subsystems of the quantum device with high fidelity; one of the primary assumptions in deriving our main result of Sec. II [28].
It is common to model such devices as systems of interacting Duffing oscillators on a planar lattice with a tunable coupling [23]. The local Hilbert space of each transmon is a Fock space, isomorphic to the local Hilbert space of a lattice site in the Bose-Hubbard model. Thus we map the model Hilbert space to the device Hilbert space by mapping each lattice site i in the model Hamiltonian Eq. (13) to a single transmon i. This allows us to write down a device Hamiltonian with the same geometry as the model Hamiltonian as, i.e. a one dimensional chain of N transmons where ω i , δ i are the frequency and anharmonicity of transmon i, respectively, and γ i,i+1 is the controllable (time-dependent) coupling of transmon i to transmon i + 1. The tunable coupling addresses pairs of transmons individually from the rest of the device. When a pair of transmons interact, the strength of the interaction is related to the the strength of the coupling γ i,i+1 and the relative detuning of the transmon frequencies ω i − ω i+1 . The strongest interaction is present when the transmons are tuned to the same frequency, ω int. , which is the case we will consider here [31]. Thus, the Hamiltonian governing dynamics of an interacting two-transmon subsystem is Moreover, quantum information processing in transmon systems is performed in reference to a pre-calibrated rotating frame associated with the idling frequency of each transmon. This transformation is given by the unitary transform which yields an effective Hamiltonian in this rotating frame, , and the corresponding time-ordered evolution operator for control time T c of a two-transmon subsystem is  1. Amplitude constrained infidelities. a, heat-map showing the minimum observed infidelities after optimization for control times (y-axis) and Trotter decomposition order q (x-axis). For each parameter pair (pixel coordinate) 10 optimizations were run with independent initial parameter choices for a control ansatz of 20 Gaussian functions. A q-dependent infidelity is observed demonstrating approximate control times of 75 ns required to obtain arbitrarily low infidelity. b, heat-map showing the average observed infidelities of all 10 optimization instances for each parameter pair. The slow, monotonic decrease in infidelity below Tc = 75 ns is due to increasing q decreasing the norm of Trotter-step unitary and therefore the norm of the infidelity.
FIG. 2. Parameterization constrained infidelities. a, heat-map showing the minimum observed infidelities after optimization for control times (y-axis) and Trotter decomposition order q (x-axis). For each parameter pair (pixel coordinate) 10 optimizations were run with independent initial parameter choices for a control time of 100 ns. We observe that the infidelity has a strong dependence on pulse complexity, as determined by the number of Gaussian basis functions. b, heat-map showing the average observed infidelities of all 10 optimization instances for each parameter pair. The slow, monotonic decrease in the mean infidelity is due to increasing q decreasing the norm of Trotter-step unitary and therefore the norm of the infidelity. the device Hamiltonian for a two-transmon system, we specify these optimization problems as and min γ(t) where we drop site labels for clarity. We remark that the first optimization problem Eq. (29) is trivial to solve in the computational subspace because the model and device Hamiltonians in this subspace are isomorphic. We use a single Gaussian function describing γ(t) to achieve infidelities below 10 −6 .
The non-trivial optimization problem given by Eq. (30) requires numerical exploration. We note that Eq. (30) is currently framed as a functional minimization problem, of γ(t). However, to minimize this function numerically requires expressing γ(t) with a finite number of parameters. We choose to do this by decomposing the coupling as the sum of multiple Gaussian functions where we define the parameter vector α = [a 0 , µ 0 , σ 0 , . . . , a M , µ M , σ M ].
This yields a numeri-  Table I. We observe a saturation of the upper bound and a sharp dip approaching the lower bound as a general topological feature of the pulse family. The control pulses are smooth and start and stops at zero with slow ramp times. b, periodograms for the family of pulses shown in a. For the optimum pulse (darkest blue) spectral density is concentrated within the range of [0, 1] GHz. This is expected due to the sharp dip in negative amplitude of the control pulse centered around 55 ns and the resonant transmon frequencies. c, The convergence of the GOAT optimization routine for various initial guesses. The most optimal pulses converge within 200 iterations however some take longer to converge. The maximum permitted iterations were 500, which was saturated by two optimizations out of the ten run. Detailed information about the optimization procedure can be found in Appendix B 2.
cal optimization task The total number of Gaussian functions M is varied to understand the impact of control complexity on the optimal infidelities. The optimization methods, including discussion of leakage levels, algorithms, and implementations are detailed in Appendix B.

C. Numerical Results
A quantum simulation enabled by optimal control requires understanding the relationship between the control time T c and the Trotter decomposition order q. We evaluate this parameter space over a grid of control times 25 ns ≤ T c ≤ 200 ns and Trotter order 1 ≤ q ≤ 10 by performing an ensemble of optimizations at each (T c , q) instance. Figure 1 plots for each pair of parameters (a) the minimum observed infidelity and (b) the mean observed infidelity obtained from each ensemble. In Fig. 1a, two distinct regimes appear that correspond to infidelities of approximately 10 −3 and 10 −12 . We refer to these regions as as high-infidelity and low-infidelity regimes, respectively.
In Fig. 1a, we observe that control times below 75 ns are unable to achieve infidelities of 10 −12 regardless of the Trotter-decomposition order. This may be related to a quantum speed limit for the controls, which represents a lower bound on the time required to implement a unitary operator with zero infidelity [32][33][34]. In fact, a crucial feature of quantum speed limits is observed here: the sharp decrease of the infidelity between regimes as the control time is increased [34]. We also observe that for T c > ∼ 75 ns, the boundary between high and low-infidelity regimes is a function of the Trotter order q, such that larger Trotter orders yield lower infidelities for shorter control times. Hence, the optimal analog quantum simulation of the EBH model on this transmon device is both fast and accurate.
In Fig. 1b, we present the mean of the infidelity. Regions with small minimum infidelities typically have small mean infidelities. We also observe a general decrease in infidelity with increasing q. This is due to the Trotter-step unitary, Eq. (17), converging to the identity operator with large q. The identity operator is a trivial solution for the optimal control problem, thus the controls which generate the identity are a type of attractor for the optimization routine, leading to lower fidelity solutions with increasing q.
Our simulations have used control pulses as linear combinations of Gaussian functions since Gaussian pulses permit control over the pulse bandwidth and are commonly used in experimental superconducting systems [35,36]. However, pulse complexity impacts how accu- is generated using a 50 ns Gaussian pulse, with parameters a = −2.07 MHz, µ = 25 ns, σ = 6.36 ns, which we found provided infidelity below 10 −6 . The EBH Trotter-step unitary and the BH Trotter-step unitary are then repeatedly applied to the initial state up to N = 60 times because the Trotterized time step is 1/6 (ns). In each column the number of applications of UE(1/6) is varied so that the equivalent ratio of V /J = 1, 2, 3 for columns 1,2, and 3, respectively. This generates the discrete evolution plotted as orange circles and is compared with the exact evolution shown in blue lines. This demonstrates that by selectively compiling the bare pulses that generate the individual Trotter steps, it is possible to reconstruct the dynamics of an EBH Hamiltonian with arbitrary V /J. rately the target unitaries are prepared [37,38]. We evaluate the dependence of infidelity on pulse complexity in Fig. 2 by fixing T c = 100 ns and varying the number of Gaussian functions from 5 to 35. We did not observe significant changes in the minimum and mean infidelities when using a control pulse with more than 35 Gaussians.
In Fig. 2a, we show a high-infidelity regime and a lowinfidelity regime with a sharp transition from ≈ 10 −4 infidelity to ≈ 10 −12 infidelity, respectively. The boundary between the regime is a function of control complexity and Trotter order. But, the dependence on these parameters is qualitatively different from Fig. 1 in that given sufficient pulse complexity, the transition between high and low-infidelity regions is only dependent on q, and dependence on the number of Gaussians vanishes. This is to be expected, as the optimal configuration of N Gaussians exist as an optimal configuration to a parameterization of N + 1 Gaussians.
In Fig. 2b, we observe that increasing the number of Gaussians leads to a decrease in the mean observed infidelity, regardless of Trotter order q. However, we observe a discrepancy in the trend in infidelity with increasing q for N = 10 and N = 35 Gaussians. These pulse complexities typically have higher mean infidelities with in-creasing q compared to other pulse complexities N , which demonstrates that the optimization landscape is a function of pulse complexity [33].
We present a family of optimal pulses in Fig. 3 that are obtained by running 10 optimizations at control time T c = 100 ns and Trotter decomposition order q = 6 with a control ansatz of 20 Gaussian functions. In Fig. 3a, we present the time-domain representation of the optimal pulses. We note that the family of optimal pulses share a number of common features: ramps to positive coupling, sharp dips into negative coupling, then back to positive coupling and ramps to zero coupling. In particular, the family of lowest infidelity pulses vary slightly in the location of the negative amplitude peak, indicating some insensitivity of the infidelity on peak location within the range 40 − 60 ns.
Periodograms of the optimal pulse family are shown in Fig. 3b. Nearly all of the spectral intensity is concentrated in frequency components ≤ 1 GHz. Finally, in Fig. 3c, we show the convergence of the pulse optimization for each pulse in the family. We observe that the optimal pulses with the best fidelity converged most quickly, typically within 200 iterations.
A final demonstration of the ability to compile opti-mized analog pulses into sequences permitting simulation with variable Hamiltonian parameters is shown in Eq. (8). For the EBH model from Eq. (15), the unitary evolution governed by an EBH Hamiltonian for different ratios of |V /J| is generated by repeated application of multiple Trotter steps of Eq. (17). We use the lowest infidelity pulse from Fig. 3a that generates evolution of U E (T s /q) and compile it with the single Gaussian pulse that generates U BH (T s /q) for variable |V /J|. In Fig. 4, we show the real and imaginary components of the probability amplitude of the state |11 as a function of the number of applied unitaries, U E (T s /q). We observe that by compiling pulses found via optimal control, dynamics can be generated for Hamiltonians with variable ratios of |V /J|. This evolution is visualized as a quantum circuit for clarity, but represents a smooth analytic pulse formed by concatenating the individual pulses that generate U E (T s /q) and U BH (T s /q).

IV. DISCUSSION AND CONCLUSION
Analog and digital quantum simulation are both promising approaches to studying Hamiltonian dynamics of quantum states. While several recent approaches to combine these ideas are driven by digital quantum computation [8,9,16], we have introduced a paradigm that uses optimal quantum control to implement Hamiltonian simulation protocols. First, we have shown how a parameterized Trotter decomposition of a target unitary propagator limits the number of local optimal control problems to a linear function of the device Hamiltonian parameters. We then compiled the resulting optimal pulses to simulate the target Hamiltonian for arbitrary model parameter values.
We have demonstrated Hamiltonian simulation by considering a dynamics simulations of the extended Bose-Hubbard model within a system of superconducting transmons. Notably, the model Hamiltonian is not isomorphic with the target Hamiltonian and the model parameter values lie outside of the device parameter regimes. Using an ensemble of numerical simulations we have explored the role control time, Trotter decomposition order, and control complexity on the obtainable infidelities using these optimal control pulses.
Our numerical results show that fast, high-fidelity control pulses are present over a range of Trotter orders and pulse complexities. The discovered optimal controls are smooth and concentrated in the low frequency domain, which supports the future realization in near-term quantum hardware. Together, these results emphasize that compilation of the optimal controls enables simulation of Hamiltonians with arbitrary parameter values.
In Section II, we decompose the global unitary time evolution of a model Hamiltonian into a sequence of local unitary evaluations using the Trotter decomposition, and demonstrate how one can further decompose the evolution into a small number of local unitaries when the Hamiltonian is appropriately parameterized. This decomposition is not unique and depending on the application, alternative decompositions may prove to be more practical. Here we consider another example of a parameterized model Hamiltonian in which each the parameters of λ lie in a finite range around zero: λ l,min ≤ 0 ≤ λ l,max and the parameter range is discretized on a fixed grid with spacing ∆λ l > 0. Then we can define an arbitrary point in the parameter space by two parameters 0 ≤ n l,− ≤ N l,− and 0 ≤ n l,+ ≤ N l,+ with an expression λ l = (n l,+ − n l,− )∆λ l .
Using the above parameterization, we consider the Trotter decomposition of the global unitary operator into a product of k-local unitary operators to produce a modified version of Eq. (7) and Eq. Comparing this decomposition to Eq. (7) and Eq. (8) we see that depending on the parameterization of the Hamiltonian, a decomposition of the time evolution operator can be chosen to limit the total number of unitary applications per Trotter step. For example, if the goal is to simulate a Hamiltonian with λ l = 2∆λ l then the above decomposition requires only two unitaries to be applied per Trotter step, whereas if one were to use Eq. (8) this same simulation would require N l,− + 2 applied unitaries per Trotter step.
In addition to selecting a preferred decomposition based on the Hamiltonian's parameterization we also re-mark on the ability to add local disorder to the quantum dynamics simulation. In this case we can consider that a particular instance of the model Hamiltonian H M ( λ l ) is perturbed by some vector l . For simplicity, we will assume that each local perturbation l is chosen from a range of perturbation strengths l,max ≤ 0 ≤ l,min that are discretized on a grid with spacing ∆ l > 0. Then any perturbation can be described as l = (n l, ,+ −n l, ,− )∆ l , using the same notation as introduced above.
Using the above decomposition we consider the Trotter decomposition of the global unitary operator into a product of k-local unitary operators to produce a modified version of Eq. (7) and Eq. (8): which still requires finding optimal controls for 2L unitary operators but permits quantum simulation of Hamiltonians with local disorder Finally, we have focused on the case of timeindependent problem Hamiltonians but in principle they may be time dependent as well. We speculate that in the case of a time-dependent Hamiltonian one could combine our strategy with digital methods for time-dependent Hamiltonians such as the one proposed in Ref. [39], but we do not consider that generalization in this work.
Appendix B: Methods

Problem formulation
As discussed in Section III in order to implement a quantum simulation requires determining a set of control parameters α such that the unitary operator generated by the device evolution, Eq. (28), is sufficiently close to a chosen Trotter-step unitary, according to some distance measure. In this work we consider a measure between two unitary operators given by the infidelity: where α is the vector of parameters to optimize over, P c is a projector onto a desired computational subspace, U M is the desired unitary operator to implement in the computational subspace, and U D ( α, T c ) is the time evolution operator over the full system Hilbert space at time T c given by Eq. (28).
In this work we demonstrate quantum simulation by considering a system of two transmons with a tunable interaction. The local Hilbert space of each transmon H i is an infinite-dimensional Fock space, and the composite space of two transmons i = 1, j = 2, is given by H 1 ⊗ H 2 . However, to perform a computation of these systems requires first choosing a finite subspace of this composite Hilbert space.
In the dispersive regime, the effective Hamiltonian between two transmons interacting via a tunable coupling transmon, is given by Eq. (5). We note that this Hamiltonian is block-diagonal in the total excitation number basis, due to the fact that the interaction term conserves the total number of excitations. Thus, if we are interested in preparing a unitary operator in the qubit subspace spanned by the orthonormal basis {|00 , |01 , |10 , |11 }, we need only consider the Hilbert space composed of the first three levels of each transmon because this includes the subspace with total excitations N = 0, 1, 2 and the device Hamiltonian does not permit population transfer between blocks of different excitation number.
By taking into account the first three levels of each transmon in our simulations any state overlapping |11 will populate the states |02 , |20 during the control time. This is due to the fact that these quantum states share a total number of excitations N = 2. This phenomena is referred to as "leakage," and is a common issue in superconducting systems. We note that the formulation of the infidelity function, Eq. (3), will penalize unitary operators that, couple |11 to |02 or |20 at the final time T c . However, it will not penalize unitaries that couple these states at intermediate control times T < T c , which is not generally possible given the form of Eq. (20).
Based on these facts we remark that the optimal controls found in this work minimize leakage outside of the computational subspace only at the final time, not intermediate times. At first glance this seems like an unwanted effect. However, due to the device Hamiltonian, dynamics in the N = 2 excitation subspace are necessary to generate the correct final-time unitary. I.e, the additional basis states |02 , |20 are actually being used for information processing in a beneficial way, and are necessary to implement the desired evolution in the computational subspace.

Optimal Control with GOAT
There are numerous quantum optimal control algorithms for optimizing a problem like the one above, here we focus on the recently developed gradient optimization of analytic controls (GOAT) algorithm [40]. This method was chosen because of it's rapid convergence to optimal solutions and the ability to determine low parametercount controls which may be easier to calibrate on real quantum devices [34,40].
The GOAT algorithm considers Hamiltonians that are a linear combination of drift Hamiltonians that are constant and control channels that are parameterized and time-dependent: and expands each driving field in a superposition of functions which need not necessarily be orthogonal or complete. For example each function might be pulled from a Fourier basis or could be Gaussian shaped pulses, which may lead to different number of parameters per basis function.
To perform a gradient based search over the set of parameters we take the derivative of the objective function w.r.t each parameter where we identify that we must find ∂ α U ( α, T c ). But there doesn't exist an analytic expression for this operator in general. Thus it was proposed to find this quantity via an alternative equation of motion derived from the Schrödinger equation. The procedure is to take the Schrödinger equation for the unitary operator and take its partial derivative w.r.t a parameter α and exchange the order of the second derivative on the left side leading to an equation of motion for the which leads to a set of coupled differential equations: .
(B10) Thus the GOAT algorithm has two primary steps: (1) solution of the coupled Schrödinger equations for every instance of the parameters, and (2) the optimization of the parameters to minimize the objective function g. By iterating these steps until a convergence threshold is reached, we are able identify a set of control parameters that are a minimum of the infidelity, Eq. (3). The main advantage to using the GOAT algorithm over other optimal control methods are that the optimized fields are analytic and thus easily interpreted, and potentially easier to calibrate because the number of parameters is independent of the control time, unlike other methods [40].

Implementation
We implement the GOAT algorithm using the programming language Julia and various open-source pack-ages. Our implementation uses the Julia package Dif-ferentialEquations.jl to numerically solve the coupled GOAT equations of motion, Eq. (B10) using a order 5/4 Runge-Kutta method with adaptive time stepping [41]. For the gradient-based optimization we use a limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with a backtracking line-search method, both implemented in the Optim.jl package [42]. We limit each optimization routine to a maximum of 500 iterations and define a stopping criteria when the infinity-norm of the gradient falls below 1e-5.
For each of the parameters in Figures 1 and 2 we run 10 optimization problems each with various initial guesses. The initial guesses for the gradient-based optimization are chosen as the guesses with minimal infidelity observed from an ensemble of 10000 guesses based on a uniform random sampling of each parameter guess range. This large number of initial guesses was chosen to mitigate traps in local minima and because evaluating the infidelity of a control is less computationally expensive than performing optimization using GOAT. For each Gaussian basis function we select an initial guess from the parameter ranges: These parameter regimes were selected to bias the optimization a priori towards optimal solutions that are more experimentally feasible, i.e. low amplitude, low frequency controls that start and end smoothly at zero amplitude.
In all of the numerical optimizations performed we constrain the control amplitude to within a range specified by experimental restrictions, shown in Table I. To implement this constraint we construct a saturation function based on a generalized logistic function defined by: where A = min(g 12 ), B = max(g 12 ), g = 4 are the lower asymptote, upper asymptote, and maximum gradient of the sigmoid, respectively. The additional, variable A is chosen such that S(0) ≈ 0. Normally, the logistic function is non-constant within a finite domain around [A, B] thus the control field is linearly re-scaled within the specified amplitude range [A, B] by the ratio of (B − A)/2 which gives the magnitude of the input amplitude relative to the desired range.
This amplitude saturation function acts to bound the field amplitude in a continuous manner and leads to a vanishing gradient component w.r.t the control amplitude when the control amplitude lies outside [A, B], thus implementing a barrier for gradient based optimization.
Then, the controls interact with the Hamiltonian as and (for the purposes of the GOAT algorithm) the partial derivative of H(t) w.r.t. a parameter α is given by