Engineering Fast High-Fidelity Quantum Operations With Constrained Interactions

Understanding how to tailor quantum dynamics to achieve a desired evolution is a crucial problem in almost all quantum technologies. We present a very general method for designing high-efficiency control sequences that are always fully compatible with experimental constraints on available interactions and their tunability. Our approach reduces in the end to finding control fields by solving a set of time-independent linear equations. We illustrate our method by applying it to a number of physically-relevant problems: the strong-driving limit of a two-level system, fast squeezing in a parametrically driven cavity, the leakage problem in transmon qubit gates, and the acceleration of SNAP gates in a qubit-cavity system.

Understanding how to tailor quantum dynamics to achieve a desired evolution is a crucial problem in almost all quantum technologies. We present a very general method for designing high-efficiency control sequences that are always fully compatible with experimental constraints on available interactions and their tunability. Our approach reduces in the end to finding control fields by solving a set of time-independent linear equations. We illustrate our method by applying it to a number of physically-relevant problems: the strong-driving limit of a two-level system, fast squeezing in a parametrically driven cavity, the leakage problem in transmon qubit gates, and the acceleration of SNAP gates in a qubit-cavity system.

I. INTRODUCTION
The success of any nascent quantum technology will ultimately be limited by our ability to manipulate relevant quantum states. Finding the required time-dependent control fields that generate with high accuracy a desired unitary evolution is in general not a trivial task: it is sufficient to consider a simple driven two-level system in the strong driving limit [1] to find an example of a complex control problem. This generic problem becomes even more complicated when including realistic constraints: unavailable control fields, bandwidth and amplitude limitations, etc. Finding new widely applicable methods to attack such problems is thus highly desirable.
There are of course many existing approaches to quantum control. Of these, the most ubiquitous is to exploit numerical algorithms (see e.g. [2][3][4][5][6][7]) based on optimal quantum control theory [8]. The methods ultimately rely on the numerical optimization of an objective function, for example the fidelity of a desired target state with the actual time-evolved state. For many problems the effective landscape of the objective function has many local minima, which can make it challenging to find the truly optimal protocol. While methods to overcome these limitations exist [9][10][11][12], they become difficult to implement as the dimension of the control space increases. An alternative approach is to use an analytical method to design effective protocols; control pulses designed in this way could then be further improved by using them to seed a numerical optimal control algorithm. Analytic methods are however often system specific (see e.g. [13,14]), or only work with a specific restricted class of dynamics (for example methods based on shortcuts to adiabaticity, which are specific to protocols based on adiabatic evolution [15][16][17][18][19][20][21]). These approaches are also generally impractical in systems with many degrees of freedom or sufficiently complex interactions.
In this work we present a general framework for constructing control fields that realize a desired evolution, in a manner that is explicitly consistent with experimental constraints. At its heart, it allows one to use the analytic solution of a simple control problem to then find a final pulse DFT basis functions initial pulse high-fidelity pulse sequence for a more complex problem where a closed-form analytic solution is not possible. Our method has many potential virtues: it is applicable to an extremely wide class of systems and protocols, produces smooth control fields, and only requires one to numerically solve a finite set of linear equations. It builds on the recently proposed Magnus-based control method introduced in Ref. [22], but greatly extends its power and applicability.
Our generic goal is to use a specific time-dependent HamiltonianĤ(t) (whose form and tunability is constrained) to produce (at time t f ) a desired unitary operation. We start by splitting the Hamiltonian into two parts asĤ(t) =Ĥ 0 (t) +V (t), where H 0 (t) is simple enough to be analytically tractable, andV (t) represents all the additional interactions that make the problem unsolvable. The basic strategy then has two parts: (1) First, choose control fields in the "simple" Hamil-arXiv:2003.12096v1 [quant-ph] 26 Mar 2020 tonianĤ 0 (t) so that in the absence ofV (t), one realizes the desired operation. This can be done analytically.
(2) Adding backV (t) will then destroy the ideal evolution. We address this by modifying available control fields so as to average out the impact of V (t). This amounts to adding a control correction to the full Hamiltonian:Ĥ(t) →Ĥ(t) +Ŵ (t) [see Fig. 1 (a)].
The question is of course how to find the desired control correctionŴ (t). We address this using the strategy described recently in Ref. [22], whereŴ (t) is found perturbatively using a Magnus expansion [23,24]. A major limitation of this approach is that it often requires terms inŴ (t) that are incompatible with the physical system at hand (e.g. interaction terms that do not exist, or that cannot be made time-dependent in the given experimental platform). This is where the present work makes a substantial contribution. We introduce a novel way to find terms in the series expansion ofŴ (t) that are always compatible with all constraints. We achieve this by expandingŴ (t) at each order as a finite sum of timedependent basis functions multiplied by free weights. Finding the required control corrections then amounts in most cases to solving time-independent linear equations for these weights.
As we demonstrate through several examples, this methodology is both extremely flexible and effective; it can also work in systems with many degrees of freedom. The examples we consider include the strong non-RWA driving of a qubit (Sec. III A), leakage errors in a superconducting qubit (Sec. III C), rapid squeezing generation in a parametrically driven bosonic mode (Sec. III B), and accelerated SNAP gates [25,26] in a coupled transmoncavity system (Sec. III D).
Note that the general idea of looking for control fields represented as a finite combination of basis functions was previously used in Refs. [27,28] to design two-qubit superconducting qubit gates that minimize leakage errors. In contrast to those works, our work is both more general and more systematic. Our approach is also complementary to a variational approach for approximately finding STA protocols in complex systems that are compatible with experimental constraints [29,30].

