Engineering Fast High-Fidelity Quantum Operations With Constrained Interactions
Thales Figueiredo Roque Aashish A. Clerk Hugo Ribeiro
quant-ph




Max Planck Institute for the Science of Light, Staudtstraße 2, 91058 Erlangen, Germany


Pritzker School of Molecular  Engineering,  University  of  Chicago, 5640  South  Ellis  Avenue,  Chicago,  Illinois  60637,  U.S.A.


Max Planck Institute for the Science of Light, Staudtstraße 2, 91058 Erlangen, Germany


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.


Engineering Fast High-Fidelity Quantum Operations With Constrained Interactions
    Hugo Ribeiro
    December 30, 2023
===============================================================================





§ 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 <cit.> 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. <cit.>) based on optimal quantum control theory
<cit.>.  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 <cit.>, 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. <cit.>),  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 <cit.>). 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 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. <cit.>, 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_) a desired unitary operation.  We start by splitting the Hamiltonian into two parts as (t) = _0 (t) +
(t), where H_0(t) is simple enough to be analytically tractable, and (t) represents all the additional interactions that
make the problem unsolvable. The basic strategy then has two parts:


    
  * First, choose control fields in the “simple" Hamiltonian _0(t) so that in the absence of (t), one realizes
	  the desired operation.  This can be done analytically.
    
  * Adding back (t) will then destroy the ideal evolution.  We address this by modifying available control fields so as
	  to  average out the impact of (t).  This amounts to adding a control correction to the full Hamiltonian:  (t)
	  →(t) + (t) [see Fig. <ref> (a)]. 


The question is of course how to find the desired control correction (t).  We address this using the strategy described
recently in Ref. <cit.>, where (t) is found perturbatively using a Magnus expansion <cit.>.
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 time-dependent 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. <ref>),
leakage errors in a superconducting qubit (Sec. <ref>), rapid squeezing generation in a parametrically driven bosonic
mode (Sec. <ref>), and accelerated SNAP gates <cit.> in a coupled transmon-cavity system
(Sec. <ref>).      

Note that the general idea of looking for control fields represented as a finite combination of basis functions was previously
used in Refs. <cit.> 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 <cit.>. 



§ THEORETICAL FRAMEWORK





 §.§ Imperfect Unitary Evolution


We consider the generic Hamiltonian:


    (t) = _0 (t) + ϵ (t).