A. Imperfect Unitary Evolution
We consider the generic Hamiltonian: The HamiltonianĤ 0 (t) generates the desired time evolution, whileV (t) is the spurious "error" Hamiltonian that disrupts the ideal dynamics and which can be treated as a perturbation. The perturbative character ofV (t) can originate, e.g., fromV (t) being proportional to a parameter 1, or becauseV (t) is a fast oscillating function. In Appendix A we show why nonresonant error-Hamiltonians can also be corrected with the method presented below. In this section, however, we consider the situation whereV (t) is proportional to a parameter 1 simply because this allows one to count the orders of the perturbative series in a straightforward way. We stress, however, that one can apply the method that we are about to introduce independently of the reason that makesV (t) a perturbation.
The time evolution operator generated byĤ(t) is given byÛ HereÛ 0 (t) represents the ideal time evolution generated byĤ 0 (t) ( = 1), whereT is the time ordering operator, and we assume that the time evolution starts at t = 0. The effect of the error HamiltonianV (t) on the dynamics is given bŷ U I (t), which is defined aŝ Here, an operatorÔ(t) in the interaction picture is given byÔ I (t) =Û † 0 (t)Ô(t)Û 0 (t). Our goal is to have the time evolution operator at t = t f match a specific desired unitary operatorÛ G ; the form of the time evolution operator at earlier times is not relevant for us. This is the case in many problems, the most prominent example being the engineering of quantum gates. We also assume thatĤ 0 (t) provides us the desired time evolution at t = t f , i.e.Û 0 (t f ) =Û G . Consequently, the presence of a non-zero error Hamiltonian V (t) disrupts the evolution and prevents us to generate the desired evolution, since in generalÛ

B. General Strategy to Correct Unitary Evolution
To obtain the ideal unitary evolution at t = t f , we wish to modify the time-dependence ofĤ(t) to cancel the deleterious effects ofV (t). This is formally accomplished by introducing the modified Hamiltonian Here,Ŵ (t) is an unknown control Hamiltonian that cancels, or at least mitigates, the effects ofV (t) on the dynamics, bringing us closer to the desired time evolution [see Fig. 1 (a)]. The unitary evolution generated bŷ H mod (t) is given byÛ mod (t) =Û 0 (t)Û mod,I (t), wherê is the unitary evolution operator generated by the modified Hamiltonian in the interaction picture with respect toĤ 0 (t). We havê The desired unitary operator at t = t f is achieved if A trivial solution to this problem is to takeŴ (t) = − V (t). This solution is almost always infeasible, as the general form ofŴ (t) will be constrained by the kinds of interactions available in the system and their tunability. Furthermore, we are only interested in generating the correct unitary at t = t f and consequently cancelling the spurious Hamiltonian at all times is in some sense demanding more than it is required. A better solution was found in Ref. [22], where one makes use of the fact that the time evolution at intermediate times is not important. This leads to relatively lax conditions that the control HamiltonianŴ (t) must satisfy. Nevertheless, finding an exactŴ (t) is a complex task and generally one needs to resort to perturbation theory to find approximated solutions.
Let us start by writingŴ (t) as a series in , In order to findŴ (t), one could work with the series expansion of the time-ordered exponential of Eq. (6), but a more convenient approach is to use the Magnus expansion [23,24]. With the Magnus expansion we can convert the complicated time-ordered exponential to a simple exponential of an operator that can be expanded in a series: The terms of the Magnus expansion,Ω k (t), are recursively defined by differential equations [23,24], with the first two terms being given by (see also Appendix B) In order to correct the dynamics up to order O( m ), one needs to find a control HamiltonianŴ (t) such that Ω k (t f ) = 0 for k = 1, . . . , m. As shown in Ref. [22] this is accomplished if one firstly truncates the series represent-ingŴ (t) [see Eq. (8)] up to order m and then requires the operatorsŴ (n) I (t) to satisfy the following equations: whereΩ (n) k (t) is the kth term of the Magnus expansion associated to the partially-corrected Hamiltonian Here, the series representing the correctionŴ (t) has been truncated at order n.
To first order (k = 1), Eq. (12) reduces to Equation (12) is the only restriction on the terms of the control HamiltonianŴ (t). This implies we have considerable latitude in how we make our specific choice ofŴ (t). In what follows we fully exploit this freedom to systematically find control Hamiltonians that are completely compatible with experimental constraints on kinds and tunability of available interactions.

C. Constrained Control Hamiltonians
To proceed, we introduce a set of N op timeindependent Hermitian operators {Â j } that form a basis forĤ 0 (t),V (t), andŴ (t). By this, we mean that these operators allow for a unique decomposition of the different Hamiltonian operators at each instant of time: Here h j (t), v j (t), and w j (t) are the real control fields (expansion coefficients) associated with the decomposition ofĤ 0 (t),V (t), andŴ (t), respectively. For instance, the elements of the set {Â j } for a two-level system are the Pauli operatorsσ j with j ∈ {1, 3}. We also introduce the Lie algebra g generated by the set of operators {−iÂ j } with the Lie bracket given by the commutation operation. Having a Lie algebra ensures that one can use the basis formed by the set {Â j } to decompose the operators generated by the Magnus expansion. Finally, we stress that N op can be finite even if the dimension of the Hilbert space is infinite. This is the case for quadratic bosonic forms that can be characterized by the special unitary groups SU(2) or SU(1, 1), which are associated to the Lie algebras su(2) or su(1, 1) [31]. Transforming Eqs. (16) and (17) to the interaction picture defined byĤ 0 (t), we havê Using the fact that {Â j } forms a basis, we can writê Here, the functions a j,l (t) fully encode the action of the interaction picture transformation on our basis operators. Substituting Eq. (19) in Eq. (18), we obtain where we use tildes to denote control fields in the interaction picture, and we havẽ Proceeding analogously forŴ (t) and using the series representation defined in Eq. (8), we get W (n) withw (n) We now return to the fundamental equations of our approach, Eqs. (12), which need to be satisfied to cancel the effects ofV (t) to the desired order. As written, these equations do not reflect any information about relevant experimental constraints. Typical examples of constraints are the inability to control the fields that couple to certainÂ j , i.e. that particular field has to obey w     j (t) that obey Eqs. (12) and simultaneously fulfill the previously mentioned constraints. This then enables the design of high fidelity control pulses that are fully compatible with experimental constraints. As we discuss below, it is enough to show how one derives equations for the first-order control fields w  We proceed by substituting Eqs. (20) and (22) into Eq. (14), which determines the first-order correction Hamiltonian. We obtain an operator equation which can be split into N op equations, one for each operatorÂ j : We stress that Eq. (24) may be ill-defined since it is possible to havew (1) j (t) = 0 whileṽ j (t) = 0 for certain values of j. We show in Sec. III D how to deal with such situations for a large class of problems. For the remainder of this section, we focus on the simpler case where we have a well-defined system of equations.
The problem still remains of how to solve for w (1) j (t); this is still a complex task since one is dealing with a system of N op coupled integral equations. This problem can be overcome by choosing an appropriate parametrization for the functions w with ω k = 2πk/t f and d (1) j0 = 0. This parametrization allows us to carry out the time integration over the duration of the protocol and use the Fourier coefficients as the free parameters to satisfy the system of equations given by Eq. (24). We stress that at this stage finding the first order correction that fulfills Eq. (14) has been reduced to determining a set of N coeffs = N op × (2k max + 1) coefficients. Note that one could use other basis functions for the decomposition, e.g., Slepian functions [32].
The sum in Eq. (25) runs from 0 to k max which allows us to limit the bandwidth of the field associated tô A j . We also note that k max can take different values for different values of j. For constrained systems where a particular field w j (t) must be time independent, we set all the coefficients in Eq. (25) to zero with the exception of c j0 . If one requires w    (25) runs from 0 to k max , but the more general case where the summation runs from k min to k max is also allowed.
We now can formulate the final basic equations of our approach. We substitute Eqs. (23) and (25) in the system of equations defined by Eq. (24). Since we know the explicit time dependence ofw (1) j (t), we can perform the time integration. This leads to a system of timeindependent N op linear equations than can be written in matrix form: Here, x (1) is a vector of coefficients (length N coeffs ) determining the first order control correction that we are trying to find. In contrast, the matrix M and the vector y (1) are known quantities: y (1) parameterizes the error HamiltonianV (t), whereas M encodes the dynamics of the ideal evolution generated byĤ 0 (t).
To be more explicit, the y (1) is a vector of length N op whose components are the spurious error-Hamiltonian elements we wish to average out, x (1) is the vector of the N coeffs unknown Fourier coefficients c (1) lk and d (1) lk that determine our control corrections, c.f. Eq. (25). We order these as follows where j 0 = N op (k max + 1), and the indices l and k in Eq. (28) are functions of j. We have and Here, a//b denotes the integer division of a by b and a%b denotes the remainder of the integer division of a by b.
Finally, M is a (N op ×N coeffs ) matrix that characterizes the evolution under the ideal HamiltonianĤ 0 (t). Recall that the interaction picture transformation generated by this Hamiltonian is described by the functions a jl (t). The matrix elements of M involve the Fourier series of these functions [see Eq. (19)]: where l and k are given by Eqs. (29) and (30), respectively. We stress that Eqs. (27) to (31) are valid when the summation in Eq. (25) runs from 0 to k max for all values of j, but they can be modified to describe other cases.
Higher orders controls are found with an identical procedure. Ultimately, each order is found by solving a system of time-independent N op linear equations similar to Eq. (26) [see Appendix C].
In principle, a set of constrained controls fields that allows one to correct the dynamics up to order n, does not necessarily allows one to correct the dynamics up to order n + 1. In such situations, namely when the obtained linear system does not have a solution [33], one usually has to choose another correction Hamiltonian for the system. There are cases, however, where an alternative solution can be found. We illustrate this situation when we discuss the SNAP problem in Sec. III D.

OI
Operator in the interaction picture with respect toĤ0(t) - jth Magnus operator associated with the modified Hamiltonian Eq. (13)

Aj
Basis operator of the Hilbert space - Eq. (25) Table I. Definition of the most important symbols.

III. APPLICATIONS
In this section, we apply our general strategy to several experimentally relevant problems. These examples highlight the fact that our method is broadly applicable (without modification) to a wide range of very diverse problems.