The Hamiltonian _0 (t) generates the desired time evolution, while (t) is the spurious “error" Hamiltonian that disrupts the
ideal dynamics and which can be treated as a perturbation. The perturbative character of (t) can originate, e.g., from (t) being proportional to a parameter ϵ≪ 1, or because (t) is a fast oscillating function. In
Appendix <ref> we show why nonresonant error-Hamiltonians can also be corrected with the method presented
below. In this section, however, we consider the situation where (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 makes (t) a perturbation. 

The time evolution operator generated by (t) is given by


    (t) = _0 (t) _ (t).


Here _0 (t) represents the ideal time evolution generated by _0 (t) (ħ=1),


    _0 (t)=exp[-i ∫_0^t t_1_0(t_1) ],


where  is the time ordering operator, and we assume that the time evolution starts at t=0. The effect of the error
Hamiltonian (t) on the dynamics is given by _ (t), which is defined as 


    _ (t)=exp[-i ϵ∫_0^t t_1_ (t_1) ].


Here, an operator (t) in the interaction picture is given by _ (t) = _0^† (t)  (t) _0 (t). 

Our goal is to have the time evolution operator at t=t_  match a specific desired unitary operator _ ; 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_, i.e.  _0 (t_) = _. Consequently, the presence of a non-zero error Hamiltonian (t) disrupts the
evolution and prevents us to generate the desired evolution, since in general _ (t_) ≠1 [see
Eq. (<ref>)]. 



 §.§ General Strategy to Correct Unitary Evolution


To obtain the ideal unitary evolution at t=t_, we wish to modify the time-dependence of (t) to cancel the deleterious
effects of (t). This is formally accomplished by introducing the modified Hamiltonian


    _mod (t) = _0 (t) + ϵ (t) +  (t).


Here, (t) is an unknown control Hamiltonian that cancels, or at least mitigates, the effects of (t) on the
dynamics, bringing us closer to the desired time evolution [see Fig. <ref> (a)]. The unitary evolution generated by
_mod (t) is given by  _mod (t) = _0 (t) _mod, (t), where 


    _mod, (t) = exp[ -i ∫_0^t t_1_mod, (t_1) ],


is the unitary evolution operator generated by the modified Hamiltonian in the interaction picture with respect to _0 (t). We
have 


    _mod, (t) = ϵ_ (t) + _ (t).


The desired unitary operator at t=t_ is achieved if _mod, (t_) = 1, i.e., _mod (t_) =
_0 (t_) = _. 

A trivial solution to this problem is to take (t) = -ϵ(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_ 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. <cit.>, 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 ϵ,


    (t) = ∑_n=1^∞ϵ^n ^(n) (t).


In order to find (t), one could work with the series expansion of the time-ordered exponential of Eq. (<ref>), but a
more convenient approach is to use the Magnus expansion <cit.>. 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:


    _mod, (t) =  exp[ ∑_n=1^∞_n(t) ].


The terms of the Magnus expansion, _k(t), are recursively defined by differential equations <cit.>,
with the first two terms being given by (see also Appendix <ref>)


    _t Ω̂_1(t)    = -i _mod, (t),
    	
    _t Ω̂_2(t)    = 1/2 [_t _1(t), _1(t)].


In order to correct the dynamics up to order 𝒪(ϵ^m), one needs to find a control Hamiltonian (t) such
that _k (t_) = 0 for k = 1, …, m. As shown in Ref. <cit.> this is accomplished if one
firstly truncates the series representing (t) [see Eq. (<ref>)] up to order m and then requires the operators
^(n)_(t) to satisfy the following equations:


    ϵ^n ∫_0^t_t ^(n)_ (t) = -i ∑ _k=1^n _k^(n-1)(t_),


where _k^(n)(t) is the kth term of the Magnus expansion associated to the partially-corrected Hamiltonian 


    _mod,^(n) (t) = ϵ_(t) + ∑_k=1^n ϵ^k ^(k)_(t).


Here, the series representing the correction (t) has been truncated at order n. 

To first order (k=1), Eq. (<ref>) reduces to 


    ∫_0^t_t ^(1)_ (t) = - ∫_0^t_t _(t).
 

Equation (<ref>) 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.   



 §.§ Constrained Control Hamiltonians


To proceed, we introduce a set of N_op time-independent Hermitian operators {_j} that form a basis for _0 (t),
(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:


    _0(t)    = ∑_j h_j (t) _j,
    	
     (t)    = ∑_j v_j (t) _j,
    	
    (t)    =  ∑_j w_j (t) _j.


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), (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 𝔤 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) <cit.>.

Transforming Eqs. (<ref>) and (<ref>) to the interaction picture defined by _0(t), we have 


    _ (t) = ∑_j v_j(t) _j, .


Using the fact that {_j} forms a basis, we can write 


    _j,  = ∑_l a_j,l(t) _l.

Here, the functions a_j,l(t) fully encode the action of the interaction picture transformation on our basis operators.

Substituting Eq. (<ref>) in Eq. (<ref>), we obtain


    _(t) = ∑_j ṽ_j(t) _j,


where we use tildes to denote control fields in the interaction picture, and we have


    ṽ_j(t) = ∑_l a_l,j(t) v_l(t).


Proceeding analogously for (t) and using the series representation defined in Eq. (<ref>), we get


    _^(n)(t) = ∑_j w̃_j^(n)(t) _j,


with


    w̃_j^(n)(t) = ∑_l a_l,j(t) w_l^(n)(t).


We now return to the fundamental equations of our approach, Eqs. (<ref>), which need to be satisfied to cancel the effects
of (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^(n)(t) = 0.  Note that in general it is possible to have h_j (t) ≠ 0 while one must work
with the condition  w_j^(n)(t) = 0.  Moreover, even if w_j^(n)(t) can be controlled, it might have restrictions, e.g.
w_j^(n) must be time independent or has bandwidth limitations. In the following we show how to derive equations for w_j^(n)
(t) that obey Eqs. (<ref>) 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_j^(1) (t), which must obey Eq. (<ref>), since the
procedure for w_j^(n) (t) is similar.

We proceed by substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), 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:


    ∫_0^t_tw̃_j^(1)(t) = -∫_0^t_tṽ_j(t).


We stress that Eq. (<ref>) may be ill-defined since it is possible to have w̃_j^(1)(t) = 0 while
ṽ_j(t) ≠ 0 for certain values of j. We show in Sec. <ref> 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_j^(1)(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_j^(1)(t). Here, since w_j^(1)(t) must only have support on the interval [0, t_], we use a finite Fourier series
decomposition, 


    w_j^(1)(t) = ∑_k=0^k_max c_jk^(1)cos(ω_k t) + d_jk^(1)sin(ω_k t),


with ω_k = 2 π k/t_ and d_j0^(1) = 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. (<ref>). We stress that at this stage finding the first order correction that fulfills Eq. (<ref>) has
been reduced to determining a set of N_coeffs = N_op× (2 k_max + 1) coefficients.  Note that one could use
other basis functions for the decomposition, e.g., Slepian functions <cit.>. 

The sum in Eq. (<ref>) runs from 0 to k_max which allows us to limit the bandwidth of the field associated to
_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. (<ref>) to zero with the
exception of c_j0. If one requires w_j^(1)(0) = w_j^(1)(t_) =0, then one finds using Eq. (<ref>) that the
coefficients c_jk^(1) must obey ∑_k=0^k_max c_jk^(1) =0. For simplicity the summation in
Eq. (<ref>) 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. (<ref>) and (<ref>) in
the system of equations defined by Eq. (<ref>).  Since we know the explicit time dependence of
w̃_j^(1)(t), we can perform the time integration.  This leads to a system of time-independent N_op linear
equations than can be written in matrix form:


    M  ^(1) = ^(1).


Here, ^(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 ^(1) are known quantities: ^(1) parameterizes the error
Hamiltonian (t), whereas M encodes the dynamics of the ideal evolution generated by _0(t).

To be more explicit, the ^(1) is a vector of length N_op whose components are the spurious error-Hamiltonian
elements we wish to average out, 


    y^(1)_j = - ∫_0^t_tṽ_j(t).


^(1) is the vector of the N_coeffs unknown Fourier coefficients c_lk^(1) and d_lk^(1) that determine our
control corrections, c.f. Eq. (<ref>).  We order these as follows


    x^(1)_j = 
    	 
    		c_lk^(1)   if  j ≤ j_0, 
    
    		d_lk^(1)   if  j > j_0,


where j_0 = N_op(k_max + 1), and the indices l and k in Eq. (<ref>) are functions of j. We have 


    l = 
    	 
    		(j-1)//(k_max + 1) + 1    if  j ≤ j_0, 
    
    		(j - j_0 - 1)//k_max + 1    if  j > j_0,


and


    k = 
    	 
    		(j - 1) % (k_max + 1)    if  j ≤ j_0, 
    
    		(j - j_0 - 1) % k_max  + 1    if  j > j_0.


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. (<ref>)]: 


    M_ij = 
    	∫_0^t_t a_l,i(t) cos(ω_k t)    if  j ≤ j_0, 
    ∫_0^t_t a_l,i(t) sin(ω_k t)    if  j > j_0,


where l and k are given by Eqs. (<ref>) and (<ref>), respectively.  We stress that
Eqs. (<ref>) to (<ref>) are valid when the summation in Eq. (<ref>) 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. (<ref>) [see Appendix <ref>]. 

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
[There are also situations where the linear system has a solution, however the solution coefficients are large in comparison with
the errors that they should correct. In such cases, the perturbative series of the correction Hamiltonian typically does not converge. We
treat this problem in Sec. <ref>.], 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. <ref>.





§ 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.



 §.§ Strong Driving of a Two-Level System
 


As a first example we consider the problem of a two-level 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 <cit.>. An ad
hoc method based on time optimal control of a two-level system <cit.> was also proposed: it consists in
realizing Bang-Bang control with imperfect square control fields <cit.>. 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 <cit.>, but this transforms a low dimensional control
problem into a high dimensional one.



The Hamiltonian of a driven two-level system is given by


    _qubit (t) = ω_/2z + f_(t) cos(ω_ t) x,


where ω_ is the qubit splitting frequency, ω_ is the driving frequency, f_(t) is the driving envelope, and
we introduce the Pauli operators:


    x   = |0⟩⟨1|+|1⟩⟨0|,
    y   = i |0⟩⟨1| - i |1⟩⟨0|,
    z   = |1⟩⟨1|-|0⟩⟨0|.


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_(t) ≪ω_ ∀ t, Eq. (<ref>) allows one to generate rotations
around the x-axis if one sets ω_ = ω_. This is best understood in the frame rotating at the drive frequency
[We have used the unitary transformation _qubit (t) →_ (t) = _^† (t) _qubit_
(t) - i _^† (t) _t _ (t) with _ (t) = exp[-i ω_ t z/2].]. In this frame the Hamiltonian is
given by _ (t) = _,0 (t) + _(t) with 


    _,0 (t) = f_(t)/2x,


and


    _ (t) = v_,x(t) x + v_,y(t) y.


The coefficients v_,j(t) are given by


    v_,x (t)    = f_(t)/2cos( 2 ω_ t ), 
    
    		 v_,y (t)    = -f_(t)/2sin( 2 ω_ t ).


Here, the driving is set on resonance with the qubit frequency, i.e., ω_ = ω_. If the system is in the weak
driving limit, the fast oscillating terms (also known as counter-rotating terms) in _(t) can be neglected as they average
themselves out over the long evolution time set by the slow varying envelope function f_(t). As a consequence, one can
approximate _ (t) by _,0 (t). This is known as the rotating wave approximation (RWA). The resulting Hamiltonian
generates a rotation of angle θ(t_) around the x-axis, where we have introduced


    θ(t) = ∫_0^tt_1 f_(t_1).


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 _ (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. (<ref>) <cit.> [see Fig. <ref> (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. <ref>, we can mitigate the effects of _ (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. (<ref>)], we look for a correction of the form 


    _qubit (t) =   ∑_n [ g_x^(n)(t) cos(ω_ t) + g_y^(n)(t) sin(ω_ t)
    		] x.


Here, g_x^(n)(t) and g_y^(n)(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 ω_ = ω_ 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 Δ = ω_ - ω_.
Therefore, we consider the following modified Hamiltonian in the rotating frame 


    _,mod (t) = _,0 (t)+ _ (t) + ∑_n _^(n)(t).
 

In terms of the Pauli operators, _^(n)(t) is given by


    _^(n)(t) = w_,x^(n)(t) σ̂_x + w_,y^(n)(t) σ̂_y + w_,z^(n)(t)
    	σ̂_z,


with


    w_,x^(n)(t)    = g^(n)(t) cos(ω_ t),
    
    		w_,y^(n)(t)    = -g^(n)(t) sin(ω_ t),
    
    		w_,z^(n)(t)    = Δ^(n).


In practice, having a control field with two-quadratures driving [see Eq. (<ref>)] 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 <ref>]. Note that the modified
detuning is given by Δ=∑_n Δ^(n) in complete analogy to having the control fields represented by a series [see
Eq. (<ref>)].

Following the general strategy presented in Sec. <ref>, we first move to the interaction picture with respect to
_,0(t) [see Eq. (<ref>)]. In the interaction picture, _(t) [see Eq. (<ref>)] and the
control Hamiltonian _^(n)(t) [see Eq. (<ref>)] are respectively given by


    _, (t) = ṽ_,x(t) x + ṽ_,y(t) y + ṽ_,z(t) z,


with


    ṽ_,x(t)    = f_(t)/2cos( 2 ω_ t ),
    ṽ_,y(t)    = -f_(t)/2sin( 2 ω_ t ) cosθ,
    ṽ_,z(t)    = f_(t)/2sin( 2 ω_ t ) sinθ,


and


    _,^(n)(t) = w̃_,x^(n)(t) x + w̃_,y^(n)(t) y + w̃_,z^(n)(t) z


with 


    w̃_,x^(n)(t)    = g^(n)(t) cos(ω_ t),
    w̃_,y^(n)(t)    = -g^(n)(t) sin(ω_ t)cosθ  + Δ^(n)sinθ,
    w̃_,z^(n)(t)    = g^(n)(t) sin(ω_ t) sinθ  + Δ^(n)cosθ.


In Eqs. (<ref>) and (<ref>), we have omitted the explicit time dependence of θ for simplicity, i.e.,
θ = θ(t) [see Eq. (<ref>)]. The next step consists in expanding the control fields
w̃_,j^(n)(t) (j∈{x,y,z}) [see Eq.(<ref>)] into a Fourier series. However, before proceeding it is
useful to notice the special form of the functions w̃_,j^(n)(t): an unknown function that multiplies a known fast
oscillating function. It is therefore more suitable to just expand the unknown functions g_x^(n)(t), g_y^(n)(t) [see
Eq. (<ref>)], 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. <ref> and imposing the necessary constraints on the Fourier series.

If we constrain g_α=x,y^(n)(t) to be zero at t=0 and t=t_, which is often the case experimentally, we obtain the
following Fourier expansions


    g_α^(n)(t) = ∑_k=1^∞ c_α, k^(n)[ 1 - cos(ω_k t) ] + d_α, k^(n)sin(ω_k t),


and


    Δ^(n) = c_z, 0^(n) + ∑_k=1^∞ c_z, k^(n)cos(ω_k t) + d_z, k^(n)sin(ω_k t),


where ω_k = 2 π/t_. Since we have a total of three equations of the form of Eq. (<ref>) to solve (one
for each Pauli operator), we need at least three free parameters. Consequently, we can set all coefficients to zero in
Eqs. (<ref>) and (<ref>) except c_x, 1^(n), c_y, 1^(n) and c_z, 0^(n) [We have chosen
this set of coefficients for simplicity. In principle one could choose another set of three coefficients.]. With this choice,
Eqs. (<ref>) and (<ref>) reduce to 


    g_α=x,y^(n)(t) = c_α, 1^(n)[ 1 - cos(ω_1 t) ],


and


    Δ^(n) = c_z, 0^(n).


The final step is to find the value of the free parameters c_x, 1^(n), c_y, 1^(n) and Δ^(n). We start by
verifying that by substituting Eqs. (<ref>) and (<ref>) in Eq. (<ref>), we obtain a correction
Hamiltonian _,^(n)(t) [see Eq. (<ref>)] that depends, as desired, linearly on the free parameters c_x,
1^(n), c_y, 1^(n) and Δ^(n). The system of equations defining the first order coefficients [n=1, see
Eq. (<ref>)], is given by


    P__^(1) = _^(1),


where _^(1) = { c_x, 1^(1), c_y, 1^(1), Δ^(1)}^T is the vector of unknown coefficients [see
Eq. (<ref>)], _^(1)= - ∫_0^t_t{ṽ_,x (t), ṽ_,y (t), ṽ_,z
(t)}^T is the vector of the spurious error-Hamiltonian elements with ṽ_,j (t) (j∈{x,y,z}) defined in
Eq. (<ref>), and P_ [The matrices P_ and M [see Eq. (<ref>)], although they fulfill the same
purpose, have different matrix elements. The difference arises because we are expanding in a Fourier series the unknown envelope
functions g_α^(n) (t) (α=x,y) and the detuning Δ^(n) [see Eqs. (<ref>) and (<ref>)]
instead of the functions w̃_,j (t) (j∈{x,y,z})) [see Eq. (<ref>)].] is the matrix that
characterizes the evolution under the ideal Hamiltonian _,0 (t) [see Eq. (<ref>)].  The explicit matrix
elements of P_ can be found in Appendix <ref>. Higher-order correction Hamiltonians can be found in a similar
way.

In Fig. <ref> (a), we plot the average fidelity error ε <cit.> for a Hadamard gate generated
with an initial envelope


    f_(t) = θ_0/t_[1 - cos(2 π t/t_)],
	

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. <ref> (a), globally increases when ω_ t_ decreases, but around ω_ t_≃ 1
the error of the corrected evolution starts decreasing again. This can be understood by considering the limit t_→ 0
(ω_ t_→ 0). In this limit, we have  ṽ_,x (t) → f_ (t) /2 and ṽ_,y (t) =
ṽ_,z (t) → 0 [See Eq. (<ref>)] which implies that _ (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. <ref> (b), where we plot the coefficients of the correction versus the gate
time t_, the modified control sequences require control fields with diverging amplitudes. Restricting ourselves to gate times
close to unity (ω_ t_≃ 1), where the modified control sequences can be experimentally realized, our strategy
improves the error ε by more than two orders of magnitude. In Figs. <ref> (c) and (d), we compare the
original and corrected pulses for ω_ t_≈ 5. One can observe that the changes to the original pulse are small.
For convenience we write the nth order modified pulse as


    f_, mod (t) = f_,x^(n)(t) cos(ω_ t) + f_,y^(n)(t) sin(ω_ t),


where f_,x^(n)(t) = f_(t) + ∑_k=1^n g_x^(n)(t) and f_,y^(n)(t) = ∑_k=1^n g_y^(n)(t). When n
= 0 we have simply the original pulse, thus f_,x^(0)(t) = f_(t) and f_,y^(0)(t) = 0.





 §.§ Strong driving of a Parametrically Driven Cavity
 


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 <cit.> or to reach the single-photon strong coupling regime with optomechanical
systems using only linear resources <cit.>. Recently, optimal control techniques have been used to achieve
squeezing of an optomechanical oscillator at finite temperature <cit.>.

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. <ref>, 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. This can be achieved, e.g.,
in the microwave regime by modulating the magnetic flux through a SQUID loop (flux-pumped Josephson parametric
amplifier) <cit.>. We have


    _PDC (t) = ω_^† + f_(t) sin(ω_ t) ( + ^†)^2,
 

with  (^†) the bosonic annihilation (creation) operator. The frequency of the mode  is ω_ and the
drive has frequency ω_. 

It is convenient to introduce the operators <cit.> 


    μ̂_x    = 1/2( â^2 + â^† 2), 
    μ̂_y    = -i/2( â^2 - â^†2), 
    μ̂_z    = 1/2( â^†â + ââ^†),


which define (multiplied by the imaginary number -i) a Lie algebra with respect to the commutation operation [see Appendix
<ref>]. As mentioned earlier, since the Hamiltonian is quadratic, the three operators defined in
Eq. (<ref>) 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 by x̂ = ( + ^†)/√(2) and ŷ
= -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 ω_/2 = ω_, the Hamiltonian becomes _,R(t) = _,0(t) +
_(t) with 


    _,0(t) =  f_(t) μ̂_y,


and


    _(t) = f_(t) [ sin(2 ω_ t) μ̂_x - cos(2 ω_ t) μ̂_y + 2 sin(ω_
    	t) μ̂_z ].


In analogy with the qubit problem (see Sec. <ref>), one can neglect the fast oscillating Hamiltonian _(t) [see
Eq. (<ref>)] in the weak driving limit (RWA), i.e., when f_(t) ≪ω_ ∀ t. This results in
_,R(t) ≈_,0(t) and the generated dynamics corresponds to squeezing along the y-axis with a degree of
squeezing depending on r(t_), with 


    r(t) = ∫_0^tt_1 f_(t_1).


As one deviates from the weak driving limit, _(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. <ref> (d)].

To mitigate the effects of _(t) [see Eq. (<ref>)], we consider a control Hamiltonian that corresponds to just
changing the initial form of the parametric modulation. This leads to the correction Hamiltonian 


    _PDC (t) = ∑_n    [ g_x^(n)(t) cos(ω_ t) + g_y^(n)(t) sin(ω_ t) ] 
       ×(â + â^†)^2.


Furthermore, we are at liberty to drive the PDC at a frequency that is detuned from that of mode , 


    ω_/2 - ω_ = Δ,


with Δ = ∑_n Δ^(n) a static detuning.

Following the general procedure introduced in Sec. <ref> (see Appendix <ref>), we can
easily determine Δ^(n), g_x^(n)(t) and g_y^(n)(t). We stress that in this example we correct the unitary
evolution generated by Eq. (<ref>), 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. <ref> (a), we plot the degree of squeezing S as a function of the total evolution time t_ for the RWA
(red trace), the uncorrected (blue trace), and the corrected (green trace) evolutions. The degree of squeezing is given by 


    S = -10 log[( ⟨ŷ^2 ⟩_f - ⟨ŷ⟩_f^2 )/(⟨ŷ^2 ⟩_i - ⟨ŷ⟩_i^2 )]


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


    f_(t) =1/t_[1 - cos(2 π t/t_)].


Within the RWA the degree of squeezing is independent of the pulse width t_, since the squeezing depends just on
r(t_)=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_. In Fig. <ref> (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. <ref> (b), we plot the coefficients of the correction
Hamiltonian as a function of the total evolution time t_. 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. 



 §.§ Transmon Qubit



As a next example we consider the problem of realizing single-qubit gates with a transmon qubit <cit.>, where the logical
qubit states are encoded in the two lowest energy states of an anharmonic oscillator with eigenstates |n⟩ [see
Fig. <ref> (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) <cit.>. 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 <cit.> and theoretically <cit.> 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


    Ĥ_TLS(t) =   ω_/2z + (3 ω_/2  + α) |2⟩⟨2|
        + f_(t) cos(ω_ t) (x + ηx,1)


as an approximation of the weakly anharmonic oscillator. Here, ω_ is the frequency splitting between the energy levels
|0⟩ and |1⟩ while the frequency splitting between |1⟩ and |2⟩ is given by ω_ + α, where
α is the anharmonicity. We have also defined the operators


    x,12   = |1⟩⟨2| + |2⟩⟨1|,  y,12 = i|2⟩⟨1| - i|1⟩⟨2|,  
    x,02   = |0⟩⟨2| + |2⟩⟨0|,  y,02 = i|2⟩⟨0| - i|0⟩⟨2|,


which describe transitions between the logical qubit states and the leakage state |2⟩. These operators together with the
Pauli operators [see Eq. (<ref>)] and the operator |2⟩⟨2| form the operator basis for this problem
[i.e. the operators _j in Eq. (<ref>)-(<ref>)]. 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 ω_ and an envelope function f_(t). As one can see from
Eq. (<ref>), 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 ω_, the Hamiltonian is given by _(t) = _, 0(t) + V̂_(t), where


    _, 0(t) = α|2⟩⟨2|  + f_(t)/2x,


and


    V̂_(t) = ηf_(t)/2x,1.


Here, we assume that the driving is on resonance with the |0⟩↔|1⟩ transition, i.e., ω_ =
ω_. The Hamiltonian _, 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 Hamiltonian
V̂_(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 ω_ in Eqs. (<ref>) and (<ref>) (RWA)[In contrast
to the examples treated in Secs. <ref> and <ref>, counter-rotating terms are not a main source of error since
there is a relatively large separation between the driving frequency and the anharmonicity, i.e., f_(t)/α≫ f_(t)/(2
ω_). As a result the error due to leakage out of the computational space is much larger than the error due to
counter-rotating terms. We stress that our framework would allow us to simultaneously deal with leakage and the counter-rotating
terms, but neglecting the latter  allows us to work with simpler expressions.].



Given the constraints of the problem [see Eq. (<ref>)], 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) as


    _TLS(t) =   ∑_n [ g_x^(n)(t) cos(ω_ t) + g_y^(n)(t) sin(ω_ t)
    		] 
       × (x  + ηx,1),


with g_x^(n)(t) and g_y^(n)(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. (<ref>), 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_x^(n)(t) and g_y^(n)(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. <cit.>, 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 <cit.> (see
Appendix <ref>).

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. <ref> (a), we plot the average fidelity error as a function of the gate time t_. 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) <cit.>. The results show that the 6th order Magnus correction reduces the average fidelity error by more than
four orders of magnitude for small α t_, greatly outperforming the DRAG correction. In Fig. <ref> (b)
we compare the original and modified pulses for |α| t_ = 5. 
For convenience we write the nth order modified pulse as


    f_, mod (t) = f_,x^(n)(t) cos(ω_ t) + f_,y^(n)(t) sin(ω_ t),


where f_,x^(n)(t) = f_(t) + ∑_k=1^n g_x^(k)(t) and f_,y^(n)(t) = ∑_k=1^n g_y^(k)(t). The case n
= 0 corresponds to the original pulse, i.e., f_,x^(0)(t) = f_(t) and f_,y^(0)(t) = 0.

A legitimate concern at this point is related to the possibility of realizing the pulses obtained with the Magnus 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 |α|.  



 §.§ 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 number-dependent
arbitrary phase gates) combined with cavity displacements <cit.>. 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 <cit.>. 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 <cit.>, 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. <ref>, we stressed that it is crucial for the Hamiltonian _ (t) describing the modification of the control
fields to have terms involving all of the basis operators _j appearing in the Magnus expansion of the unitary
evolution generated by the error-Hamiltonian _ (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. (<ref>) 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) = _ + _(t), with


    _ = 1/2(ω_ + χ^†) z + ω_^†,


and


    _(t) = [f_x(t) cos(ω_ t) + f_y(t) sin(ω_ t)] x.


The Pauli operators α act on the Hilbert space of the qubit and have been defined in Eq. (<ref>). 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
ω_ but with different phases.

In the interaction picture with respect to _, the Hamiltonian becomes


    _(t)    = 1/2∑_n  ( f_x(t) [ cos(δω_n t) x - sin(δω_n t) y]  
       -f_y(t)  [ sin(δω_n t) x + cos(δω_n t) y] ) |n⟩⟨n|,


where δω_n = ω_ + χ n - ω_, |n⟩ is a bosonic number state, and we have neglected fast
oscillating terms. If the drive is now chosen to fulfil ω_ = ω_ + χ n_0, so that the drive is resonant for
a particular number-selected qubit transition, the Hamiltonian defined in Eq. (<ref>) can be written as _(t) =
_,0(t) + _(t).  Here


    _,0(t) = 1/2[ f_x(t) x - f_y(t) y] |n_0⟩⟨n_0|,


is the resonant part of the Hamiltonian defined in Eq. (<ref>) and allows one to generate a unitary operation in the
subspace spanned by {|g, n_0⟩, |e, n_0⟩}.  In contrast 


    _(t) =   1/2∑_n ≠ n_0( f_x(t) [ cos(δω_n t) x - sin(δω_n t)
    		y]  
         -f_y(t) [ sin(δω_n t) x + cos(δω_n t) y] )
    		|n⟩⟨n|


is the non-resonant part of Eq. (<ref>). 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 of _(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
_,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_ the state of the qubit is back to |g⟩. This leads to the accumulation of a Berry phase γ at
t=t_ for the qubit which conditioned on the state of the cavity being |n_0⟩. In other words, 



    _, 0 (t_) |g, n⟩ = 
    	
    		e^i γ|g, n⟩   if  n = n_0, 
    |g, n⟩   if   n ≠ n_0,
   

where _, 0 (t_) = exp[ -i ∫_0^t_t _,0 (t) ] is the unitary evolution
generated by the ideal Hamiltonian in Eq. (<ref>).  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 Hamiltonian [see
Eq. (<ref>)] by


    _(t) = ∑_n=0^N-1 [f_x,n(t) cos(ω_,n  t) + f_y,n(t) sin(ω_,n t)] x,
 

where ω_,n = ω_ + χ 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 off-resonant error interaction given by the generalization of  _(t)
[c.f. Eq. (<ref>)] 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 <cit.>, 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:


    _SNAP(t) =∑_n=0^N-1 [g_x,n(t) cos(ω_,n t) + g_y,n(t)
    		sin(ω_,n t)] x


where ω_,n = ω_ + χ n. Moving to the interaction picture with respect to _ [see
Eq. (<ref>)] and neglecting non-resonant terms, we obtain


    _(t) = 1/2∑_n=0^N-1 [g_x,n(t) x - g_y,n(t) y] |n⟩⟨n|.


In the interaction picture defined by _,0(t) [see Eq. (<ref>)], we find that the form of the
non-resonant error Hamiltonian is unchanged: 


    _,(t) = _(t),


since _,0(t) commutes with _(t); _,0(t) and _(t) act on orthogonal subspaces. On the other
hand, _(t) acts on the whole Hilbert space, and is transformed when moving to the interaction picture.  We find:


    _,   (t) = _,0^† (t) 1/2 [g_x,n_0(t) x - g_y,n_0(t) y]
    		|n_0⟩⟨n_0|_,0 (t) 
       + 1/2∑_n=0^N-1 (1 - δ_n, n_0) [g_x,n(t) x - g_y,n(t) y] |n⟩⟨n|.


The first term of Eq. (<ref>) 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. (<ref>) cannot correct
errors proportional to z|n⟩⟨n| (in the interaction picture) and which appear at 2nd order in the Magnus expansion of
_,(t) [see Eq. (<ref>)]. Unfortunately, an analysis of the Magnus expansion generated by Eq. (<ref>)
shows that these terms are by far the dominant source of errors which corrupt the ideal dynamics.  We are left with no choice but
to modify the general strategy of Sec. <ref> 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. (<ref>). 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 _,(t)
is missing important terms, these can nonetheless be dynamically generated. In the same way that _,(t) generates
problematic terms proportional to z|n⟩⟨n| at second order in the Magnus expansion, so can _,(t). Thus,
we look for a correction Hamiltonian ^(1)(t) that cancels the sum of the first two terms of the Magnus expansion:


    _1^(1)(t_) + _2^(1)(t_) = 0.


This is in contrast of the general strategy in Sec. <ref>, where we would just cancel the first term in the Magnus
expansion.  

We can use Eqs. (<ref>) and (<ref>) to write Eq. (<ref>) in terms of integrals involving
_,(t) [see Eq. (<ref>)] and _,^(1)(t) [see Eq. (<ref>) with n=1]. The explicit
equation can be found in the Appendix <ref>. Here, to be able to correct for the term proportional to z|n⟩⟨n|, we cannot discard the higher order term generated by _,^(1)(t) in _2^(1)(t_) (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 [The
number of free parameters has a lower limit corresponding to the number of equations, but it is typically useful to have
more free parameters than equations. In such a case, one can use Lagrange multipliers to find solutions that minimize the sum of
the modulus squared of the free parameters.]. Following the strategy presented in Sec. <ref>, 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. (<ref>) is quadratic in _,^(1)(t). This system of equations can be
solved numerically; see Appendix <ref>, 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. <ref>. In this case, since one is driving all
frequencies resonantly, the ideal unitary _,0(t) acts on the whole truncated Hilbert space of the cavity. As a
consequence, transforming the correction Hamiltonian _ (t) [see Eq. (<ref>)] to the interaction picture will
generate terms proportional to z|n⟩⟨n| for all values of n (see Appendix <ref>). 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 <cit.>, e.g., |0⟩_ = (|0⟩ + |4⟩)/√(2) and |1⟩_ = |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. (<ref>)].

In Fig. <ref> (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⟩_ = (|0⟩ + |4⟩)/√(2) and |1⟩_ = |2⟩. The envelope functions for n = 0 and n=4 are given by


    f_x, n(t) = 
    	π/2 t_[1 - cos(4 πt/t_)],      t < t_/2, 
    
    	0,     t≥t_/2,


and 


    f_y, n(t) = 
    	
    	0,    t < t_/2, 
    π/2 t_[1 - cos(4 πt/t_)],      t≥t_/2.


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. (<ref>)]. 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_ the difference can reach
almost four orders of magnitude. 

In Fig. <ref> (b) and (c) we show the spectrum of the original and modified pulses for a gate time of χ t_ = 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.



§ 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. 



§ ACKNOWLEDGEMENTS


AC acknowledges partial support from the Center for Novel Pathways to Quantum Coherence in Materials, an Energy Frontier Research
Center funded by the Department of Energy, Office of Science, Basic Energy Sciences.

61
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Fuchs et al.(2009)Fuchs,
  Dobrovitski, Toyli, Heremans, and Awschalom]fuchs2009
  author author G. D. Fuchs, author V. V. Dobrovitski, author D. M. Toyli, author F. J. Heremans, and author D. D. Awschalom, title title Gigahertz dynamics of a
  strongly driven single quantum spin, https://doi.org/10.1126/science.1181193 journal journal Science volume 326, pages
  1520 (year 2009)NoStop
[Khaneja et al.(2005)Khaneja, Reiss, Kehlet, Schulte-Herbrüggen, and Glaser]khaneja2005
  author author N. Khaneja, author T. Reiss,
  author C. Kehlet, author T. Schulte-Herbrüggen, and  author S. J. Glaser, title title Optimal control of coupled spin dynamics: design
  of nmr pulse sequences by gradient ascent algorithms, https://doi.org/http://dx.doi.org/10.1016/j.jmr.2004.11.004 journal journal Journal of Magnetic Resonance volume 172, pages 296  (year
  2005)NoStop
[Krotov and Feldman(1983)]krotov1983
  author author V. F. Krotov and author I. N. Feldman, title title An iterative method for
  solving problems of optimal control, @noop  journal
  journal Engineering Cybernetics volume
  21, pages 123  (year 1983)NoStop
[Somlói et al.(1993)Somlói, Kazakov, and Tannor]somloi1993
  author author J. Somlói, author V. A. Kazakov, and author D. J. Tannor, title title Controlled dissociation of
  I_2 via optical transitions between the X and B electronic
  states, https://doi.org/https://doi.org/10.1016/0301-0104(93)80108-L journal journal Chemical Physics volume 172, pages 85  (year
  1993)NoStop
[Doria et al.(2011)Doria,
  Calarco, and Montangero]doria2011
  author author P. Doria, author T. Calarco,  and author S. Montangero,  title title Optimal control technique for
  many-body quantum dynamics, https://doi.org/10.1103/PhysRevLett.106.190501 journal
  journal Phys. Rev. Lett. volume 106,  pages 190501 (year 2011)NoStop
[Glaser et al.(2015)Glaser,
  Boscain, Calarco, Koch,
  Köckenberger, Kosloff, Kuprov, Luy, Schirmer, Schulte-Herbrüggen, Sugny, and Wilhelm]glaser2015
  author author S. J. Glaser, author U. Boscain,
  author T. Calarco, author C. P. Koch, author
  W. Köckenberger, author
  R. Kosloff, author I. Kuprov, author B. Luy, author S. Schirmer, author T. Schulte-Herbrüggen, author D. Sugny, and author F. K. Wilhelm, title title Training
  schrödinger's cat: quantum optimal control, https://doi.org/10.1140/epjd/e2015-60464-1 journal journal The European Physical Journal D volume
  69, pages 279 (year 2015)NoStop
[Machnes et al.(2018)Machnes, Assémat, Tannor, and  Wilhelm]machnes2018
  author author S. Machnes, author E. Assémat,
  author D. Tannor, and author F. K. Wilhelm, title title Tunable, flexible, and efficient optimization of
  control pulses for practical qubits, https://doi.org/10.1103/PhysRevLett.120.150401 journal
  journal Phys. Rev. Lett. volume 120,  pages 150401 (year 2018)NoStop
[Werschnik and Gross(2007)]werschnik2007
  author author J. Werschnik and author E. K. U. Gross, title title Quantum optimal control
  theory, https://doi.org/10.1088/0953-4075/40/18/r01 journal journal Journal of Physics B: Atomic, Molecular and
  Optical Physics volume 40, pages R175
  (year 2007)NoStop
[Kirkpatrick et al.(1983)Kirkpatrick, Gelatt, and Vecchi]kirkpatrick1983
  author author S. Kirkpatrick, author C. D. Gelatt, and author M. P. Vecchi, title title Optimization by simulated
  annealing, https://doi.org/10.1126/science.220.4598.671 journal journal Science volume
  220, pages 671 (year 1983)NoStop
[Swendsen and Wang(1986)]swendsen1986
  author author R. H. Swendsen and author J.-S. Wang, title title Replica monte carlo
  simulation of spin-glasses, https://doi.org/10.1103/PhysRevLett.57.2607 journal journal Phys. Rev. Lett. volume 57, pages 2607 (year 1986)NoStop
[Whitley(1994)]whitley1994
  author author D. Whitley, title title A genetic algorithm
  tutorial, https://doi.org/10.1007/BF00175354 journal
  journal Statistics and Computing volume
  4, pages 65 (year 1994)NoStop
[Wenzel and Hamacher(1999)]wenzel1999
  author author W. Wenzel and author K. Hamacher, title title Stochastic tunneling
  approach for global minimization of complex potential energy landscapes,  https://doi.org/10.1103/PhysRevLett.82.3003 journal
  journal Phys. Rev. Lett. volume 82,  pages 3003 (year 1999)NoStop
[Motzoi et al.(2009)Motzoi,
  Gambetta, Rebentrost, and Wilhelm]motzoi2009
  author author F. Motzoi, author J. M. Gambetta, author P. Rebentrost, and author F. K. Wilhelm, title title Simple pulses for
  elimination of leakage in weakly nonlinear qubits, https://doi.org/10.1103/PhysRevLett.103.110501 journal
  journal Phys. Rev. Lett. volume 103,  pages 110501 (year 2009)NoStop
[Economou and Barnes(2015)]economou2015
  author author S. E. Economou and author E. Barnes, title title Analytical approach to
  swift nonleaky entangling gates in superconducting qubits, https://doi.org/10.1103/PhysRevB.91.161405 journal journal Phys. Rev. B volume 91, pages 161405 (year 2015)NoStop
[Demirplak and Rice(2003)]demirplak2003
  author author M. Demirplak and author S. A. Rice, title title Adiabatic Population
  Transfer with Control Fields, http://dx.doi.org/10.1021/jp030708a
  journal journal J. Phys. Chem. A volume 107, pages 9937 (year
  2003)NoStop
[Demirplak and Rice(2008)]demirplak2008
  author author M. Demirplak and author S. A. Rice, title title On the consistency,
  extremal, and global properties of counterdiabatic fields, https://doi.org/http://dx.doi.org/10.1063/1.2992152 journal
  journal The Journal of Chemical Physics volume 129, pages 154111 (year
  2008)NoStop
[Berry(2009)]berry2009
  author author M. V. Berry, title title Transitionless quantum
  driving, http://stacks.iop.org/1751-8121/42/i=36/a=365303
  journal journal Journal of Physics A:
  Mathematical and Theoretical volume 42, pages 365303 (year 2009)NoStop
[Ibáñez et al.(2012)Ibáñez, Chen, Torrontegui,
  Muga, and Ruschhaupt]ibanez2012
  author author S. Ibáñez, author X. Chen,
  author E. Torrontegui, author J. G. Muga, and author A. Ruschhaupt, title
  title Multiple schrödinger pictures and dynamics in shortcuts
  to adiabaticity, https://doi.org/10.1103/PhysRevLett.109.100403
  journal journal Phys. Rev. Lett. volume 109, pages 100403 (year
  2012)NoStop
[Chen and Muga(2012)]chen2012
  author author X. Chen and author J. G. Muga, title title Engineering of fast
  population transfer in three-level systems, https://doi.org/10.1103/PhysRevA.86.033405 journal journal Phys. Rev. A volume 86, pages 033405 (year 2012)NoStop
[Baksic et al.(2016)Baksic,
  Ribeiro, and Clerk]baksic2016
  author author A. Baksic, author H. Ribeiro,  and author A. A. Clerk,  title title Speeding up adiabatic quantum state
  transfer by using dressed states, http://link.aps.org/doi/10.1103/PhysRevLett.116.230503 journal journal Phys. Rev. Lett. volume 116, pages 230503 (year
  2016)NoStop
[Ribeiro and Clerk(2019)]ribeiro2019
  author author H. Ribeiro and author A. A. Clerk, title title Accelerated adiabatic
  quantum gates: Optimizing speed versus robustness, https://doi.org/10.1103/PhysRevA.100.032323 journal journal Phys. Rev. A volume 100, pages 032323 (year 2019)NoStop
[Ribeiro et al.(2017)Ribeiro, Baksic, and Clerk]ribeiro2017
  author author H. Ribeiro, author A. Baksic,  and author A. A. Clerk,  title title Systematic magnus-based approach for
  suppressing leakage and nonadiabatic errors in quantum dynamics, https://doi.org/10.1103/PhysRevX.7.011021 journal journal Phys. Rev. X volume 7, pages
  011021 (year 2017)NoStop
[Magnus(1954)]magnus1954
  author author W. Magnus, title title On the exponential
  solution of differential equations for a linear operator, https://doi.org/10.1002/cpa.3160070404 journal journal Communications on Pure and Applied Mathematics volume 7, pages 649 (year
  1954)NoStop
[Blanes et al.(2009)Blanes,
  Casas, Oteo, and Ros]blanes2009
  author author S. Blanes, author F. Casas,
  author J. A. Oteo, and  author J. Ros, title title The magnus expansion and some of its
  applications, https://doi.org/http://dx.doi.org/10.1016/j.physrep.2008.11.001 journal journal Physics Reports volume
  470, pages 151  (year 2009)NoStop
[Heeres et al.(2015)Heeres,
  Vlastakis, Holland, Krastanov, Albert, Frunzio, Jiang, and Schoelkopf]heeres2015
  author author R. W. Heeres, author B. Vlastakis,
  author E. Holland, author S. Krastanov, author
  V. V. Albert, author
  L. Frunzio, author L. Jiang, and author R. J.  Schoelkopf, title title Cavity
  state manipulation using photon-number selective phase gates, https://doi.org/10.1103/PhysRevLett.115.137002 journal
  journal Phys. Rev. Lett. volume 115,  pages 137002 (year 2015)NoStop
[Krastanov et al.(2015)Krastanov, Albert, Shen, Zou, Heeres, Vlastakis, Schoelkopf, and Jiang]krastanov2015
  author author S. Krastanov, author V. V. Albert, author C. Shen,
  author C.-L. Zou, author R. W. Heeres, author
  B. Vlastakis, author
  R. J. Schoelkopf, and author
  L. Jiang, title title Universal control of an oscillator with dispersive coupling to a
  qubit, https://doi.org/10.1103/PhysRevA.92.040303 journal journal Phys. Rev. A volume
  92, pages 040303 (year 2015)NoStop
[Barends et al.(2014)Barends, Kelly, Megrant, Veitia, Sank, Jeffrey, White, Mutus, Fowler, Campbell, Chen, Chen, Chiaro, Dunsworth, Neill, O'Malley, Roushan, Vainsencher,
  Wenner, Korotkov, Cleland,  and Martinis]barends2014
  author author R. Barends, author J. Kelly,
  author A. Megrant, author A. Veitia, author
  D. Sank, author E. Jeffrey, author T. C.  White, author J. Mutus, author A. G. Fowler,
  author B. Campbell, author Y. Chen, author
  Z. Chen, author B. Chiaro, author A. Dunsworth, author C. Neill, author P. O'Malley, author P. Roushan, author A. Vainsencher, author J. Wenner, author A. N.  Korotkov, author A. N.  Cleland, and author J. M.  Martinis, title title
  Superconducting quantum circuits at the surface code threshold for fault
  tolerance, https://doi.org/10.1038/nature13171 journal journal Nature volume 508,  pages 500 (year 2014)NoStop
[Martinis and Geller(2014)]martinis2014
  author author J. M. Martinis and author M. R. Geller, title title Fast adiabatic qubit gates
  using only _z control, https://doi.org/10.1103/PhysRevA.90.022307 journal journal Phys. Rev. A volume 90, pages 022307 (year 2014)NoStop
[Sels and Polkovnikov(2017)]sels2017
  author author D. Sels and author A. Polkovnikov, title title Minimizing
  irreversible losses in quantum systems by local counterdiabatic driving,  journal journal Proceedings of the National
  Academy of Sciences https://doi.org/10.1073/pnas.1619826114
  10.1073/pnas.1619826114 (year 2017)NoStop
[Claeys et al.(2019)Claeys,
  Pandey, Sels, and Polkovnikov]claeys2019
  author author P. W. Claeys, author M. Pandey,
  author D. Sels, and author A. Polkovnikov, title title Floquet-engineering counterdiabatic protocols in
  quantum many-body systems, https://doi.org/10.1103/PhysRevLett.123.090602 journal
  journal Phys. Rev. Lett. volume 123,  pages 090602 (year 2019)NoStop
[Yurke et al.(1986)Yurke,
  McCall, and Klauder]yurke1986
  author author B. Yurke, author S. L. McCall,  and author J. R. Klauder,  title title SU(2) and SU(1,1)
  interferometers, https://doi.org/10.1103/PhysRevA.33.4033
  journal journal Phys. Rev. A volume 33, pages 4033 (year
  1986)NoStop
[Slepian and Pollak(1961)]slepian1961
  author author D. Slepian and author H. O. Pollak, title title Prolate spheroidal wave
  functions, fourier analysis and uncertainty — i, https://doi.org/10.1002/j.1538-7305.1961.tb03976.x journal
  journal The Bell System Technical Journal volume 40, pages 43 (year 1961)NoStop
[Note1()]Note1
  note There are also situations where the linear system has a
  solution, however the solution coefficients are large in comparison with the
  errors that they should correct. In such cases, the perturbative series of
  the correction Hamiltonian typically does not converge. We treat this problem
  in Sec. <ref>.Stop
[Scheuer et al.(2014)Scheuer, Kong, Said, Chen,
  Kurz, Marseglia, Du,
  Hemmer, Montangero, Calarco,
  Naydenov, and Jelezko]scheuer2014
  author author J. Scheuer, author X. Kong,
  author R. S. Said, author J. Chen, author
  A. Kurz, author L. Marseglia, author J. Du, author P. R.  Hemmer, author S. Montangero, author T. Calarco,
  author B. Naydenov, and  author F. Jelezko, title title Precise qubit control beyond the rotating wave
  approximation, https://doi.org/10.1088/1367-2630/16/9/093022
  journal journal New Journal of Physics volume 16, pages 093022 (year
  2014)NoStop
[Boscain and Mason(2006)]boscain2006
  author author U. Boscain and author P. Mason, title title Time minimal trajectories
  for a spin 1∕2 particle in a magnetic field, https://doi.org/10.1063/1.2203236 journal journal Journal of Mathematical Physics volume
  47, pages 062101 (year 2006)NoStop
[Garon et al.(2013)Garon,
  Glaser, and Sugny]garon2013
  author author A. Garon, author S. J. Glaser,  and author D. Sugny, title title Time-optimal control of su(2) quantum
  operations, https://doi.org/10.1103/PhysRevA.88.043422 journal journal Phys. Rev. A volume
  88, pages 043422 (year 2013)NoStop
[Hirose and Cappellaro(2018)]hirose2018
  author author M. Hirose and author P. Cappellaro, title title Time-optimal control
  with finite bandwidth, https://doi.org/10.1007/s11128-018-1845-6
  journal journal Quantum Information Processing  volume 17, pages 88 (year
  2018)NoStop
[Deng et al.(2015)Deng,
  Orgiazzi, Shen, Ashhab, and  Lupascu]deng2015
  author author C. Deng, author J.-L. Orgiazzi,
  author F. Shen, author
  S. Ashhab, and author
  A. Lupascu, title title Observation of floquet states in a strongly driven artificial
  atom, https://doi.org/10.1103/PhysRevLett.115.133601 journal journal Phys. Rev. Lett. volume 115, pages 133601 (year
  2015)NoStop
[Deng et al.(2016)Deng,
  Shen, Ashhab, and Lupascu]deng2016
  author author C. Deng, author F. Shen, author S. Ashhab, and author
  A. Lupascu, title title Dynamics of a two-level system under strong driving: Quantum-gate
  optimization based on floquet theory, https://doi.org/10.1103/PhysRevA.94.032323 journal journal Phys. Rev. A volume 94, pages 032323 (year 2016)NoStop
[Note2()]Note2
  note We have used the unitary transformation hatA5EH_qubit (t) →hatA5EH_R(t) = hatA5ES_d^† (t) hatA5EH_qubithatA5ES_d(t) - i hatA5ES_d^† (t) ∂ _t hatA5ES_d(t) with hatA5ES_d(t) =  oexp[-i
  ω _dt hatA5Eσ_z/2].Stop
[Bloch and Siegert(1940)]bloch1940
  author author F. Bloch and author A. Siegert, title title Magnetic resonance for
  nonrotating fields, https://doi.org/10.1103/PhysRev.57.522
  journal journal Phys. Rev. volume 57, pages 522 (year 1940)NoStop
[Note3()]Note3
  note We have chosen this set of coefficients for simplicity. In
  principle one could choose another set of three coefficients.Stop
[Note4()]Note4
  note The matrices P_q and M [see
  Eq. ], although they
  fulfill the same purpose, have different matrix elements. The difference
  arises because we are expanding in a Fourier series the unknown envelope
  functions g_α ^(n) (t) (α =x,y) and the detuning Δ
  ^(n) [see Eqs.  and
  ] instead of the
  functions tildeA7Ew_q,j (t)
  (j∈{x,y,z})) [see Eq. ].Stop
[Pedersen et al.(2007)Pedersen, Møller, and Mølmer]pedersen2007
  author author L. H. Pedersen, author N. M. Møller, and author K. Mølmer, title title Fidelity of quantum
  operations, https://doi.org/10.1016/j.physleta.2007.02.069
  journal journal Physics Letters A volume 367, pages 47  (year
  2007)NoStop
[Burd et al.(2019)Burd,
  Srinivas, Bollinger, Wilson,
  Wineland, Leibfried, Slichter, and Allcock]burd2019
  author author S. C. Burd, author R. Srinivas,
  author J. J. Bollinger, author A. C. Wilson, author
  D. J. Wineland, author
  D. Leibfried, author
  D. H. Slichter, and author
  D. T. C. Allcock, title
  title Quantum amplification of mechanical oscillator motion,  https://doi.org/10.1126/science.aaw2884 journal
  journal Science volume 364, pages 1163 (year 2019)NoStop
[Lemonde et al.(2016)Lemonde, Didier, and Clerk]lemonde2016
  author author M.-A. Lemonde, author N. Didier,  and author A. A. Clerk,  title title Enhanced nonlinear interactions in
  quantum optomechanics via mechanical amplification, https://doi.org/10.1038/ncomms11338 journal journal Nature Communications volume 7,  pages 11338 (year 2016)NoStop
[Basilewitsch et al.(2019)Basilewitsch, Koch, and Reich]basilewitsch2019
  author author D. Basilewitsch, author C. P. Koch, and author D. M. Reich, title title Quantum optimal control for
  mixed state squeezing in cavity optomechanics, https://doi.org/10.1002/qute.201800110 journal journal Advanced Quantum Technologies volume 2,  pages 1800110 (year 2019)NoStop
[Ojanen and Salo(2007)]ojanen2007
  author author T. Ojanen and author J. Salo,  title title Possible scheme for on-chip element
  for squeezed microwave generation, https://doi.org/10.1103/PhysRevB.75.184508 journal journal Phys. Rev. B volume 75, pages 184508 (year 2007)NoStop
[Yamamoto et al.(2008)Yamamoto, Inomata, Watanabe, Matsuba, Miyazaki, Oliver, Nakamura, and Tsai]yamamoto2008
  author author T. Yamamoto, author K. Inomata,
  author M. Watanabe, author K. Matsuba, author
  T. Miyazaki, author
  W. D. Oliver, author
  Y. Nakamura, and author
  J. S. Tsai, title title Flux-driven josephson parametric amplifier, https://doi.org/10.1063/1.2964182 journal journal Applied Physics Letters volume 93,  pages 042510 (year 2008)NoStop
[Koch et al.(2007)Koch,
  Yu, Gambetta, Houck,
  Schuster, Majer, Blais,
  Devoret, Girvin, and Schoelkopf]koch2007
  author author J. Koch, author T. M. Yu,
  author J. Gambetta, author A. A. Houck, author
  D. I. Schuster, author
  J. Majer, author A. Blais, author M. H.  Devoret, author S. M.  Girvin, and author R. J.  Schoelkopf, title title
  Charge-insensitive qubit design derived from the cooper pair box, https://doi.org/10.1103/PhysRevA.76.042319 journal journal Phys. Rev. A volume 76, pages 042319 (year 2007)NoStop
[Gambetta et al.(2011)Gambetta, Motzoi, Merkel, and  Wilhelm]gambetta2011
  author author J. M. Gambetta, author F. Motzoi,
  author S. T. Merkel, and  author F. K. Wilhelm,  title title Analytic control methods for
  high-fidelity unitary operations in a weakly nonlinear oscillator, https://doi.org/10.1103/PhysRevA.83.012308 journal journal Phys. Rev. A volume 83, pages 012308 (year 2011)NoStop
[Chen et al.(2016)Chen,
  Kelly, Quintana, Barends,
  Campbell, Chen, Chiaro,
  Dunsworth, Fowler, Lucero,
  Jeffrey, Megrant, Mutus,
  Neeley, Neill, O'Malley,
  Roushan, Sank, Vainsencher,
  Wenner, White, Korotkov,  and Martinis]chen2016
  author author Z. Chen, author J. Kelly,
  author C. Quintana, author R. Barends, author
  B. Campbell, author
  Y. Chen, author B. Chiaro, author A. Dunsworth, author A. G.  Fowler, author E. Lucero, author E. Jeffrey,
  author A. Megrant, author J. Mutus, author
  M. Neeley, author C. Neill, author P. J. J.  O'Malley, author P. Roushan, author D. Sank,
  author A. Vainsencher, author J. Wenner, author
  T. C. White, author
  A. N. Korotkov, and author
  J. M. Martinis, title
  title Measuring and suppressing quantum state leakage in a
  superconducting qubit, https://doi.org/10.1103/PhysRevLett.116.020501 journal
  journal Phys. Rev. Lett. volume 116,  pages 020501 (year 2016)NoStop
[Note5()]Note5
  note In contrast to the examples treated in Secs. <ref>
  and <ref>, counter-rotating terms are not a main source of error
  since there is a relatively large separation between the driving frequency
  and the anharmonicity, i.e., f_T(t)/ |α
  | ≫ f_T(t)/(2 ω _d). As a
  result the error due to leakage out of the computational space is much larger
  than the error due to counter-rotating terms. We stress that our framework
  would allow us to simultaneously deal with leakage and the counter-rotating
  terms, but neglecting the latter allows us to work with simpler
  expressions.Stop
[Schutjens et al.(2013)Schutjens, Dagga, Egger, and  Wilhelm]schutjens2013
  author author R. Schutjens, author F. A. Dagga, author D. J. Egger,  and author F. K. Wilhelm,  title title Single-qubit gates in
  frequency-crowded transmon systems, https://doi.org/10.1103/PhysRevA.88.052330 journal journal Phys. Rev. A volume 88, pages 052330 (year 2013)NoStop
[Moore(1920)]moore1920
  author author E. H. Moore, title title On the reciprocal of the
  general algebraic matrix, https://projecteuclid.org:443/euclid.bams/1183425340 journal
  journal Bull. Amer. Math. Soc. volume
  26, pages 395 (year 1920)NoStop
[Bjerhammar(1951)]bjerhammar1951
  author author A. Bjerhammar, @noop  title Application of
  calculus of matrices to method of least squares; with special references to
  geodetic calculations, Vol. volume 49 (publisher Trans. Roy. Inst. Tech. Stockholm, year
  1951)NoStop
[Penrose(1955)]penrose1955
  author author R. Penrose, title title A generalized inverse for
  matrices, https://doi.org/10.1017/S0305004100030401 journal journal Mathematical Proceedings of the Cambridge
  Philosophical Society volume 51, pages
  406–413 (year 1955)NoStop
[Heeres et al.(2017)Heeres,
  Reinhold, Ofek, Frunzio,
  Jiang, Devoret, and Schoelkopf]heeres2017
  author author R. W. Heeres, author P. Reinhold,
  author N. Ofek, author
  L. Frunzio, author L. Jiang, author M. H.  Devoret, and author R. J.  Schoelkopf, title title
  Implementing a universal gate set on a logical qubit encoded in an
  oscillator, https://doi.org/10.1038/s41467-017-00045-1 journal journal Nature Communications volume 8, pages 94 (year 2017)NoStop
[Reinhold et al.(2019)Reinhold, Rosenblum, Ma,
  Frunzio, Jiang, and Schoelkopf]reinhold2019
  author author P. Reinhold, author S. Rosenblum, author W.-L. Ma, author L. Frunzio,
  author L. Jiang, and author R. J. Schoelkopf, title title Error-corrected gates on an encoded qubit,  @noop  journal journal arXiv e-prints  , eid arXiv:1907.12327 (year 2019), https://arxiv.org/abs/1907.12327 arXiv:1907.12327 [quant-ph] NoStop
[Note6()]Note6
  note The number of free parameters has a lower limit
  corresponding to the number of equations, but it is typically useful to have
  more free parameters than equations. In such a case, one can use Lagrange
  multipliers to find solutions that minimize the sum of the modulus squared of
  the free parameters.Stop
[Michael et al.(2016)Michael, Silveri, Brierley, Albert, Salmilehto, Jiang, and  Girvin]michael2016
  author author M. H. Michael, author M. Silveri,
  author R. T. Brierley, author V. V. Albert, author
  J. Salmilehto, author
  L. Jiang, and author
  S. M. Girvin, title
  title New class of quantum error-correcting codes for a bosonic
  mode, https://doi.org/10.1103/PhysRevX.6.031006 journal journal Phys. Rev. X volume
  6, pages 031006 (year 2016)NoStop
















empty



Appendices for: Engineering Fast High-Fidelity Quantum Operations With Constrained Interactions





§ NONRESONANT PERTURBATIONS



In Sec. <ref>, where we introduced the general strategy to correct unitary evolution, we have taken as a starting point
the Hamiltonian defined in Eq. (<ref>), where the spurious error-Hamiltonian (t) is multiplied by a small parameter
ϵ, such that (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,


    _(t) = 1/t_[ e^-i ω_ tq̂(t/t_) + e^i ω_ tq̂^†(t/t_) ],


where ω_ is the frequency at which the unwanted Hamiltonian _(t) oscillates and
q̂(t/t_) is a bounded operator. Considering the Magnus expansion of the evolution operator _(t) = exp[-i ∫_0^tt _(t)], where  is the time-ordering operator, we have to first order, at t = t_,


    _1(t_)    = ∫_0^t_t _(t)
    	
       = 1/ω_ t_[ i e^-i ω_ tq̂(t/t_) |_0^t_ + H.c. ] 
    	+ 1/ω_ t_∫_0^t_t [ i e^-i ω_ td/tq̂(t/t_) 
    	+ H.c. ],


where Eq. (<ref>) was obtained by integrating Eq. (<ref>) by parts. One can perform the integration by
parts repeatedly and obtain a series expansion for _1(t_) in powers of (ω_ t_)^-1. As one can see
from Eq. (<ref>), the leading order of _1(t_) is (ω_ t_)^-1. This suggests that
(ω_ t_)^-1 plays the role of the small parameter ϵ. Taking this analogy one step further, one would
naively assume that the leading order of _n(t_) is (ω_ t_)^-n, since this is the case when the
error-Hamiltonian is proportional to ϵ. This is, however, not true for oscillating Hamiltonians [see
Eq. (<ref>)]. To show that, let us calculate the second order term of the Magnus expansion,


    _2(t_) = -∫_0^t_t_1∫_0^t_1t_2  [_(t_1), _(t_2)].


Substituting Eq. (<ref>) in Eq. (<ref>) and performing the integration in t_2 by parts, we obtain


    _2(t_) = -1/t_^2∫_0^t_t_1 [ e^-i ω_ t_1q̂(t_1/t_) + H.c.,
    	-e^-i ω_ t_2/i ω_q̂(t_2/t_)|_0^t_1 + H.c. ] + 𝒪[(ω_
    	t_)^-2],


where we have omitted terms that are 𝒪[(ω_ t_)^-2] and higher. Further simplifying
Eq. (<ref>), we find [Note that the terms calculated at t_2 = 0 give rise to terms proportional to
(ω_ t_)^-2. Therefore we omit these terms in Eq. (<ref>).]


    _2(t_)    = -1/t_^2∫_0^t_t_1 {i/ω_[ q̂(t_1/t_),
    	q̂^†(t_1/t_) ] + H.c. } + 𝒪[(ω_ t_)^-2]
    	
       = -1/ω_ t_∫_0^1 x { i [ q̂(x),  q̂^†(x) ] + H.c. }
    	+ 𝒪[(ω_ t_)^-2],


which shows that _2(t_) scales to leading order like (ω_ t_)^-1. 

Nevertheless, since we know the Magnus expansion converges, we must have that the leading order of _n(t) scales with
higher powers of (ω_ t_)^-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^(0)(t_) is given by


    _j^(0)(t_) =
    	
    		0 + 𝒪[ (ω_ t_)^-j/2 ]    if  j  is even, 
    
    		0 + 𝒪[ (ω_ t_)^-(j+1)/2 ]    if  j  is odd.


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 by (t). 



§ THE MAGNUS EXPANSION



In the main text we showed the expression only for the first two terms of the Magnus expansion. For the bosonic system, however,
we have obtained a sixth order correction. The equations for the first four terms of the Magnus expansion are:


    ∂_t Ω̂_1 =    -i _,
    	
    ∂_t Ω̂_2 =    -1/2 [Ω̂_1, ∂_t Ω̂_1],
    	
    ∂_t Ω̂_3 =    -1/2 [Ω̂_2, ∂_t Ω̂_1] 
    	+ 1/12 [ Ω̂_1, [Ω̂_1, ∂_t Ω̂_1 ]],
    	
    ∂_t Ω̂_4 =    -1/2 [ Ω̂_3, ∂_t Ω̂_1] 
    	+ 1/12 [ Ω̂_2, [Ω̂_1, ∂_t Ω̂_1]]
    	+ 1/12 [ Ω̂_1, [Ω̂_2, ∂_t Ω̂_1]],


A generator for arbitrary order terms, which is convenient to obtain high order terms, can be found in the subsection 2.3 of
Ref. <cit.>.

When trying to calculate the Magnus expansion terms, one might be tempted to calculate the terms iteratively, i.e. firstly
integrate eq. (<ref>) to obtain Ω̂_1, then use this result to integrate Eq. (<ref>) and obtain
Ω̂_2, and so on. It is nonetheless much more efficient to treat all the terms that one intends to calculate as a
system of differential equations and solve them simultaneously. In this work we solved the differential equations using the
DifferentialEquations.jl package <cit.> from the Julia programming language <cit.>.	



§ ARBITRARY ORDER CORRECTIONS



A nth order correction, that generalizes Eq. (<ref>) of the main text, must satisfy the following
relation <cit.>:


    ∫_0^t_t _^(n)(t) = -i ∑_k=1^n_k^(n-1)(t_).


Since the set of operators {_j} forms a basis and {-i _j} generates a Lie algebra (see main text), we can write
_k^(l) as a linear combination of the operators {_j},


    Ω̂_k^(l)(t) = ∑_j Ω_k,j^(l)(t) Â_j.


Substituting Eqs. (<ref>) of the main text and (<ref>) in Eq. (<ref>), we obtain


    ∫_0^t_t w̃_j^(n)(t) = -i ∑_k=1^nΩ_k,j^(n-1)(t_) .


The next steps are very similar to what was done for the first order correction. First, we expand w_j^(n)(t) in a Fourier
series [see Eq. (<ref>) of the main text]. Since w̃_j^(n)(t) = ∑_l w_j^(n)(t) a_l,j (t), we can substitute
Eq. (<ref>) of the main text in Eq. (<ref>), and we obtain


    M  ^(n) = ^(n),


where M is the same known (N_op× N_coeffs) matrix obtained for n=1 [see Eq. (<ref>) of the main
text] and which encodes the dynamics of the ideal evolution generated _0 (t), ^(n) is the vector of the
N_coeffs unknown Fourier coefficients c_lk^(n) and d_lk^(n) [see Eq. (<ref>)], and ^(n) is
the known vector of spurious elements we wish to average out. In the case where the summation in Eq. (<ref>) runs from
0 to k_max for all values of j, the explicit expressions for the elements of the matrix M are given by
Eq. (<ref>). The elements of the vector ^(n) are


    ^(n)_j = 
    	 
    		c_l,k^(n)   if  j ≤ j_0, 
    
    		d_l,k^(n)   if  j > j_0.


Here j_0 = N_op(k_max + 1), and l and k are given by the Eqs. (<ref>) and (<ref>) of
the main text. The elements of ^(n) are given by


    ^(n)_j = -i  ∑_k=1^nΩ_k,j^(n-1)(t_).
	



§ STRONG DRIVING OF A TWO-LEVEL SYSTEM




 §.§ Derivative-Based Correction



While discussing the problem of strong driving of a two-level system, we mention that other choices for (t) were possible,
but that they all require more resources to be implemented experimentally. In this appendix, we illustrate this by considering the 
derivate-based control method introduced in Ref. <cit.>.

Let us first summarize the principle on which the derivate-based control method is built on. According to Eq. (<ref>) of the
main text, the first order correction term must satisfy


    ∫_0^t_t V̂_,(t) = - ∫_0^t_t Ŵ_I^(1)(t)


where V̂_,(t) is given by Eq. (<ref>) of the main text. Integrating the left hand side of Eq. (<ref>)
by parts, we find


    ∫_0^t_t V̂_,(t) =   f_ (t)/4 ω_[ sin( 2 ω_ t ) x + cos( 2 ω_ t ) 
    		cosθy - cos( 2 ω_ t ) sinθz] |_0^t_
       - 1/4 ω_∫_0^t_t {ḟ_ (t) sin( 2 ω_ t ) x 
    		+ [ḟ_ (t) cosθ - f_^2(t)sinθ] cos( 2 ω_ t)  y
                      + [ḟ_ (t) sinθ + f_^2(t)cosθ] cos( 2 ω_ t)  
    		z},


where we have omitted the explicit time dependence of θ for simplicity, i.e., θ = θ (t).

If the envelope function f_(t) vanishes at t = 0 and t = t_, Eq. (<ref>) becomes 


    ∫_0^t_t V̂_,(t) =    0
    		- 1/4 ω_∫_0^t_t {ḟ_ (t) sin( 2 ω_ t ) x 
    		+ [ḟ_ (t) cosθ - f_^2(t)sinθ] cos( 2 ω_ t)  y
                      + [ḟ_ (t) sinθ + f_^2(t)cosθ] cos( 2 ω_ t)  
    		z}.


Identifying Eq. (<ref>) with Eq. (<ref>), we find that a possible choice for _I^(1)(t) is 


    _^(1)(t) = w̃_1^(1)(t) x + w̃_2^(1)(t) y + w̃_3^(1)(t) z,


where


    w̃_1^(1)(t)    = 1/4 ω_ḟ_ (t) sin( 2 ω_ t),
    	
    w̃_2^(1)(t)    = 1/4 ω_[ḟ(t) cosθ - f^2(t)sinθ] cos( 2 ω_ t),
    	
    w̃_3^(1)(t)    = 1/4 ω_[ḟ(t) sinθ + f^2(t)cosθ] cos( 2 ω_ t).


The second order control Hamiltonian _I^(2)(t) can be found by considering the Magnus expansion associated to the unitary
evolution generated by _mod,I^(1) (t) = _ (t) + _^(1) (t). By construction, the first term of the Magnus
expansion is given by 


    Ω̂_1^(1) (t) = -f_ (t)/4 ω_[ sin( 2 ω_ t ) x + cos( 2 ω_ t ) 
    	cosθy - cos( 2 ω_ t ) sinθz],


and vanishes at both t=0 and t=t_ because f_ (0) = f_ (t_) = 0.

According to Eq. (<ref>) of the main text, the second term in the Magnus expansion is given by,


    Ω_2^(1)(t) - Ω_2^(1)(0) = -i/2∫_0^t t_1  [_ (t_1) + _^(1) (t_1),
    	Ω̂_1^(1) (t_1) ].


A natural choice for _I^(2)(t) is then (see Eq. (<ref>) of the main text)


    Ŵ_I^(2)(t) = -1/2[V̂_(t) + Ŵ_I^(1)(t), Ω̂_1(t) ],


Transforming _I^(1)(t) and _I^(2)(t) back to the original frame, we obtain the correction Hamiltonian (t) =
^(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 P_



In this appendix, we give explicitly the matrix elements of the matrix P_ introduced in Eq. (<ref>) of the
main text when discussing the strong driving of a qubit. We have


    P_,11   = ∫_0^t_t  [1 - cos(ω_1 t)] cos^2(ω_ t),
    	
    
    	P_,12   = ∫_0^t_t  [1 - cos(ω_1 t)] sin(ω_ t) cos(ω_ t),
    	
    
    	P_,13   = 0,
    	
    
    	P_,21    = -∫_0^t_t  [1 - cos(ω_1 t)] sin(ω_ t) cos(ω_ t) cosθ,
    	
    
    	P_,22    = -∫_0^t_t  [1 - cos(ω_1 t)] sin^2(ω_ t) cosθ,
    	
    
    	P_,23    = ∫_0^t_t sinθ,
    	
    
    	P_,31    = ∫_0^t_t  [1 - cos(ω_1 t)] sin(ω_ t) cos(ω_ t) sinθ,
    	
    
    	P_,32    = ∫_0^t_t  [1 - cos(ω_1 t)] sin^2(ω_ t) sinθ,
    	
    
    	P_,33    = ∫_0^t_t cosθ,


where once more we have omitted the explicit time dependence of θ.



§ STRONG DRIVING OF A PARAMETRICALLY DRIVEN CAVITY




 §.§ The Operators μ̂_x, μ̂_y, and μ̂_z



In this section, we give the commutation relations for the operators μ̂_x, μ̂_y, and μ̂_z introduced in
Eq. (<ref>) 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]    = 2 i μ̂_z, 
    
    	[μ̂_x, μ̂_z]    = 2 i μ̂_y, 
    
    	[μ̂_z, μ̂_y]    = 2 i μ̂_x.


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.<ref> of the main text), we start by writing
the full modified Hamiltonian in the frame rotating at the drive frequency ω_:


    _, mod(t) = _,0(t) + _(t) + ∑_n _^(n)(t),


where _,0(t) and _(t) are, respectively, given by Eq. (<ref>) and Eq. (<ref>) of the main text, and 


    _^(n)(t) = 2 g^(n)(t) cos(ω_ t)μ̂_x + 2 g^(n)(t) sin(ω_ t) μ̂_y + [Δ^(n) + 2 g^(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 _,0(t). The Hamiltonian _(t) is
then given by


    _,(t) = ṽ_,x(t) μ̂_x + ṽ_,y(t) μ̂_y + ṽ_,z(t) μ̂_z,


where


    ṽ_,x(t)    = f_(t) [ sin(2 ω_ t) cosh(2 r) + 2 sin(ω_ t) sinh(2 r) ],
    ṽ_,y(t)    = - f_(t) cos(2 ω_ t), 
    ṽ_,y(t)    =  f_(t) [ sin(2 ω_ t) sinh(2 r) + 2 sin(ω_ t) cosh(2 r) ].


Similarly, we find that the correction Hamiltonian in the interaction picture is given by


    _,^(n)(t) = w̃_,x^(n)(t) μ̂_x + w̃_,y^(n)(t) μ̂_y + w̃_,z^(n)(t) μ̂_z,


where


    w̃_,x^(n)(t)    = 2 g^(n)(t) cos(ω_ t) cosh(2 r) + [Δ^(n) + 2 g^(n)(t)] sinh(2 r), 
    w̃_,y^(n)(t)    = 2 g^(n)(t) sin(ω_ t), 
    w̃_,z^(n)(t)    = 2 g^(n)(t) cos(ω_ t) sinh(2 r) + [Δ^(n) + 2 g^(n)(t)] cosh(2
    		r).


For simplicity we have omitted the explicit time dependence of r, i.e.,  r = r(t) [see Eq. (<ref>) of the main text], in 
Eqs. (<ref>) and (<ref>).



 §.§ The Linear System of Equations


In this section we give explicitly the matrix elements of P_. As explained in Sec. <ref>, the matrix P_ is
analog to the matrix M introduced in Eq. (<ref>) of the main text and encodes the dynamics of the ideal evolution
generated by _,0 (t) (see Eq. (<ref>) of the man text). The difference comes from our choice of expanding in a
Fourier series the unknown envelope functions g_α^(n) (t) (α =x,y, see Eq. (<ref>) of the main
text) and the static detuning Δ^(n) instead of the functions w̃_,α (t) [α =x,y,z, see
Eq. (<ref>)]. Within this framework, the system of linear equations that allows one to determine the coefficients
defining the nth order control correction is 


    P__^(n) = _^(n),


where _^(n) = { c_x, n^(1), c_y, n^(1), Δ^(n)}^T is the vector of unknown coefficients and _^(n) is
the vector of spurious error-Hamiltonian elements. To first order we have _^(1)= - ∫_0^t_t{ṽ_,x
(t), ṽ_,y (t), ṽ_,z(t) }^T, with ṽ_,α (t) (α =x,y,z) given in
Eq. (<ref>).

Finally, the matrix elements of P_ are 


    P_,11   = ∫_0^t_t  [1 - cos(ω_1 t)] {[1 + cos(4 t)] cosh[2 r] + 2 cos(2 t) sinh[2 r] },
    	
    
    	P_,12   = ∫_0^t_t  [1 - cos(ω_1 t)] {sin(4 t) cosh[2 r] + 2 sin(2 t) sinh[2 r] },
    	
    
    	P_,13   = ∫_0^t_tsinh[2 r],
    	
    
    	P_,21    = -∫_0^t_t  [1 - cos(ω_1 t)] sin(4 t),
    	
    
    	P_,22    = -∫_0^t_t  [1 - cos(ω_1 t)] [1 - cos(4 t)],
    	
    
    	P_,23    = 0,
    	
    
    	P_,31    = ∫_0^t_t  [1 - cos(ω_1 t)] { [1 + cos(4 t)] sinh[2 r] + 2 cos(2 t) cosh[2 r] },
    	
    
    	P_,32    = ∫_0^t_t  [1 - cos(ω_1 t)] {sin(4 t) sinh[2 r] + 2 sin(2 t) cosh[2 r] },
    	
    
    	P_,33    = ∫_0^t_t cosh[2 r],


where ω_1 = 2 π/t_ and once more we have omitted the explicit time dependence of r.



§ TRANSMON QUBIT





 §.§ Interaction Picture


In this appendix we show some steps of the general method applied to the transmon qubit that were omitted in the main text for
brevity. 

We first write the full modified Hamiltonian in a frame rotating with the drive frequency: 


    _, mod(t) = _, 0(t) + V̂_(t) + ∑_n Ŵ_T^(n)(t),


where _, 0(t) is given by Eq. (<ref>), V̂_(t) is given by Eq. (<ref>), and


    Ŵ_^(n) = 1/2Δ^(n)( z + 3 |2⟩⟨2|) 
    	+ 1/2 g_x^(n)(t) ( x + ηx,12) 
    	+ 1/2 g_y^(n)(t) ( -y + ηy,12).


As we previously did for the two-level system (see Sec. <ref>) and the DPA (see Sec. <ref>), we use the detuning
as yet another free parameter in the control Hamiltonian.

Before we move to the interaction picture with respect to _, 0(t), let us adopt, for convenience, the following notation: 


    1, 2, 3  = x, y, z;  
       4, 5, 6, 7 = x,12, y,12, x,02, y,02; 
       8 = |2⟩⟨2|.


Moving to the interaction picture with respect to _, 0(t), the Hamiltonian V̂_(t) is given by


    V̂_, (t) = ṽ_, 1(t) 1 + ṽ_, 2(t) 2 + ṽ_, 3(t) 3 
    	+ ṽ_, 4(t) 4 + ṽ_, 5(t) 5 + ṽ_, 6(t) 6 
    	+ ṽ_, 7(t) 7 + ṽ_, 8(t) 8,


where

    ṽ_, 1(t)    = ṽ_, 2(t) = ṽ_, 3(t) = 0, 
    ṽ_, 4(t)    = η/2 f_(t) cos(θ/2) cos(α t), 
    ṽ_, 5(t)    = η/2 f_(t) cos(θ/2) sin(α t), 
    ṽ_, 6(t)    = η/2 f_(t) sin(θ/2) sin(α t), 
    ṽ_, 7(t)    = -η/2 f_(t) sin(θ/2) cos(α t), 
    ṽ_, 8(t)    = 0,


where for simplicity we have omitted the explicit time dependence of θ, i.e.,


    θ(t) = ∫_0^t t_1 f_(t_1).


Proceeding similarly we find 


    _, ^(n)(t) =  w̃_, 1^(n)(t) 1 + w̃_, 2^(n)(t) 2 
    	+ w̃_, 3^(n)(t) 3 + w̃_, 4^(n)(t) 4 + w̃_, 5^(n)(t) 5 
    	+ w̃_, 6^(n)(t) 6 + w̃_, 7^(n)(t) 7   + w̃_, 8^(n)(t) 8


where 


    w̃_, 1^(n)(t)    = 1/2 g_x^(n)(t), 
    w̃_, 2^(n)(t)    = 1/2 g_y^(n)(t) cosθ + 1/2Δ^(n)sinθ, 
    w̃_, 3^(n)(t)    = -1/2 g_y^(n)(t) sinθ + 1/2Δ^(n)cosθ, 
    w̃_, 4^(n)(t)    = η/2[ g_x^(n)(t) cos(α t) cos(θ/2) - g_y^(n)(t)
    		sin(α t) cos(θ/2) ], 
    w̃_, 5^(n)(t)    = η/2[ g_x^(n)(t) sin(α t) cos(θ/2) + g_y^(n)(t)
    		cos(α t) cos(θ/2) ], 
    w̃_, 6^(n)(t)    = η/2[ g_x^(n)(t) sin(α t) sin(θ/2) + g_y^(n)(t)
    		cos(α t) sin(θ/2) ],
    w̃_, 7^(n)(t)    = η/2[ - g_x^(n)(t) cos(α t) sin(θ/2) + g_y^(n)(t)
    		sin(α t) sin(θ/2) ],
    w̃_, 8^(n)(t)    = 3/2Δ^(n).


It is convenient to use the Gell-Mann λ̂ operators to calculate commutators. The Gell-Mann λ̂ operators
are given by Eq. (<ref>), except for 8, which is given by


    8 = (|0⟩⟨0| + |1⟩⟨1| - 2 |2⟩⟨2|)/√(3)
 

The Gell-Mann operators satisfy the following commutation relations:


    [a, b] = 2 i ∑_c f^abcc,


where the structure constants f^abc are completely antisymmetric in the three indices, and are given by 


    f^123 = 1,  f^147 = f^165 = f^246 = f^257 = f^345 = f^376 = 1/2,  f^458 = f^678 = √(3)/2.


The commutation relations of the Gell-Mann matrices are very convenient, specially when evaluating the Magnus expansion for this problem. 



 §.§ Choice of Free Parameters


In the main text we showed that one has seven equations to fulfill for the transmon qubit problem, and this requires at least
seven free parameters. We also commented that it is important that the envelope functions g_x^(n)(t) and g_y^(n)(t) of the
correction Hamiltonian (see Eq. (<ref>) of the main text) have a bandwidth comparable to α, so that one can
access transitions between the levels |1⟩ and |2⟩. This becomes more clear if one considers the expressions of
w̃_, j^(n)(t) in Eqs. (<ref>). One can see that  w̃_, j^(n)(t) oscillates with
frequency α for j = 4, …, 7, while w̃_, j^(n)(t) is a slowly varying function for other values
of j. Since the effect of the correction Hamiltonian on the dynamics at t = t_ is given by the integral of Ŵ_,
(t), the terms proportional to 4, …, 7 average out unless g^(n)(t) has a bandwidth comparable to
α. As a consequence g_x^(n)(t) and g_y^(n)(t) must have a bandwidth comparable to α.

Practically, this means that the envelope functions g_x^(n)(t) and g_y^(n)(t) associated to the correction Hamiltonian
(see Eq. (<ref>) of the main text) need to have a certain number of non-zero coefficients such that the condition on the
bandwidth can be satisfied. A systematic way of determining which coefficients are non-zero is to choose the coefficients of the
harmonics between k=1 and k ≃α t_ / 2 π in the Fourier expansion of the envelopes to be non-zero and set
all the other coefficients to zero. Furthermore, assuming that the detuning is time-independent, all coefficients of its Fourier
series except c_z,0 are zero. This typically gives us more than seven free coefficients in total, and we end up with an
underdetermined system of linear equations.

As mentioned in the main text, we can use the Moore–Penrose pseudo-inverse <cit.> to solve
this underdetermined system of linear equations. Importantly, the pseudo-inverse always exists, which guarantees that the linear
system always has a solution, and the pseudo-inverse also enforces that the solution has the smallest possible norm, which results
in having only a fraction of the free coefficients to actually be non-zero.



§ SNAP GATES




 §.§ Correction Hamiltonian



As discussed in the main text, the correction Hamiltonian for SNAP gates (see Eq. (<ref>) of the main text) does not
allow one to correct terms proportional to z|n⟩⟨n| using the general strategy presented in Sec. <ref>. As
we argue in the main text, the most important source of errors are precisely those originating from terms in the error-Hamiltonian
proportional to z|n⟩⟨n|. This makes it necessary to find another strategy to correct those errors. 

A correction Hamiltonian with terms proportional to z|n⟩⟨n| in the rotating frame is out of question, since it would
require a dispersive coupling constant dependent on n. We must, therefore, abandon the general strategy used so far and look
for an alternative approach. 

Let us write explicitly the Magnus expansion, up to the second order, of the evolution operator associated to the modified
Hamiltonian _mod,^(1)(t) = _,(t) + _,^(1)(t). We have 


    _1^(1)(t_) + _2^(1)(t_) =   _1^(0)(t_) + _2^(0)(t_) - i
    		∫_0^t_t_1_,^(1)(t_1) 
        - 1/2∫_0^t_t_1∫_0^t_1t_2 {[ _,(t_1),
    		_,^(1)(t_2) ] + [ _,^(1)(t_1), _,(t_2) ] + [
    		_,^(1)(t_1), _,^(1)(t_2) ] }


In the general strategy presented in Sec. <ref>, we neglect the term originating from the double integral with the
argument that it is a high order term in the perturbative series. However, if one calculates the commutators [
_mod,^(0)(t_1), _,^(1)(t_2) ] and [ _,^(1)(t_1), _,^(1)(t_2) ] (cf.
Eqs. (<ref>) and (<ref>) of the main text), one finds terms proportional to z|n⟩⟨n|. This gives us a
path forward to find a correction Hamiltonian: we will keep the higher order term originating from the double integral in
Eq. (<ref>) and try to find a correction Hamiltonian _,^(1)(t) that guarantees that _1^(1)(t_)
+ _2^(1)(t_) = 0. 

Substituting the expression for the correction Hamiltonian in the interaction picture (see Eq. (<ref>) of the main text) in
Eq. (<ref>) and expanding the envelope functions  g_x,n(t) and g_y,n(t) in a Fourier series that we truncate at
k=k_max, we get a quadratic system of equations in the free parameters that allows us to satisfy the condition
_1^(1)(t_) + _2^(1)(t_) = 0.

Solving such a system of equations is still a difficult thing to do, since we have a system of 3 N_trunc quadratic
equations depending on 4 k_max N_trunc free parameters. Here, N_trunc is the dimensionality of the truncated
cavity Hilbert space. There is, however, a convenient approximation one can do to simplify the problem: one can assume that the
effect of g_x,n(t) and g_y,n(t) on cavity levels other than |n⟩ is small and can be neglected. This allow us to
break the initial system of 3 N_trunc equations in N_trunc independent systems of 3 equations each, depending on
4 k_max free parameters only. These systems of equations have, however, several solutions since they are nonlinear. To
choose the “best” solution, it is convenient to work with more free variables than equations and use Lagrange multipliers to
find solutions that minimize the norm of the vector of free parameters. Such systems
can easily be solved with numerical methods. In this work we have solved the system of quadratic equations using the package
HomotopyContinuation.jl <cit.> available for the Julia programming language <cit.>.

Note that even in the case where solutions for Eq. (<ref>) exist, it is not guaranteed that we will be able to mitigate
the effects of the unwanted Hamiltonian _(t). This is the case when the envelopes of the correction Hamiltonian are not
“small”. In such a case, even though the first two terms of the Magnus expansion are zero at t = t_, higher order terms
become larger, spoiling the dynamics. In fact, one could ask what a “small” correction Hamiltonian means in this context, and we
do not have a straightforward answer to that. In the example that we showed in Fig. <ref> of the main text, the spectral
weight of the corrected pulse was about 1/3 of the original pulse. This is not what usually one calls “small” for a perturbative
series, but this solution does mitigate the effects of the error-Hamiltonian.

One can use the above strategy to also correct higher order errors. By including higher order terms in Eq. (<ref>), one
gets a system of higher order polynomials to solve, with substantially more free parameters. It is, however, cumbersome to handle
these higher order terms. A more convenient approach is to assume that the correction Hamiltonian is small enough such that its effect on higher
order terms in the Magnus expansion can be neglected. If this assumption is true, we can attack this problem in the following
way: assume that we have found a correction Hamiltonian ∑_n _^(n-1)(t), that corrects errors up to order
ϵ_^n-1 = (χ t_)^n-1. Thus,


    ∑_k=1^∞Ω_k^(n-1)(t_) = 0 + 𝒪(ϵ_^n).


We want to add a term _^(n)(t) to the correction Hamiltonian such that


    ∑_k=1^∞Ω_k^(n)(t_) = 0 + 𝒪(ϵ_^n+1).


We can write Ω_k^(n)(t_) in terms of Ω_k^(n-1)(t_) and _^(n)(t). If we assume that the effect of
_^(n)(t) on Ω_k^(n)(t_) for k > 2 is negligible, then 


    Ω_k^(n)(t_) ≈Ω_k^(n-1)(t_)  for  k > 2.


Substituting Eqs. (<ref>) and (<ref>) in Eq. (<ref>), we obtain (up to order
ϵ_^n)


    _1^(n)(t_) + _2^(n)(t_) + ∑_j>2_j^(n-1)(t_) = 0.


The above equation can be rewritten as


    ∑_j _j^(n-1)(t_)    = i ∫_0^t_t_1Ŵ_^(n)(t_1) 
        + 1/2∫_0^t_t_1∫_0^t_1t_2 {[ _mod,^(n-1)(t_1),
    		Ŵ_^(n)(t_2) ] + [ Ŵ_^(n)(t_1), _mod,^(n-1)(t_2) ] + [
    		Ŵ_^(n)(t_1), Ŵ_^(n)(t_2) ] }.


The sum on the left hand side runs over the Magnus terms whose leading order is ϵ_^n. Here it is useful to simply
replace the left hand side sum by a sum running from j=1 to j=2n (see Appendix <ref>). We can
then find Ŵ_^(n)(t) using the methods discussed previously. We used this method to find fourth order corrections for
the SNAP problem shown in the main text.



 §.§ Correction Hamiltonian when Adding Phases on All Bosonic Number States



In the situation in which one wants to have arbitrary phases imprinted in all bosonic number states of the truncated bosonic
Hilbert space, one can use the general method described in Sec. <ref>. 

Since we want to have phases on all bosonic number states of the truncated bosonic Hilbert space, we need to have driving
components resonant with all the frequencies ω_,n = ω_ + n χ, for n = 0, 1, …, N-1. Therefore, the
driving Hamiltonian in the lab frame is given by Eq. (<ref>) of the main text. Thus, in the interaction picture with
respect to _ (see Eq. (<ref>) of the main text), the Hamiltonian describing the dispersive coupling between
the qubit and cavity is given by _(t) = _,0(t) + _(t), where


    _,0(t) = ∑_n=0^N-1[ f_x, n(t) x - f_y,n(t) y] |n⟩⟨n|,


and


    _(t) = ∑_n=0^N-1   ∑_n'=0^N-1( f_x,n(t) {cos[χ (n' - n) t] x - sin[χ (n' - n_0) t] y}
            -f_y,n_0(t) {sin[χ (n' - n) t] x + cos[χ (n' - n) t] y}) (1 -
    		δ_n',n) |n'⟩⟨n'|,


Considering the same correction Hamiltonian as in the main text [see Eqs. (<ref>) and (<ref>)], we find that in
the interaction picture with respect to _,0(t) [see Eq. (<ref>)] _,(t) has terms proportional to
z|n⟩⟨n|. In contrast to the case we considered in the main text, where _,0(t) acts only on the subspace of
|n_0⟩ (see Eq. (<ref>) of the main text), when we drive resonantly all bosonic number states the Hamiltonian
_,0(t) [see Eq. (<ref>)] 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. 




9
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Note7()]Note7
  note Note that the terms calculated at t_2 = 0 gives rise to
  terms proportional to (ω _Vt_f)^-2. Therefore we omit these terms in Eq. .Stop
[Blanes et al.(2009)Blanes,
  Casas, Oteo, and Ros]Sblanes2009
  author author S. Blanes, author F. Casas,
  author J. A. Oteo, and  author J. Ros, title title The magnus expansion and some of its
  applications, https://doi.org/http://dx.doi.org/10.1016/j.physrep.2008.11.001 journal journal Physics Reports volume
  470, pages 151  (year 2009)NoStop
[Rackauckas and Nie(2010)]Srackauckas2010
  author author C. Rackauckas and author Q. Nie, title title Differentialequations.jl –
  a performant and feature-rich ecosystem for solving differential equations in
  julia, https://doi.org/10.5334/jors.151 journal
  journal Journal of Open Research Software volume 5, pages 15 (year 2010)NoStop
[Bezanson et al.(2017)Bezanson, Edelman, Karpinski, and  Shah]Sbezanson2017
  author author J. Bezanson, author A. Edelman,
  author S. Karpinski, and  author V. B. Shah, title title Julia: A fresh approach to numerical computing,  https://doi.org/10.1137/141000671 journal journal SIAM review volume 59, pages
  65 (year 2017)NoStop
[Ribeiro et al.(2017)Ribeiro, Baksic, and Clerk]Sribeiro2017
  author author H. Ribeiro, author A. Baksic,  and author A. A. Clerk,  title title Systematic magnus-based approach for
  suppressing leakage and nonadiabatic errors in quantum dynamics, https://doi.org/10.1103/PhysRevX.7.011021 journal journal Phys. Rev. X volume 7, pages
  011021 (year 2017)NoStop
[Moore(1920)]Smoore1920
  author author E. H. Moore, title title On the reciprocal of the
  general algebraic matrix, https://projecteuclid.org:443/euclid.bams/1183425340 journal
  journal Bull. Amer. Math. Soc. volume
  26, pages 395 (year 1920)NoStop
[Bjerhammar(1951)]Sbjerhammar1951
  author author A. Bjerhammar, @noop  title Application of
  calculus of matrices to method of least squares; with special references to
  geodetic calculations, Vol. volume 49 (publisher Trans. Roy. Inst. Tech. Stockholm, year
  1951)NoStop
[Penrose(1955)]Spenrose1955
  author author R. Penrose, title title A generalized inverse for
  matrices, https://doi.org/10.1017/S0305004100030401 journal journal Mathematical Proceedings of the Cambridge
  Philosophical Society volume 51, pages
  406–413 (year 1955)NoStop
[Breiding and Timme(2018)]Sbreiding2018
  author author P. Breiding and author S. Timme, title title Homotopycontinuation.jl: A
  package for homotopy continuation in julia, in @noop  booktitle International Congress on Mathematical Software  (organization Springer, year 2018) pp. pages 458–465NoStop