A. Strong Driving of a Two-Level System
As a first example we consider the problem of a twolevel system (qubit) in the strong driving limit. As we discuss below, this regime generates a complex dynamics that renders precise control of the qubit hard to achieve. Several techniques were used to predict control schemes which generate high-fidelity gates. Optimal control methods have been used, but the resulting control fields are not bandwidth limited and cannot be accurately reproduced by an arbitrary wave form generator [34]. An ad hoc method based on time optimal control of a twolevel system [35,36] was also proposed: it consists in realizing Bang-Bang control with imperfect square control fields [37]. However, to achieve a gate with a reasonably low error the imperfect square pulse must still have a relatively large bandwidth. A method based on analyzing the dynamics of the system using Floquet theory has also been put forward [38,39], but this transforms a low dimensional control problem into a high dimensional one.
The Hamiltonian of a driven two-level system is given byĤ where ω q is the qubit splitting frequency, ω d is the driving Envelopes Gate time ! q t f Coe↵.
α,2 and ∆ = ∆ (1) + ∆ (2) . (c-d) Original pulse and corrected pulse for frequency, f q (t) is the driving envelope, and we introduce the Pauli operators: We label by |0 and |1 the ground and excited states of the system, respectively. We note that the Pauli operators (multiplied by the imaginary number −i) define a Lie algebra with respect to the commutation operation.
In the weak driving limit, i.e. f q (t) ω d ∀t, Eq. (32) allows one to generate rotations around the x-axis if one sets ω d = ω q . This is best understood in the frame rotating at the drive frequency [40]. In this frame the Hamiltonian is given byĤ andV The coefficients v q,j (t) are given by Here, the driving is set on resonance with the qubit frequency, i.e., ω q = ω d . If the system is in the weak driving limit, the fast oscillating terms (also known as counterrotating terms) inV q (t) can be neglected as they average themselves out over the long evolution time set by the slow varying envelope function f q (t). As a consequence, one can approximateĤ R (t) byĤ q,0 (t). This is known as the rotating wave approximation (RWA). The resulting Hamiltonian generates a rotation of angle θ(t f ) around the x-axis, where we have introduced However, when one deviates from the weak driving limit, the counter-rotating terms cannot be neglected anymore since they do not average themselves out on short evolution times. As a result, the dynamics generated byĤ R (t) describes a complex rotation around a time-dependent axis evolving in the xy-plane of an angle which is no more accurately described by Eq. (37) [41] [see Fig. 2 (e)]. To this day there is no known exact solution to this problem, which makes the strong-driving limit impractical to control a qubit with high-fidelity. However, using the general framework laid out in Sec. II, we can mitigate the effects ofV q (t) in situations where the RWA breaks down. This allows us to generate any high-fidelity single-qubit gate beyond the RWA regime. Given the constraints of the original problem, i.e., we only have temporal control over a field coupling toσ x [see Eq. (32)], we look for a correction of the form  x (t) and g y (t) are unknown envelope functions. In addition to the driving field, we also have the liberty to choose the driving frequency; nothing tells us that having ω d = ω q is the best thing to do in terms of control beyond the RWA. In the rotating frame, this is equivalent to have a non-zero detuning ∆ = ω q −ω d . Therefore, we consider the following modified Hamiltonian in the rotating framê In terms of the Pauli operators,Ŵ (n) q (t) is given bŷ with In practice, having a control field with two-quadratures driving [see Eq. (38)] and introducing a detuning has given us the ability to implement 3-axes control. We stress that there are other possible choices forŴ (t), but they all require more resources to be implemented experimentally [see Appendix D 1]. Note that the modified detuning is given by ∆ = n ∆ (n) in complete analogy to having the control fields represented by a series [see Eq. (8)].
Following the general strategy presented in Sec. II, we first move to the interaction picture with respect tô In Eqs. (43) and (45), we have omitted the explicit time dependence of θ for simplicity, i.e., θ = θ(t) [see Eq. (37)]. The next step consists in expanding the control fieldsw (n) q,j (t) (j ∈ {x, y, z}) [see Eq. (45)] into a Fourier series. However, before proceeding it is useful to notice the special form of the functionsw (n) q,j (t): an unknown function that multiplies a known fast oscillating function. It is therefore more suitable to just expand the unknown functions g (n) y (t) [see Eq. (38)], and ∆ (n) in a Fourier series and use the corresponding Fourier coefficients as the free parameters to satisfy the system of equations generated by the Magnus-based approach. We stress, however, that one obtains exactly the same results using the general procedure of Sec. II and imposing the necessary constraints on the Fourier series.
If we constrain g α=x,y (t) to be zero at t = 0 and t = t f , which is often the case experimentally, we obtain the following Fourier expansions (46) and (47) where ω k = 2π/t f . Since we have a total of three equations of the form of Eq. (24) to solve (one for each Pauli operator), we need at least three free parameters. Consequently, we can set all coefficients to zero in Eqs. (46) and (47) y,1 and c (n) z,0 [42]. With this choice, Eqs. (46) and (47) reduce to and The final step is to find the value of the free parameters c (n) x,1 , c (n) y,1 and ∆ (n) . We start by verifying that by substituting Eqs. (48) and (49) in Eq. (44), we obtain a correction HamiltonianŴ where x x,1 , c y,1 , ∆ (1) } T is the vector of unknown coefficients [see Eq. (28)], y ,ṽ q,y (t),ṽ q,z (t)} T is the vector of the spurious error-Hamiltonian elements with v q,j (t) (j ∈ {x, y, z}) defined in Eq. (43), and P q [43] is the matrix that characterizes the evolution under the ideal HamiltonianĤ q,0 (t) [see Eq. (34)]. The explicit matrix elements of P q can be found in Appendix D 2. Higher-order correction Hamiltonians can be found in a similar way.
In Fig. 2 (a), we plot the average fidelity error ε [44] for a Hadamard gate generated with an initial envelope with θ 0 = π/2. Other gates can be realized by choosing θ 0 ∈ [0, 2π]. The blue trace shows the error for the uncorrected evolution while the green trace shows the error of the corrected evolution up to second order. The latter, as one can observe in Fig. 2 (a), globally increases when ω q t f decreases, but around ω q t f 1 the error of the corrected evolution starts decreasing again. This can be understood by considering the limit t f → 0 (ω q t f → 0). In this limit, we haveṽ q, (43)] which implies thatV I (t) commutes with itself at all times. As a consequence, one can find exact modifications to the control fields since only the first order of the Magnus expansion is non-zero. However, as one can see in Fig. 2 (b), where we plot the coefficients of the correction versus the gate time t f , the modified control sequences require control fields with diverging amplitudes. Restricting ourselves to gate times close to unity (ω q t f 1), where the modified control sequences can be experimentally realized, our strategy improves the error ε by more than two orders of magnitude. In Figs. 2 (c) and (d), we compare the original and corrected pulses for ω q t f ≈ 5. One can observe that the changes to the original pulse are small. For convenience we write the nth order modified pulse as where f (n) q,x (t) = f q (t) + n k=1 g (n) x (t) and f As a second example, we consider the problem of fast generation of squeezed states using a parametrically driven cavity (PDC). The ability to generate squeezed states with quantum oscillators is of particular interest since it allows one, among others, to enhance sensing capabilities [45] or to reach the single-photon strong coupling regime with optomechanical systems using only linear resources [46]. Recently, optimal control techniques have been used to achieve squeezing of an optomechanical oscillator at finite temperature [47].
Here, we are interested in generating squeezing on a relatively short time scale by using a pulsed drive. As for the qubit problem discussed in Sec. III A, this turns out to be a complex task due to fast counter-rotating terms that prevent the preparation of the desired squeezed state.
The Hamiltonian of a PDC corresponds to having a harmonic oscillator with a modulated spring constant. corr. RWA uncorr.

7 11
Pulse width ! a t f 10 5 This can be achieved, e.g., in the microwave regime by modulating the magnetic flux through a SQUID loop (flux-pumped Josephson parametric amplifier) [48,49]. We havê withâ (â † ) the bosonic annihilation (creation) operator. The frequency of the modeâ is ω a and the drive has frequency ω d . It is convenient to introduce the operators [31] which define (multiplied by the imaginary number −i) a Lie algebra with respect to the commutation operation [see Appendix E 1]. As mentioned earlier, since the Hamiltonian is quadratic, the three operators defined in Eq. (54) are enough to completely describe the full dynamics in spite of having an infinite Hilbert space. The action of these operators is best understood in the phase space defined byx = (â +â † )/ √ 2 and y = −i(â −â † )/ √ 2:μ x generates squeezing along the x-axis,μ y generates squeezing along the y-axis, andμ z generates a rotation around the origin of the phase space.
In a frame rotating at a frequency ω d /2 = ω a , the Hamiltonian becomesĤ D,R (t) =Ĥ D,0 (t) +V D (t) witĥ and ω d ∀t. This results inĤ D,R (t) ≈Ĥ D,0 (t) and the generated dynamics corresponds to squeezing along the y-axis with a degree of squeezing depending on r(t f ), with As one deviates from the weak driving limit,V D (t) cannot be neglected anymore. The generated dynamics becomes then more complex with the counter-rotating terms changing the direction along which the squeezing is generated as well as degrading the final degree of squeezing [see Fig. 3

(d)].
To mitigate the effects ofV D (t) [see Eq. (56)], we consider a control Hamiltonian that corresponds to just changing the initial form of the parametric modulation. This leads to the correction Hamiltonian W PDC (t) = n g (n) x (t) cos(ω d t) + g (n) y (t) sin(ω d t) × â +â † 2 .
(58) Furthermore, we are at liberty to drive the PDC at a frequency that is detuned from that of modeâ, with ∆ = n ∆ (n) a static detuning. Following the general procedure introduced in Sec. II (see Appendix E 2), we can easily determine ∆ (n) , g (n) x (t) and g (n) y (t). We stress that in this example we correct the unitary evolution generated by Eq. (53), which allows us to generate the ideal squeezing dynamics for any initial state. This is in contrast to optimizing the dynamics to get optimal squeezing of the vacuum state only.
In Fig. 3 (a), we plot the degree of squeezing S as a function of the total evolution time t f for the RWA (red trace), the uncorrected (blue trace), and the corrected (green trace) evolutions. The degree of squeezing is given by whereŷ = (â −â † )/i √ 2, and ŷ i,f = ψ i,f |ŷ|ψ i,f is the quantum average of the operatorŷ with respect to the initial and final states, respectively. Here, the initial state is the vacuum state |0 . The initial pulse envelope is given by Within the RWA the degree of squeezing is independent of the pulse width t f , since the squeezing depends just on r(t f ) = c. In the regime where the fast oscillating terms cannot be neglected, it is clear that the corrected evolution gives substantially better results (closer to the RWA evolution), specially for small values of t f . In Fig. 3 (c), we compute the deviation angle ϕ in the phase space (with respect to the y-axis) where the maximum squeezing is obtained. Ideally, the maximum squeezing should be in the direction of the y-axis and ϕ should be zero.
With the correction Hamiltonian ϕ is much closer to the ideal value. In Fig. 3 (b), we plot the coefficients of the correction Hamiltonian as a function of the total evolution time t f . As for the qubit case, we observe that the modified control fields can be seen as adding a small correction to the original control fields.

C. Transmon Qubit
As a next example we consider the problem of realizing single-qubit gates with a transmon qubit [50], where the logical qubit states are encoded in the two lowest energy states of an anharmonic oscillator with eigenstates |n [see Fig. 4 (c)]. Since the oscillator is only weakly anharmonic, driving the |0 ↔ |1 transition unavoidably leads to transitions to higher energy states outside of the computational subspace (leakage). Several strategies have been put forward to suppress leakage while implementing a gate, with perhaps the most well-known approach being DRAG (Derivative Removal by Adiabatic Gate) [13,51]. However, the correction predicted by DRAG cannot be fully implemented experimentally as it also requires one to drive the |0 ↔ |2 transition. There is no charge matrix element connecting these states, hence it cannot be driven by an extra tone at the transition frequency. While neglecting this unrealizable control field is the simplest thing to do, this is a somewhat uncontrolled approximation; further, it has been demonstrated experimentally [52] and theoretically [22] that this is indeed not the optimal approach. In the rest of this section, we demonstrate how our general strategy allows one to systematically find control sequences that are fully compatible with the constraints of the problem (i.e. no direct |0 ↔ |2 drive, no time-dependent detuning), and also are highly efficient in suppressing both leakage and phase errors.
As in the original DRAG paper, we consider the three-level Hamiltonian as an approximation of the weakly anharmonic oscillator.
Here, ω T is the frequency splitting between the energy levels |0 and |1 while the frequency splitting between |1 and |2 is given by ω T + α, where α is the anharmonicity. We have also defined the operatorŝ  (17)]. This set of eight operators (multiplied by the imaginary number −i) also form a Lie algebra with respect to the commutation operation, thus this set of eight operators can also be used to uniquely decompose the operators generated by the Magnus expansion. The control pulse consists of a drive at frequency ω d and an envelope function f T (t). As one can see from Eq. (62), driving the |0 ↔ |1 transition also results in the |1 ↔ |2 being driven with a relative strength given by η, which unavoidably generates leakage out the qubit subspace.
In a frame rotating with frequency ω d , the Hamiltonian is given byĤ T (t) =Ĥ T,0 (t) +V T (t), wherê andV Here, we assume that the driving is on resonance with the |0 ↔ |1 transition, i.e., ω T = ω d . The Hamilto-nianĤ T,0 (t) gives us the desired interaction: it couples the levels |0 and |1 , allowing one to perform unitary operations in the computational space, while leaving the level |2 isolated. The HamiltonianV T (t) couples levels |1 and |2 leading to leakage out of the computational subspace. Note that we have neglected the terms oscillating at frequencies close to 2ω d in Eqs. (64) and (65) (RWA) [53]. Given the constraints of the problem [see Eq. (64)], we want to find a correction that only involves modifying the drive envelope we use, and possibly changing the detuning in a static manner. We thus write the control Hamiltonian (in the rotating frame) aŝ W TLS (t) = n g (n) x (t) cos(ω d t) + g (n) y (t) sin(ω d t) with g (n) x (t) and g y (t) the unknown envelope functions. Furthermore, we allow the drive frequency to be detuned with respect to the base frequency of the transmon, As for the envelope functions, the detuning is parametrized as a series: ∆ = n ∆ (n) . Within our framework, we would in principle need a total of eight free parameters to satisfy Eqs. (24), which determine the first-order correction; this is because there are eight operators in the basis. Taking into account that |2 , which is outside the computational space, is of no interest to us, the equation associated to the operator |2 2| can be neglected. More generally, the equations originating from operatorsÂ j that act strictly outside of the computational space do not need to be fulfilled, and one can simply neglect them to arrive at the relevant system of equations for the given order.
We are therefore left with seven equations to fulfill, and we need at least seven coefficients. However, we cannot select seven coefficients at random: we must retain enough harmonics to ensure that g (n) x (t) and g (n) y (t) have a bandwidth comparable to |α|; this is needed to easily correct leakage transitions. This was also identified in an earlier work by Schutjens et al. [54], which also aims at finding modified pulses to mitigate leakage errors in a transmon. Their strategy consists in suppressing the spectral weight associated to leakage transitions from the control fields. A systematic way of ensuring that the control fields have enough bandwidth is to keep more non-zero coefficients than necessary in their Fourier expansion. This choice leads to an underdetermined linear system of equations which can be solved using the Moore-Penrose pseudo-inverse [55][56][57] (see Appendix F).
To show the performance of our strategy, we considered the situation where one wants to perform a Hadamard gate in the computational subspace. In Fig. 4 (a), we plot the average fidelity error as a function of the gate time t f . We compare the results obtained in the absence of any correction (blue trace) with the results for a 2nd order Magnus-based correction (green trace), a 6th order Magnus-based correction (red trace), and the DRAG correction (purple trace) [13]. The results show that the 6th order Magnus correction reduces the average fidelity error by more than four orders of magnitude for small |αt f |, greatly outperforming the DRAG correction. In Fig. 4 (b) we compare the original and modified pulses for |α|t f = 5. For convenience we write the nth order modified pulse as x (t) and f  Envelopes T,y (t) f (6) T,x (t) f (6) T,y (t) (c) formalism, since arbitrary waveform generators (AWG) have bandwidth limitations. We remind the reader, however, that our method allows direct control over the bandwidth of the pulse through truncation of the Fourier series. If a stricter limitation over the bandwidth of the correction pulse is needed, one can make use of Lagrange multipliers to look for solutions of the linear system. As a rule of thumb, the minimum requirement of our method is that the AWG bandwidth should approximately be comparable to or larger than the anharmonicity |α|.

D. SNAP Gates
We now turn to an example that combines both qubit and bosonic degrees of freedom. The general problem is to use a qubit coupled dispersively to a cavity to achieve control over the bosonic cavity mode. A method for doing this was recently proposed and implemented experimentally in a superconducting circuit QED architecture: the so-called SNAP gates (selective numberdependent arbitrary phase gates) combined with cavity displacements [25,26]. Our goal will be to use our general method to accelerate SNAP gates without degrading their overall fidelity.
An optimal control approach based on GRAPE has been used to accelerate the manipulation of the bosonic cavity mode [58]. There is, however, a major advantage in using SNAP gates in combination with cavity displacements: the SNAP gate can be made robust against qubit errors [59], i.e., noise acting on the qubit will not affect the quantum state of the cavity.
As we will see, this problem involves an interesting technical subtlety. When introducing our general method in Sec. II, we stressed that it is crucial for the Hamil-tonianŴ I (t) describing the modification of the control fields to have terms involving all of the basis operatorŝ A j appearing in the Magnus expansion of the unitary evolution generated by the error-HamiltonianV I (t). If this was not true, it would seemingly be impossible to correct errors proportional to these basis operators. Surprisingly, there are cases where this conclusion is overly pessimistic. In certain cases, one can still use a modified version of our Magnus-based strategy which uses an alternate method for finding an appropriateŴ (t). As we show below, correcting SNAP gates is an example of this kind of situation. The general price we pay is that now, to find an appropriate set of control corrections, we need to solve a nonlinear set of equations (instead of the linear equations in Eq. (26) that we used in all the previous examples).
The basic setup for SNAP gates involves a driven qubit that is dispersively coupled to a cavity mode. The Hamiltonian isĤ SNAP (t) =Ĥ qc +Ĥ D (t), witĥ andĤ The Pauli operatorsσ α act on the Hilbert space of the qubit and have been defined in Eq. (33). We also introduce the annihilation (creation) operatorâ (â † ) destroying (creating) an excitation of the oscillator. The qubit is driven by two independent pulses, f x (t) and f y (t), which couple both toσ x with the same frequency ω d but with different phases.
In the interaction picture with respect toĤ qc , the Hamiltonian becomeŝ where δω n = ω q + χn − ω d , |n is a bosonic number state, and we have neglected fast oscillating terms. If the drive is now chosen to fulfil ω d = ω q + χn 0 , so that the drive is resonant for a particular number-selected qubit transition, the Hamiltonian defined in Eq. (71) can be written asĤ S (t) =Ĥ S,0 (t) +V S (t). Herê is the resonant part of the Hamiltonian defined in Eq. (71) and allows one to generate a unitary operation in the subspace spanned by {|g, n 0 , |e, n 0 }. In contrast − f y (t) [sin(δω n t)σ x + cos(δω n t)σ y ] |n n| (73) is the non-resonant part of Eq. (71). This error Hamiltonian is responsible for the unwanted dynamics in the subspace spanned by {|g, n , |e, n }, for n = n 0 . While in principle the effects ofV S (t) on the dynamics cannot be avoided, they are minimal in the weak-driving regime where f x (t), f y (t) χ. In this limit, we can useĤ S,0 (t) to generate a dynamics that imprints a phase on |n 0 while leaving all other states |n (n = n 0 ) unchanged. Our general goal will be to relax this weak-driving constraint, allowing for a faster overall gate.
For concreteness, we assume that the qubit is initially in the state |g and the driving pulses f x (t) and f y (t) are chosen such that the qubit undergoes a cyclic evolution, i.e., the trajectory on the Bloch sphere encloses a finite solid angle and at t = t f the state of the qubit is back to |g . This leads to the accumulation of a Berry phase γ at t = t f for the qubit which conditioned on the state of the cavity being |n 0 . In other words, is the unitary evolution generated by the ideal Hamiltonian in Eq. (72). This approach can be generalized so that the ideal evolution yields different qubit phase shifts for a set of different cavity photon numbers. One simply replaces the driving where ω d,n = ω q + χn. The pulse envelopes f x,n (t) and f y,n (t) are chosen such that one gets the desired phase in the nth energy level.
Of course, the above ideal evolution requires that f x,n (t), f y,n (t) χ, constraining the overall speed of the gate. Without this assumption, the effects of the offresonant error interaction given by the generalization of V S (t) [c.f. Eq. (73)] cannot be neglected, and will compromise the ideal SNAP gate evolution. Again, our goal is to mitigate these errors, allowing for faster gates.
In the following, we consider for simplicity the situation where one wants to imprint a phase on a single energy level of the oscillator. The extension to the more general situation where one imprints arbitrary phases in different levels is straightforward. We truncate the bosonic Hilbert space and work only within the subspace formed by the N trunc first number states. This procedure is justified by the fact that SNAP gates are typically used to manipulate "kitten" states [25,26], which are themselves restricted to a truncated subspace of the original bosonic Hilbert space.
As we did for the previous examples, we start by choosing a correction HamiltonianŴ SNAP (t) that one can realize experimentally. Here, this corresponds to a modification of the qubit drive amplitudes: [g x,n (t) cos(ω d,n t) + g y,n (t) sin(ω d,n t)]σ x (76) where ω d,n = ω q + χn. Moving to the interaction picture with respect toĤ qc [see Eq. (69)] and neglecting nonresonant terms, we obtain In the interaction picture defined byĤ S,0 (t) [see Eq. (72)], we find that the form of the non-resonant error Hamiltonian is unchanged: sinceĤ S,0 (t) commutes withV S (t);Ĥ S,0 (t) andV S (t) act on orthogonal subspaces. On the other hand,Ŵ S (t) acts on the whole Hilbert space, and is transformed when moving to the interaction picture. We find: n0 (t)σ x − g y,n0 (t)σ y ]|n 0 n 0 |Û S,0 (t) (1 − δ n,n0 )[g x,n (t)σ x − g y,n (t)σ y ]|n n|.
(79) The first term of Eq. (79) acts on the {|g, n 0 , |e, n 0 } subspace only and has terms proportional to all three Pauli matrices. While the explicit expression is too lengthy to be displayed here, it can be readily found using the group properties of the Pauli operators. The second term, which acts on the orthogonal subspace, has only terms proportional toσ x |n n| andσ y |n n|. This means that the correction Hamiltonian in Eq. that we have used successfully in all of the previous examples.
The naive thing to do would be to find an alternative control Hamiltonian that directly provides terms proportional toσ z |n n| in the interaction picture. However, in the lab frame this translates into a Hamiltonian with a dispersive coupling constant dependent on photon number n, i.e., we would need a term n χ n |n n| in Eq. (69). This is extremely difficult to achieve experimentally, hence we do not pursue this approach further.
A more promising approach is to use the fact that even though our original (constrained) correction Hamiltonian W S,I (t) is missing important terms, these can nonetheless be dynamically generated. In the same way thatV S,I (t) generates problematic terms proportional toσ z |n n| at second order in the Magnus expansion, so canŴ S,I (t). Thus, we look for a correction HamiltonianŴ (1) (t) that cancels the sum of the first two terms of the Magnus expansion:Ω (1) This is in contrast of the general strategy in Sec. II, where we would just cancel the first term in the Magnus expansion. We can use Eqs. (10) and (11)  . The explicit equation can be found in the Appendix G 1. Here, to be able to correct for the term proportional toσ z |n n|, we cannot discard the higher order term generated byŴ (1) S,I (t) inΩ (1) 2 (t f ) (as the standard procedure prescribes). Once this is done, we proceed as usual: we expand the pulse envelopes g x,n (t) and g y,n (t) in a Fourier series, and we truncate the series keeping a sufficiently large number of free parameters [60]. Following the strategy presented in Sec. II, we derive a system of equations for the free parameters, but instead of obtaining a linear system of equations we get a system of quadratic equations for the free parameters, since Eq. (80) is quadratic inŴ (1) S,I (t). This system of equations can be solved numerically; see Appendix G 1, where we also show how one can find higher order corrections.
In the situation where one wishes to imprint non-zero phases to all energy levels of the truncated Hilbert space, one can actually solve the problem following the general (linear) strategy shown in Sec. II. In this case, since one is driving all frequencies resonantly, the ideal unitarŷ U S,0 (t) acts on the whole truncated Hilbert space of the cavity. As a consequence, transforming the correction HamiltonianŴ S (t) [see Eq. (77)] to the interaction picture will generate terms proportional toσ z |n n| for all values of n (see Appendix G 2). We stress, however, that SNAP gates are most often used to manipulate logical qubit states encoded in a finite superposition of same parity bosonic number states [61], e.g., |0 L = (|0 + |4 )/ √ 2 and |1 L = |2 . Accelerating SNAP gates that act on such logical qubit states requires one to use the strategy that cancels the sum of the first terms of the Magnus expansion [see Eq. (80)].
In Fig. 5 (a), we show the fidelity error when one tries to implement a fast SNAP gate that imprints a π/2 phase in the cavity energy levels |0 and |4 simultaneously. This is similar to implement a Z-gate for a logical qubit encoded in the states |0 L = (|0 +|4 )/ √ 2 and |1 L = |2 . The envelope functions for n = 0 and n = 4 are given by and f y,n (t) = For any other values of n, we have f x,n (t) = f y,n (t) = 0. Here, we kept only ten energy levels for the cavity, i.e., the highest bosonic number state is |10 , and the fidelity error was calculated using only the states within the truncated Hilbert space. We have plotted the fidelity as a function of gate time for the unmodified Hamiltonian (blue trace), for the second-order (green trace), and for the fourth order (red trace) modified Hamiltonians. Since we are only manipulating the cavity energy levels |0 and |4 , we need to use the strategy that cancels the sum of the first terms of the Magnus expansion [see Eq. (80)]. The fourth order modified Hamiltonian achieves fidelity errors that are at least one order of magnitude smaller than the fidelity error of the original Hamiltonian. For larger values of t f the difference can reach almost four orders of magnitude. In Fig. 5 (b) and (c) we show the spectrum of the original and modified pulses for a gate time of χt f = 50. The original pulse has only peaks located at ω = 0 and ω = 4χ, since these are the frequencies of the levels being driven. The modified pulse, however, has peaks located at frequencies ω = 0, χ, 2χ, . . . , 9χ. This shows that the corrected pulse undoes residual rotations caused by the non-resonant interaction in the different bosonic number state subspaces in order to bring the final state close to the target state. It is important to note that the modified pulse corrects the dynamics only within the truncated Hilbert space. If the initial state of the cavity, i.e. the state before the SNAP operation is performed, is not confined to the truncated Hilbert space, the corrected pulse will not bring any improvement in terms of fidelity error, since the states lying outside the truncated Hilbert space will still be affected by the correction pulse.

IV. CONCLUSION
We have developed a method that allows one to design high-fidelity control protocols which are always fully compatible with experimental constraints (available interactions and their tunability, bandwidth, etc.). At its core, our method uses the analytic solution of a simple control problem as a starting point to solve perturbatively a more complex problem, for which it is impossible to find closed-form analytic solutions. At the end of the day, the complex control problem is converted into solving a simple linear system of equations. We have applied our method to a range of problems, including the leakage problem in a transmon qubit and SNAP gates. We have shown how the control sequences predicted by our strategy allow one to substantially decrease the error of unitary operations while simultaneously speeding up the time require to complete the protocols. Finally, we note that the protocols generated by our method could be further improved by using them to seed a numerical optimal control algorithm. In Sec. II, where we introduced the general strategy to correct unitary evolution, we have taken as a starting point the Hamiltonian defined in Eq. (1), where the spurious error-HamiltonianV (t) is multiplied by a small parameter , such thatV (t) can be considered as a perturbation. However, in all the examples discussed in the main text, the spurious error-Hamiltonian is not multiplied by a small parameter , but by an oscillating function of time. Since integrating this oscillating function over a sufficiently long time interval leads to self-averaging, one can still view this class of error-Hamiltonians as a perturbation. The Magnus expansion provides a good way to formally show that.
Let us consider the generic oscillating Hamiltonian in the interaction picture, where ω V is the frequency at which the unwanted HamiltonianV I (t) oscillates andq(t/t f ) is a bounded operator.
Considering the Magnus expansion of the evolution operatorÛ , whereT is the time-ordering operator, we have to first order, at t = t f , where Eq. (A3) was obtained by integrating Eq. (A2) by parts. One can perform the integration by parts repeatedly and obtain a series expansion forΩ 1 (t f ) in powers of (ω V t f ) −1 . As one can see from Eq. (A3), the leading order of This suggests that (ω V t f ) −1 plays the role of the small parameter . Taking this analogy one step further, one would naively assume that the leading order ofΩ Substituting Eq. (A1) in Eq. (A4) and performing the integration in t 2 by parts, we obtain which shows thatΩ 2 (t f ) scales to leading order like (ω V t f ) −1 . Nevertheless, since we know the Magnus expansion converges, we must have that the leading order ofΩ n (t) scales with higher powers of (ω V t f ) −1 for increasing n, but as we show above this dependence is not trivial. Numerical tests suggest that for nonresonant perturbations the leading order ofΩ j (t f ) is given bŷ Because of this property of nonresonant perturbations, one often needs to find the correction HamiltonianŴ (t) up to second order to mitigate substantially the errors generated byV (t).
Here j 0 = N op (k max + 1), and l and k are given by the Eqs. (29) and (30) of the main text. The elements of y (n) are given by According to Eq. (11) of the main text, the second term in the Magnus expansion is given by, A natural choice forŴ (1)

TransformingŴ
(1) I (t) back to the original frame, we obtain the correction HamiltonianŴ (t) = W (1) (t) +Ŵ (2) (t). To implement this correction Hamiltonian, one would not only need to control in time the fields coupling toσ x andσ z , but one would also need an extra time-dependent control field that couples toσ y .

Matrix Elements of Pq
In this appendix, we give explicitly the matrix elements of the matrix P q introduced in Eq. (50) of the main text when discussing the strong driving of a qubit. We have where once more we have omitted the explicit time dependence of θ. In this section, we give the commutation relations for the operatorsμ x ,μ y , andμ z introduced in Eq. (54) of the main text when discussing the strong driving of a parametrically driven cavity. These operators behave as generators of the group SU(1, 1) and consequently generate the su(1, 1) Lie algebra, which one can readily verify by computing the commutation relations. We have [μ x ,μ y ] = 2iμ z , [μ x ,μ z ] = 2iμ y , [μ z ,μ y ] = 2iμ x .
(E1) Therefore these three operators are enough to fully characterize the dynamics of the parametrically driven cavity in spite of having an infinite Hilbert space.

Correction Hamiltonian
In this appendix we give some more details about the steps of the general method applied to the problem of strong driving of a parametrically driven cavity.
Following the general procedure described in the general theory (see Sec.II of the main text), we start by writing the full modified Hamiltonian in the frame rotating at the drive frequency ω d : whereĤ D,0 (t) andV D (t) are, respectively, given by Eq. (55) and Eq. (56) of the main text, and W (n) D (t) = 2g (n) (t) cos(ω d t)μ x + 2g (n) (t) sin(ω d t)μ y + [∆ (n) + 2g (n) (t)]μ z .
Once more, we stress that the final detuning is given by ∆ = n ∆ (n) . Following our recipe, we now move to the interaction picture with respect toĤ D,0 (t). The HamiltonianV D (t) is then given byV D,I (t) =ṽ D,x (t)μ x +ṽ D,y (t)μ y +ṽ D,z (t)μ z ,

The Linear System of Equations
In this section we give explicitly the matrix elements of P D . As explained in Sec. III A, the matrix P D is analog to the matrix M introduced in Eq. (26) of the main text and encodes the dynamics of the ideal evolution generated byĤ D,0 (t) (see Eq. (55) of the man text). The difference comes from our choice of expanding in a Fourier series the unknown envelope functions g α (t) (α = x, y, see Eq. (58) of the main text) and the static detuning ∆ (n) instead of the functionsw D,α (t) [α = x, y, z, see Eq. (E7)]. Within this framework, the system of linear equations that allows one to determine the coefficients defining the nth order control correction is where x (n) D = {c (1) x,n , c y,n , ∆ (n) } T is the vector of unknown coefficients and y where ω 1 = 2π/t f and once more we have omitted the explicit time dependence of r.
Hamiltonian in the lab frame is given by Eq. (75) of the main text. Thus, in the interaction picture with respect tô H qc (see Eq. (69) of the main text), the Hamiltonian describing the dispersive coupling between the qubit and cavity is given byĤ S (t) =Ĥ S,0 (t) +V S (t), wherê Considering the same correction Hamiltonian as in the main text [see Eqs. (76) and (77)], we find that in the interaction picture with respect toĤ S,0 (t) [see Eq. (G7)]Ŵ S,I (t) has terms proportional toσ z |n n|. In contrast to the case we considered in the main text, whereĤ S,0 (t) acts only on the subspace of |n 0 (see Eq. (72) of the main text), when we drive resonantly all bosonic number states the HamiltonianĤ S,0 (t) [see Eq. (G7)] acts on the whole (truncated) Hilbert space. This allow us to find a correction Hamiltonian by simply solving a system of linear equations, like we did for the other problems discussed in the Applications sections of the main text.