Engineering Effective Hamiltonians
Holger Haas, Daniel Puzzuoli, Feihao Zhang, and David G. Cory
TL;DR
This paper introduces a comprehensive framework for designing effective Hamiltonians in quantum control, enabling optimized manipulation of quantum systems through gradient-based methods and demonstrating practical experimental applications.
Contribution
It provides a general method for computing time-dependent perturbation terms and their gradients, framing effective Hamiltonian engineering as a bilinear control problem with broad applicability.
Findings
Framework for arbitrary time-dependent perturbation calculations
Application of gradient methods for Hamiltonian optimization
Experimental demonstration of control engineering in quantum systems
Abstract
In the field of quantum control, effective Hamiltonian engineering is a powerful tool that utilises perturbation theory to mitigate or enhance the effect that a variation in the Hamiltonian has on the evolution of the system. Here, we provide a general framework for computing arbitrary time-dependent perturbation theory terms, as well as their gradients with respect to control variations, enabling the use of gradient methods for optimizing these terms. In particular, we show that effective Hamiltonian engineering is an instance of a bilinear control problem - the same general problem class as that of standard unitary design - and hence the same optimization algorithms apply. We demonstrate this method in various examples, including decoupling, recoupling, and robustness to control errors and stochastic errors. We also present a control engineering example that was used in experiment,…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Engineering Effective Hamiltonians
Holger Haas
Department of Physics, University of Waterloo, Canada
Institute for Quantum Computing, University of Waterloo, Canada
These authors contributed equally to this work.
Daniel Puzzuoli
Department of Mathematics and Statistics, University of Ottawa, Canada
School of Mathematics and Statistics, Carleton University, Canada
These authors contributed equally to this work.
Feihao Zhang
State Key Laboratory of Low-Dimensional Quantum Physics and Department of Physics, Tsinghua University, China
David G. Cory
Department of Chemistry, University of Waterloo, Canada
Institute for Quantum Computing, University of Waterloo, Canada
Perimeter Institute for Theoretical Physics, Canada
Canadian Institute for Advanced Research, Canada
Abstract
In the field of quantum control, effective Hamiltonian engineering is a powerful tool that utilises perturbation theory to mitigate or enhance the effect that a variation in the Hamiltonian has on the evolution of the system. Here, we provide a general framework for computing arbitrary time-dependent perturbation theory terms, as well as their gradients with respect to control variations, enabling the use of gradient methods for optimizing these terms. In particular, we show that effective Hamiltonian engineering is an instance of a bilinear control problem – the same general problem class as that of standard unitary design – and hence the same optimization algorithms apply. We demonstrate this method in various examples, including decoupling, recoupling, and robustness to control errors and stochastic errors. We also present a control engineering example that was used in experiment, demonstrating the practical feasibility of this approach.
1 Introduction
Efficient tools for engineering control sequences that drive a quantum system to undertake desired evolution are critical for quantum computing, sensing, and spectroscopy. In the case of quantum computing [1, 2, 3], it is imperative that the effective evolution corresponds, as closely as possible, to that of the experimenter’s best characterization of the system Hamiltonian, as only this can reliably lead to high fidelity unitary operations. In realistic settings this requires successful suppression of numerous unwanted, yet unavoidable, physical effects: couplings to uncharted or unaccountable external degrees of freedom [4, 5, 6], leakage out of the computational subspace [7, 8, 9], as well as uncertainties and stochastic variations in the system’s internal and control Hamiltonians [10, 11, 12]. In the case of sensing and spectroscopy [13, 14, 15, 16, 17], the experimenter is interested in letting the system evolve under some Hamiltonian that is not fully characterized, while removing the effects of other, potentially unknown and potentially much stronger Hamiltonian terms that interfere with the effects of the Hamiltonian of interest, limiting sensing ability or spectroscopic resolution. In most cases, the unwanted and wanted effects both arise from Hamiltonian terms that are either not fully characterized or cannot be fully accounted for, and hence methodologies for suppressing undesired effects – while potentially retaining detectability of others – have to rely on perturbation theory analysis.
To formalize the above, we say that our quantum system is controlled over a period , and denote the unitary evolution over this time period, as is generated by the experimenter’s best characterization of the system’s internal and control Hamiltonians, by . Successful engineering of the desired effective evolution boils down to ensuring that and variations of it with respect to particular Hamiltonian variations, take a desired form. Such variations can generally be expressed as time-dependent perturbation theory expressions of the following form:
[TABLE]
where is a scalar function and is a set of, possibly time-dependent, operators. These integrals arise in a variety of existing treatments [13, 4, 5, 18, 19] and are further discussed in the upcoming paragraphs. Control design for quantum computing implementations often requires ensuring that some list of such nested integrals are minimized or, better yet, equal to zero. This demand can occasionally be fulfilled somewhat incidentally; by ensuring that the control fields are as strong as possible the experimenter tends to minimize the control period and thereby the effect of some perturbations. Conversely, sensing and spectroscopy applications typically need control sequences that minimize some set of the nested integrals above while maximizing others, hence, the fastest control approach does not suffice.
Analytical perturbative tools for engineering effective Hamiltonians were introduced by Haeberlen and Waugh [13] with their average Hamiltonian theory (AHT). AHT prescribed a systematic approach for setting perturbation theory integrals of the kind in Equation (1) with to some desired values. AHT immediately proved an indispensable tool for the development of a vast number of magnetic resonance control sequences, e.g., dipolar sequences [20, 21, 22, 23, 24, 25], composite pulses for control and internal Hamiltonian variations [26], imaging sequences [27] and many more. Another analytical treatment was given by dynamical decoupling (DD) [4, 5, 28] and dynamically corrected gates [29, 30] introduced in the context of quantum computing. Perturbation theory terms with in Equation (1) appear when solving for the ensemble averaged evolution of a quantum system under stochastic operators, as in stochastic Liouville theory [18]. In such cases, will be composed of correlation functions that characterize the stochastic operators. Analytic control design seeking to minimize nested integrals of that kind was performed in [31].
In addition to the above considerations, achieving the most efficient and accurate control of any quantum system – or an ensemble of quantum systems – requires tailoring of control sequences for the particular experimental setup and physical system at hand. When it comes to flexible tailored control design, numerical control optimization has a number of advantages over analytical control design: (i) it can easily deal with simultaneous control of an ensemble [32, 33], (ii) it is not specific to any Hilbert space dimension, (iii) it can accommodate any experimental constraints present for the specific hardware configuration, e.g., amplitude and bandwidth constraints for the control waveform [34], (iv) it can account for deterministic control distortions due to control hardware [35], and (v) it stands a better chance of yielding control sequences that are closer to being time optimal given (iii) and (iv). Furthermore, recent technical advances such as the use of graphics processing units and automatic differentiation [36] hold promise of significantly improving efficiency, and streamlining the implementation of numerical control engineering routines.
Given the benefits of perturbative tools and numerical control design, there has been increasing interest in numerical optimization of control sequences that implement effective Hamiltonians. A filter function formalism for mitigating the effect of stochastic noise in quantum control was introduced by Green et al [19, 37], and has been combined with gradient free numerical optimization, leading to experimental advancements [38]. The filter function approach was fully generalized to be applicable to general classical and quantum noise in [39], and furthermore, a set of fundamental filter functions, out of which all other filter functions can be constructed, was identified. Although there are individual, problem specific, numerical approaches that have previously been taken [40, 41], a complete framework for numerical control optimization that would yield a desired value for simultaneously with values for an arbitrary set of perturbation terms has so far been lacking.
With this manuscript, we provide a general method for the numerical evaluation of simultaneously with the evaluation, or arbitrarily close approximation, of any number of nested integrals of the kind in Equation (1). Furthermore, this method also enables straightforward computation of gradients of these integrals, which is crucial for efficiently searching large control landscapes. We accomplish this by generalizing the work of Van Loan [42], Carbonell et al [43] and, more recently, Goodwin and Kuprov [44], who showed that certain nested integrals involving matrix exponentials can be evaluated via exponentiation of a single block matrix. This method has found application in unitary engineering [45], as it provides an accurate and efficient tool for evaluating partial derivatives of . We also note that the first order version the method outlined here has been observed in [46, 47] for evaluating functional derivatives of with respect to control amplitudes.
Our method is aimed at complementing the existing quantum control tools dealing with open quantum systems and systems interacting with non-Markovian environments which have been reviewed in [48] and [49]. We generalize the pre-existing work by extending the block-matrix methods to compute the perturbation theory terms in Equation (1) to arbitrary order, and also develop tools for approximating nested integrals involving arbitrary scalar functions . This is done by showing that the perturbation theory terms may themselves be written as parts of solutions to first order matrix differential equations, which we call the Van Loan equations. The Van Loan equations have the same form as the Schrödinger equation, and in particular depend on the control amplitudes in the same way. The immediate benefit of this formulation is that control and optimization of the perturbation theory terms is now a computational problem of the same kind as standard unitary design, and as such the same optimization methods, including those that use gradient information, can be employed. Most of this manuscript is devoted to attempting to clearly demonstrate how to exploit the differential equation formulation for the purpose of numerical control searches that involve various perturbation expressions.
Some of the authors of this manuscript have successfully employed the methods presented for nanoscale magnetic resonance imaging experiments [50]. These experiments posed a very challenging control setting – we were dealing with an ensemble of strongly dipolar coupled proton spins that experienced a vast Rabi (control) field strength () variation of ., while the phase coherence time () of the coupled spins was 11 s. Our numerical tools helped us to find control sequences that yielded a unitary rotation that was insensitive to first order perturbations due to dipolar and chemical shift Hamiltonians for the entire spin ensemble simultaneously. Even though the rotation took 7.5 s to implement, it enabled an increase of spin by a factor of . We strongly believe that such coherence time enhancements would not have been possible without the numerical tools developed here.
In this manuscript, we first give some background for matrix differential equations and effective Hamiltonians in Section 2. We then specify our general approach for tackling control problems in Section 3. With Section 4, we present a solution for a general time dependent upper triangular block matrix differential equation and highlight how it can be used for calculating nested integrals in Equation (1). Subsequently, we exemplify the construction of Van Loan block matrix differential equations and numerical control optimizations with five examples in Section 5, which include the control sequence we engineered for the aforementioned nanoscale magnetic resonance experiments that was optimized to be implemented in the presence of a non-trivial transfer function for the control hardware.
2 Effective Hamiltonians
We denote the set of complex matrices by . The starting point for effective Hamiltonian analysis is an initial value problem (IVP) of the form:
[TABLE]
where are matrix valued functions, the initial value is , and denotes the time derivative of . In the context of quantum control, we will typically have , for a time dependent Hamiltonian, but could also represent the generator for a master equation, and in any case it is notationally convenient to consider a general . We call the generator of the above IVP, and the propagator. Under assumptions on the generator (which we will not explicitly state or worry about) this IVP has a unique solution [51], which we will write using the time-ordered exponential notation:
[TABLE]
In this manuscript we will not work with the “time-ordering” operator, we simply regard the above expression as a choice of notation for the solution of the above IVP.
The goal of any effective Hamiltonian treatment such as AHT, DD or filter function formalism is to analyze the effect that a variation in generator has on the propagator. Formally, for two functions we want to analyze how the evolution of a system with generator is different from a system with generator , where we are viewing as a variation of the generator .
The toggling frame provides a way of writing the propagator of a system evolving under in a way that clearly separates out the deviation caused by .111The toggling frame concept was utilized in [52], though our presentation more closely follows that of [13]. Neither of these references use the terminology of toggling frame, which appeared later (see, e.g., the presentation of [14]). We denote the propagator under alone:
[TABLE]
the propagator under both:
[TABLE]
and the toggling frame propagator, defined as:
[TABLE]
where . With these definitions, it holds that
[TABLE]
which may be verified by differentiating both sides of the equation and verifying that they are solutions to the same IVP.
The decomposition packages all variation of as an effect of into , and hence, the deviation of from caused by may be analyzed by studying . Operating in the perturbative limit, effective Hamiltonian schemes analyze via series expansion, either through the Dyson series [53]:
[TABLE]
or via the Magnus expansion [54, 55], which under certain conditions gives for
[TABLE]
where denotes the matrix commutator.
How robust a control sequence is to a variation is then analyzed perturbatively using one of the above expansions. Furthermore, robust control sequences are designed specifically to optimize the above terms.
2.1 General Form of Perturbation Terms
In this manuscript, we will be concerned with integrals of the form
[TABLE]
for , and where is some scalar valued function. As a shorthand, we denote the above integral as , and when (i.e. is a constant), we will write . should be read as Dyson term.
Terms arising from either the Dyson series or Magnus series may be constructed out of integrals of the above form. In application, the function will often be a correlation function of a time-dependent stochastic noise source.
2.2 The Dyson Series and Directional Derivatives
In this manuscript, we will use the Dyson series expansion, as its terms have a direct interpretation as directional derivatives. As an example, we consider the directional derivative of as a result of variation in in the direction , given by
[TABLE]
If we expand via the Dyson series, the result is a power series for in :
[TABLE]
and from this we may directly read off the directional derivative as the matrix corresponding to the term:
[TABLE]
Similarly, the second derivative is
[TABLE]
The same analysis applies with respect to multiple variations:
[TABLE]
from which we may conclude that
[TABLE]
Within the context of quantum control, directional derivatives with respect to variations in the generator are typically viewed in two ways: (1) when arises from a variation in an underlying control sequence, the directional derivatives are used to make informed decisions about how to modify the control sequence to make it better, and (2) when represents uncertainty in , or even known but unwanted terms in the generator, these terms represent how robust the propagator under is to the variation .
3 Setup of the Control Problem
In this section, we give a full, albeit abstract, description of the control problems the block matrix Van Loan differential equation framework is capable of addressing. The description that will be outlined maps almost one to one to our implementation. In fact, a lot of our treatment and notation has been chosen specifically to simplify the implementation process while retaining full generality. In Figure 1, we illustrate the general control setting addressed. We say that we have a finite set of quantum systems labelled by a single compound label . In principle, this is true for any control setting, although in many practical cases one approximates macroscopic ensembles of quantum systems as being parametrized by some set of continuous variables – Rabi field strengths, resonance offsets, etc – in such cases, we think of as a representative sample of the real ensemble. The ensemble could in some cases denote the same quantum system under different conditions for distinct experimental realizations, i.e., it might stand for an ensemble in time rather than a spatial ensemble of physical systems.
The ensemble of quantum systems is controlled by a control sequence – a real vector valued function specified over an interval and delivered by some control signal source. The source should be thought of as a physical device which outputs , ; usually, we think of it as an arbitrary waveform generator. In the context of this manuscript, we regard as the waveform generated by our numerical pulse search routines. Each quantum system labelled by has an associated transfer function , , where we use to denote the space of real vector valued functions . is an analytic deterministic map which transforms the control sequence to system specific control amplitudes , . The components of are the real valued functions that appear in the matrix differential equation determining the evolution of system . In Appendix B, we demonstrate how to construct for piecewise constant control sequences and control amplitudes in the case of linear transfer functions.
All quantum control problems boil down to engineering quantum state trajectories with certain desired properties. Mathematically this corresponds to generating a system propagator , , which satisfies some set of conditions. We emphasize that the properties wanted from need not be merely its value at time , they could also be various integral expressions of over , which describe the trajectory of . Here, a quantum system should be understood simply as a finite level system or one that can be treated as such; the time dependent state of the quantum system is determined by . The time dependent value of itself is determined by a first order linear matrix differential equation which we refer to as the system differential equation:
[TABLE]
where is the system generator, while . The value of at each instant is determined by the control amplitudes
[TABLE]
where is a constant matrix for all . This implies that the problem is a bilinear control theory problem [56]. The system differential equation should be understood as the Schrödinger equation or some generalization of it, e.g., the Liouville-von Neumann equation for vectorized density matrices.
We note we assume all to have identical generators for all . Even though this may not be the case for all quantum control problems, one can always use our problem description by employing a direct sum of different sets of . Such an approach is computationally not the most efficient, but it does substantially simplify implementing the algorithm while retaining total generality.
As we said in Section 2, Equation (17) has a formal solution
[TABLE]
We are typically interested in finding an , , that yields the wanted final unitary operations as well as some desired values for nested integral expressions of the following form:
[TABLE]
In this manuscript, we designate as a set of block matrices composed of complex matrices, hence, an element of is an element of . The main and the most significant result of this manuscript is demonstrating that control problems of this kind can still be written as bilinear control theory problems [56] that involve the same control amplitudes that appear in the system differential equation. In order to find that yields the desired and the desired values for any set of integral expressions for , , we can always construct a block matrix differential equation – called the Van Loan differential equation – which comprises the system generators and the objects that appear in the integral expressions for . The Van Loan differential equation is expressed as
[TABLE]
where is the Van Loan propagator and is the Van Loan generator
[TABLE]
for . It will be shown that the integral expressions of interest appear as various blocks of . The benefit of such block matrix methods is two-fold: it enables an accurate and efficient way for evaluating the integral expressions for piecewise constant and is readily deployable within control finding routines that take advantage of the linear differential equation structure of the problem.
Having constructed the Van Loan differential equation that enables the evaluation of all terms of interest, we can always define a target function for each system in the ensemble. being a function of the final Van Loan propagator for system , i.e., , where corresponds to having the desired properties from the system propagator and from any number of nested integral terms of interest. Finally, we combine into a target function for the whole ensemble : , where are the relative weights assigned to each member of . We have assumed that , for all , and that . Of course, the linear form of is not necessary but it does simplify the implementation. It is clear that is a functional of and its derivatives with respect to the control sequence are given as
[TABLE]
Throughout this manuscript, we will only deal with piecewise constant control amplitudes , for which, we split the interval into subintervals with respective durations such that , for all , and . We say that the control sequence takes a constant value over each of the subintervals, i.e.,
[TABLE]
where is a real valued matrix that contains all piecewise control elements . Given a non-identity transfer function , which is likely the case for any experimental setting, we have to use the chain rule to evaluate . For piecewise constant control settings, we define matrices that specify the piecewise constant amplitudes of just like for above. We split the interval into subintervals with respective durations such that , for all , and . We note that does not necessarily have to match . We can now specify the components of which correspond to the piecewise constant values of :
[TABLE]
Because we treat the control sequence and the control amplitudes as piecewise constant functions, we will also, from now on, regard the transfer functions as matrix valued functions, i.e., , such that . We write . Solving the Van Loan differential equations enables us to find the values of , as is shown in Appendix A. Hence, we can evaluate . Since we assume to be an analytic matrix valued function, we can always evaluate the elements of its Jacobian in order to implement the maximization of . For all examples considered in this manuscript are taken to be linear, meaning that their Jacobians are trivial.
Finally, in most practical cases the experimentalist needs to adhere to some constraints, e.g., pulse waveform bandwidth and amplitude constraints or periods for which . As constrained optimization is more technically challenging than unconstrained optimization, it is common to try to enforce the constraints on in a way that keeps the overall optimization problem unconstrained. Our preferred method for doing this is an optimization transfer function , , where is the number of time steps for the piecewise constant optimization control sequence . The idea is that the mapping is constructed in a way that ensures all control sequences in its output space adhere to either all or some of the constraints. Having constructed , the numerical pulse searches are then carried out over its input , which need not adhere to all constraints on . We represent the piecewise-constant by a matrix , just as we did for and above. For finding a suitable control sequence using gradient based algorithms, one then needs to evaluate
[TABLE]
which means evaluating the elements of the Jacobian for as well as . After finding an that yields a high enough value, the control sequence that is to be implemented experimentally is calculated simply as . In Appendix B, we demonstrate explicitly how to construct that introduces zero pulse amplitudes to the beginning and the end of the control sequence and how to construct that limits the bandwidth of the waveform in the frequency (Fourier) domain. We note that our use of optimization transfer function is similar to the method in [57]; i.e. it implements the mapping from the parameters of a control sequence to the control sequence itself.
4 Computational Methods
In this section, we outline the framework for computing Dyson terms of general form, , defined in Section 2.1. The general idea is that these terms may be written as solutions of first order matrix differential equations of the same form as the base differential equation. The generators for the new differential equations are block matrices with blocks consisting of pieces from the original differential equation.
To illustrate the approach, we consider the simplest case, a first order integral:
[TABLE]
for . We may further simplify this by assuming that and , i.e. they are time-independent. In this case, the expression reduces to
[TABLE]
A priori, computing the above for a particular time requires an integral approximation method. However, it was originally observed by Van Loan [42] that this expression can be computed using a single matrix exponential:
[TABLE]
Van Loan showed [42] more generally how nested integrals up to order involving matrix exponentials can be computed by exponentiating a single upper triangular block matrix, and [43] extended this to arbitrary order.
This has found application in physics where such expressions often arise [58, 44], and in particular it has been used to compute directional derivatives for pulse finding [45]. Here, we extend this idea to the case of time-dependent matrices, to compute integrals involving time-ordered exponentials. The simplest case of this extension is the time-dependent version of Equation (29). For two matrix-valued functions , and , it holds that
[TABLE]
The above formula may be verified by differentiating both sides of the equation, and verifying they are both solutions to the same initial value problem. Hence, we may compute via propagation of the differential equation
[TABLE]
We note that this particular formula has been observed in [46] in the context of pulse finding for derivative evaluation.
In this section, we present a generalization of the previous work to arbitrary order in the time-dependent case, including scalar functions . The general idea is the same as for the first order example; a Dyson term may be rephrased as part of the solution to a first order matrix differential equation. In the context of control, this rephrases controlling Dyson terms as a bilinear control theory problem [56]. In all cases, the generator of the differential equation is an upper triangular block matrix, and hence we also develop general tools for analyzing the structure of the time-ordered exponential of arbitrary upper triangular block matrices.
This section is organized as follows.
- •
In Section 4.1, we generalize the theorems in [42, 43], giving the general structure of time-ordered exponentials of upper triangular block matrices. As described therein, Appendix C describes code for symbolically simplifying this structure
- •
In Section 4.2, we give a differential equation computing , i.e. the case when no scalar function appears in the integral.
- •
In Section 4.3, we provide a similar construction for terms when is either a linear combination of exponentials, or is a polynomial.
We note that in this manuscript we are concerned specifically with terms arising in effective Hamiltonian treatments, but Section 4.1 describes a much more general class of integrals involving time-ordered exponentials that this approach may be applied to. Hence, this method may find application in control design beyond optimization of Dyson terms.
4.1 Integrals Involving Time-Ordered Matrix Exponentials
Here, we present a full time-ordered generalization of the theorems of Van Loan [42] and Carbonell et al. [43]. First, we introduce some notation. Let for . For , denote
[TABLE]
For and indices , denote
[TABLE]
and for a single index , denote
[TABLE]
Note that for and indices , these definitions satisfy the recursive expression
[TABLE]
Theorem 1**.**
Let for . Define implicitly by the equation
[TABLE]
and assume that the are such that the solution to the IVP associated to the above time-ordered exponential exists and is unique222We omit an explicit statement of conditions under which existence and uniqueness holds. As we are concerned only with applications of these expressions in physical settings, finding the most general and exact technical statements for which this assumption holds is of no real interest; most reasonable physical assumptions, such as piecewise continuity, suffice (see, e.g. Theorem 3.1 in Section I.3 in [51])..
For all , , and , it holds that
[TABLE]
and
[TABLE]
where the inner sum is over all indices satisfying the relations, and and Int are as defined before the theorem. Alternatively, these matrices can be given recursively as
[TABLE]
The proof is given in Appendix D. In Appendix C we describe code that symbolically simplifies the structure arising from this theorem. That is, in general the above expressions are quite complicated, but in the constructions we will see in the following sections, many of the blocks will be [math] or proportional to the identity, in which case the expressions of the above theorem can simplify dramatically.
4.2 The Case
First, we show how to compute expressions of the form
[TABLE]
where , i.e., perturbation theory terms without a time-dependent scalar function. In this case, we may observe that for
[TABLE]
it holds that
[TABLE]
That is, the generator for this system is in , where all blocks are [math] except:
- •
All diagonal blocks are , and
- •
The first off-diagonal is given by .
The time ordered exponential has upper triangular structure with:
- •
All diagonal blocks are , and
- •
For , the block is given by .
Hence, propagating the differential equation associated with the generator computes the desired term , as well as many other terms that will likely be of interest.
To see this, one may simply apply Theorem 1 to the generator . Alternatively, one may purposefully construct this differential equation using the procedure described in the next section.
4.3 Including Scalar Functions
Next, we consider integrals of the form
[TABLE]
for , where is a scalar valued function, which may represent, for example, a correlation function for stochastic noise.
For a given , it is not immediately clear how to write the integral in Equation (50) as a part of the solution to a linear matrix differential equation, in the way we have done in the case. Certainly, it will not be possible for most functions. However, we will show here that it is possible for a very large class of functions; in particular satisfying the following properties:
- •
is a linear sum in product form: , with
- •
Each function is drawn from a finite dimensional vector space of functions closed under differentiation.
Note that polynomials and linear combinations of products of exponentials fall into this class, and we will cover these particular cases in this section333In general, the second point restricts the functions to be linear combinations of functions of the form , for a natural number and an arbitrary complex constant, i.e., linear combinations of products of polynomials with exponentials.. These special cases have the benefit that they can approximate arbitrary continuous functions. In experiment, we will take them to approximate correlation functions, or, as the correlation functions themselves arise from fits of experimental data, we could simply fit a function in these classes to the data directly.
Here, we outline a procedure for constructing a Van Loan differential equation to compute for drawn from the above class. Note that, for functions of product form, one approach is to simply absorb into the definition of and apply the method from the original case. This is valid, however it will generally introduce explicit time-dependence into . In quantum control problems, will usually only depend on time as a function of the control amplitudes, and it is computationally preferable that the generators in the newly constructed Van Loan differential equations also depend on time only through the control amplitudes.
For satisfying the above conditions, we construct a Van Loan differential equation to compute using the following algorithm.
Define a variable . This is the term we wish to compute. 2. 2.
Differentiate with respect to time. The result will be a linear combination of expressions of the same form as the original integral. 3. 3.
Add any newly appearing expressions into the list of variables. 4. 4.
Differentiate the new variables from the previous step. 5. 5.
Repeat steps and until no new expressions appear.
The assumption that the function pieces are drawn from a finite dimensional vector space of functions closed under differentiation ensures that this procedure terminates after a finite number of steps. Once the procedure terminates, we write down the resulting coupled differential equation for the defined variables. The generator for this differential equation will be an upper triangular block matrix, i.e., the generator is of the form amenable to analysis via Theorem 1.
We do this procedure when the are exponentials, and when they are polynomials.
4.3.1 Products of Exponentials
First, consider the case
[TABLE]
where . That is, is a product of exponentials in each time variable. Hence, the goal is to write
[TABLE]
where , as part of the solution to a linear matrix differential equation.
To do this, we follow the algorithm constructing Van Loan differential equations given at the beginning of Section 4.3. First, we denote the function:
[TABLE]
Differentiating, we find
[TABLE]
The new expression appearing here is the second term in the brackets. Hence, we define this as a new variable:
[TABLE]
Next, differentiating , we find:
[TABLE]
and again we define a new variable for the newly appearing term:
[TABLE]
Continuing this procedure until no new variables appear results in the following family of functions:
[TABLE]
which evolve according to the coupled differential equations:
[TABLE]
with initial conditions , and . Note that the generator for this system has upper triangular block form. In particular, the generator lies in and has all blocks equal to [math] except:
- •
The diagonal is given by , and
- •
The first off diagonal is .
For example, when , we have:
[TABLE]
Hence, if we take the time ordered exponential of the above generator, the desired integral will be in the top right block. Denoting the generator as , we may also explicitly determine the blocks of the time ordered exponential using Theorem 1:
[TABLE]
4.3.2 Polynomials for Second Order Integrals
Next, we consider polynomials, and in particular exhibit the procedure for second order integrals involving polynomials. That is, second order integrals of the form
[TABLE]
where is a polynomial in and with either real or complex coefficients, and again with . Let and be the respective highest powers of and occurring in , so that it may be decomposed as:
[TABLE]
With respect to this decomposition, the integral becomes the linear combination
[TABLE]
Here, we show how the terms for all and can be computed using a single Van Loan differential equation. Hence, all terms with a polynomial of degree at most in and in may computed using this single generator. We construct this generator by applying the procedure to the highest order term , and find by chance that the solution contains all terms of lower order. In particular, the solution to the Van Loan differential equation for computing contains a basis for the vector space of expressions
[TABLE]
Note that we have not proven this, but we conjecture it to be true, and have computationally verified this conjecture for all pairs .
Applying the procedure for generating the Van Loan differential equation to the term , we arrive at the following set of functions:
[TABLE]
This set of variables evolves in time according to the following first order coupled differential equation:
[TABLE]
with initial conditions of all variables being [math] at other than . Again, note that the derivative of each function only depends on the functions coming before it in the ordering
[TABLE]
and hence the generator for this system, which we denote is an upper triangular block matrix in . An explicit description of how to construct is:
- •
Every diagonal block is ,
- •
For the first off diagonal, the first blocks are [math], and the remaining blocks are
[TABLE]
and
- •
The off diagonal is given by repetitions of , then repetitions of , followed by zeros.
By construction, the upper right block of is . Furthermore, we conjecture the following.
Conjecture 1**.**
It holds that the top right blocks of
[TABLE]
is a basis for the vector space of expressions
[TABLE]
To get a sense for this claim, examine the special case . By applying Theorem 1, we find the top blocks of are given by:
[TABLE]
and it can be checked that these blocks form a basis for the desired set.
We have computationally verified this conjecture for all pairs the details of which can be found in Appendix C.
5 Examples
With this section, we give five examples of increasing complexity for numerical engineering of effective Hamiltonians using the Van Loan differential equation framework. First, we set up two rather standard decoupling problems with known analytical solutions and arrive at control sequences which resemble ones that have been known for some time. Our aim is not to reiterate these solutions, rather it is to demonstrate that the length of the control sequences found using block matrix numerical tools does not significantly exceed that of the sequences that have been derived analytically based on physical insights. This demonstration provides an encouraging starting point for employing the same tools to tackle far harder control problems, for which the search of analytical solutions is intractable. With the third example we provide an illustration for a problem that demands simultaneous minimization of some Dyson terms, while preserving or maximizing other Dyson terms. Such control problems are very common in many sensing and spectroscopy applications. For the first three examples we apply no transfer functions, ensemble effects nor pulse waveform bandwidth constraints etc. Only maximum amplitude constraints are used which generally yield pulses of rather jagged form, however, such constraints are typically the only constraints considered when deriving analytical solutions. With the fourth and the fifth example, we give two experimentally realistic control search examples for which we demand that the pulse ends go smoothly to zero and that the spectral width of the pulse waveform be limited. The fifth example also employs a non-trivial set of experimentally determined transfer functions .
First, we introduce the matrix norm and the matrix fidelity function that are used throughout this section. We take to stand for the Frobenius norm [59] for a matrix , defined as . We also define a fidelity function for a pair of matrices : . Furthermore, we introduce a shorthand that makes the definitions of our target functions in this section more concise. For any nested operator integral , we denote its maximum Frobenius norm maximized over all permissible control sequences , , as . For the numerical pulse searches, we always use a target function that is a linear combination of different matrix norms for various Dyson terms and the fidelity of with the target unitary . A typical target function takes the following form:
[TABLE]
where , for all , such that . denote various weights of constituent optimization targets. It is easy to see that , for all , and equal to if and only if as well as , for all . The linear form of the target function in Equation (73) is, of course, not strictly necessary, but it does greatly simplify some calculations.
All control searches were undertaken using the modified GRAPE algorithm for evaluating partial derivatives with respect to piecewise constant control amplitudes, the details of which are given in Appendix A. Our general procedure for finding a working control sequence is to first pick a pulse length and thereafter a number of time steps . We kept all subintervals of of equal length and always picked an for which , where is the length of a Rabi cycle. For each example we picked , appearing in Equation (73), that yielded relatively equal convergence rates for each constituent of the total target function during the control optimization procedure. Since this was control problem specific we generally determined the particular weights by running the optimization for a large number () of random initial seeds and monitoring the convergence of each constituent target for a given . We adjusted these weights until we observed roughly equal convergence from all constituent targets. For the examples in this section, this procedure yielded relatively equal weights for the Dyson terms in Equation (73) and a value slightly lower than . We state the specific target functions used for our optimizations in the upcoming subsections.
Given some and and having determined the suitable , we generated seed waveforms , the pulse amplitudes of which were drawn from independent uniform distributions. We used Mathematica’s FindMaximum function for multivariate conjugate-gradients optimization on these seeds with the maximum number of target function evaluations set to 1000. If none of the searches yielded we increased and and repeated the same procedure. For the seeds that reached in under 1000 evaluations, we let the optimization run until FindMaximum was terminated by machine precision.
5.1 Dipolar Decoupling
The simplest numerical effective Hamiltonian example that we consider is the problem of dipolar decoupling. For such a problem, we imagine a pair of spins coupled via the dipolar Hamiltonian . Here, we use no optimization transfer function and we assume the ensemble size to be equal to one with the transfer function acting as an identity, i.e., . We take the spin control to be global, such that the Rabi fields for either spin are identical. Also, for notational convenience, we assign and such that the two/̄spin system generator becomes
[TABLE]
which generates the following system propagator . The control task is to engineer a sequence that enables the spins to evolve effectively uncoupled. First order perturbative solution to the problem dictates setting the first derivative of in the direction of to zero, i.e., we want . The dipolar decoupling problem is a simple yet non-trivial problem. It is easy to show that the desired , , that yields cannot be generated with either for all or for all .
Here, we are not concerned about the final unitary generated on either of the spins, as our first and foremost aim is to demonstrate an ability to engineer control sequences that yield desired values for various Dyson terms. With the subsequent examples we tackle the problem of simultaneous engineering of various Dyson terms and final unitaries. Accordingly, our target function for the optimization is
[TABLE]
where the denominator is a normalization factor ensuring that . We now set up a block matrix differential equation for , that will be used for evaluating . It easy to show either by differentiation or by employing Theorem 1 that
[TABLE]
Consequently, our target could also be given as
[TABLE]
Following the procedure outlined in Appendix A, we can evaluate as a function of along with the partial derivatives with respect to the piecewise constant control amplitudes .
Given and , we can write
[TABLE]
for all and . We employ these partial derivatives in the optimization protocol described above. For the elements of we use amplitude constraints: , for all and , to ensure that , for all . We find a pulse with a total length of consisting of subintervals shown in Figure 2(a) which yields . We point out the rough similarity to the dipolar sequence proposed by Mehring [23] consisting of two square pulses with orthogonal phases depicted in Figure 2(b). The similarity is noticeable after a particular basis rotation, described in the figure caption, is performed on our numerically found sequence which is only times longer than the known analytical solution.
5.2 Universal Decoupling with Control Variations
With this example we turn to single spin control. Again, we use no optimization transfer function and we assume with the transfer function acting as an identity. The system generator is
[TABLE]
and generates a system propagator . We consider a universal decoupling sequence which would decouple all non-identity operators acting on a single spin, which translates to setting the respective lowest order directional variations of to zero, i.e., . In addition, we also demand that the sequence is robust against variations in the control control amplitudes, such that and . A sequence know to have such properties is called an XY8 sequence [60], which implements an identity operation. To demonstrate the ability of our numerical control searches in incorporating averaging of time dependent operators, we set up a search that would simultaneously set while implementing an identity operation. We search for a pulse with the following target:
[TABLE]
A suitable Van Loan generator for such a target function is
[TABLE]
that generates the following Van Loan propagator:
[TABLE]
Note that, we have not specified the elements that are not relevant for our control problem. Given the above, can be determined as
[TABLE]
For the elements of we again use amplitude constraints: , for all , ; and using our gradient optimization scheme, we find a pulse with a total length of consisting of subintervals with its pulse metrics given as: , , , , , . We present the pulse waveform in Figure 3(a). We note that, while the pulse found is times the length of the XY8 sequence, it does display definite similarities to the latter shown in Figure 3(b) after a particular basis rotation, described in the figure caption, is performed.
5.3 Exchange Interaction Recoupling
With the following example, we wish to highlight that the block matrix method does not only enable the removal of unwanted Hamiltonian terms; it is equally easy to set up optimization targets which retain or reshape parts of the Hamiltonian, while possibly removing others. A problem which arises in many situations involving ensembles of spins, is removing pairwise dipolar interactions between all members of the ensemble as well as inhomogeneities in their energy level splittings, while retaining exchange interaction with some other system or systems the spins are interacting with [61, 62, 63].
Such a situation would be described by the following Hamiltonian:
[TABLE]
where and denotes a sum over all spin pairs. The first sum in Equation (88) specifies these spin dependent energy level splitting variations , the second sum gives all dipolar interaction strengths that correspond to the dipolar Hamiltonian for a pair of spins. The third sum contains the aforementioned exchange interactions , that can often be substantially weaker than the other two terms, yet this is frequently the term in the Hamiltonian that leads to desirable spin dynamics.
Once again, we use no optimization transfer function and we assume with the transfer function acting as an identity. Here, we consider two system propagators and that are generated by and defined by Equation (83) and (74), respectively. Our block matrix tools enable us to search for control sequences that would effectively remove the spin-spin dipolar couplings and variations in level splittings, while retaining the form of the exchange interaction and performing an identity operation. To achieve this, we need to set , and , where is a constant. In order to set up a target function that can reach its maximum value, we do not set up the optimization with any specific value for . Instead, we merely enforce that the integral is proportional to . This is achieved by demanding that is orthogonal to and , i.e., . Correspondingly, the optimization target for this problem is
[TABLE]
Staying completely consistent with the formalism of Section 3 in the case of current example would require expressing , for , where denotes partial trace over the (identity) operator acting on the second Hilbert space. Such identification would enable evaluating all terms in Equation (89) by propagating a single , the diagonal elements of which are all whereas the first off diagonal elements are , and . Nevertheless, in order to reduce the dimension of the Van Loan generator , and hence the computational cost of evaluating , instead we define a block matrix generator , which decomposes into a direct sum of and and helps us evaluate :
[TABLE]
which generates
[TABLE]
Simplifications and reductions of matrix dimension of this kind are frequent and can provide significant speed ups for the numerical pulse engineering routine. We will employ a similar simplification for the example in Section 5.5. In the following, when writing our target function in terms of , we will slightly abuse our notation for specifying the components of a block matrix. We take to mean the th row and th column of as it is specified above. However, note that not all blocks of are of the same size, e.g, while . It can now be seen that
[TABLE]
For the elements of we again use amplitude constraints: , for all , ; and using our optimization scheme, we arrive at a control sequence with a total length of and subintervals with the following characteristics: , , , , , . The sequence is presented in Figure 4. The pulse does virtually remove the dipolar and Hamiltonians while rescaling the exchange coupling by a factor of with extremely small () unwanted orthogonal components.
5.4 1/f Noise Decoupling
In this section we demonstrate our ability to engineer control sequences that are designed to decouple stochastic noise characterised by its spectral density function. For that, we employ the tools for that were developed in the previous section for evaluating generalized nested integrals of Equation (50) kind. We pick 1/f noise due to its ubiquity in solid state devices, including superconducting qubits. Noise spectroscopy experiments on flux qubits [64] have clearly revealed a 1/f-like spectral density function for the level splitting variations. We proceed by evaluating the first non-zero term in the perturbative cumulant expansion for the Liouville-von Neumann equation. We then demonstrate how such a toggling frame term can be approximated and consequently minimized using Van Loan differential equations.
When treating the evolution of quantum systems under stochastic operators one has to consider ensemble behaviour averaged over many realizations of the noise process, consequently it is necessary to work with some form of Liouville-von Neumann equation [18, 10, 19, 31, 19, 39], which determines the evolution of density matrices, as the system differential equation rather than the Schrödinger equation. Accordingly, we start with a single qubit Liouville-von Neumann generator which includes a stochastic noise term and dictates the evolution of the system
[TABLE]
where , for , and is a stationary, zero mean, Gaussian stochastic function capturing the fluctuations in the qubit level spacing. This implies that , where the angle brackets denote an ensemble average over noise realizations. Here, we take the power spectral density of to be given by , where and are the smooth low and high frequency cutoffs for , respectively. It is easy to show that . For this example we use and . According to the Wiener-Khinchin theorem:
[TABLE]
where stands for the exponential integral function.
Using stochastic Liouville theory [18, 31], we can treat the noise perturbatively and show that the toggling frame propagator is given as
[TABLE]
where and we have made use of the fact that
[TABLE]
We remark that is not a unitary matrix, in fact, it is precisely the operator that encapsulates the non-unitary decoherence effects induced by . In order to reduce such decoherence at its lowest perturbative order, we need to minimize the nested double integral term in Equation (95). Because the generator is Liouville-von Neumann generator, it also holds that , where the bar denotes entry-wise matrix conjugation, and .
To employ the method developed in the previous section, we first note that a linear combination of exponential functions provides a good approximation for over a region of integration , with . Using a least squares fit to the correlation function over that region, we find an approximation that combines seven exponential functions
[TABLE]
where , , , , , , , , , , , , and . Given the approximation, we can now write
[TABLE]
It is important to realize that in the case of actual experimental scenarios either the noise correlation function or its power spectral density would be characterized before embarking on control engineering. In such cases, is extremely unlikely to fit to some simple and specific analytic function. Therefore, our procedure, for fitting a set of functions to a set of data – in this case an analytic function – in order to approximate the noise correlation, closely matches a real control engineering protocol.
We search for a pulse implementing a Y gate, i.e., we wish to set . Consequently, we use the following target function:
[TABLE]
Combining Equation (98) with Equation (61) in the previous section, we set up a Van Loan differential equation for . will be generated by
[TABLE]
such that
[TABLE]
We can now approximate the target as a function of :
[TABLE]
where is evaluated numerically for any particular . The partial derivatives of with respect to the control amplitudes are given as
[TABLE]
Here, we attempt to closely mimic a control search procedure that would be undertaken when searching for an experimentally implementable sequence. Hence, we impose three distinct constraint on the pulse waveform: maximum amplitude constraint, bandwidth limitations for the pulse waveform frequency components and zero amplitude periods at the beginning and at the end of the sequence. The last two constraints are implemented by introducing an optimization transfer function , as was described in Section 3, and the explicit construction of is given in Appendix B. We do not consider any ensemble effects, i.e., , and we take the only experimental transfer function to act as an identity, such that . Accordingly, the numerical control searches are conducted for , with .
For the searches, we limit the Rabi frequency to be less than or equal to by enforcing that
[TABLE]
for . We take the pulse length to be Rabi cycles or , which is divided into intervals of equal length , whereas the zero amplitude periods at the beginning and at the end of the pulse have a length of , corresponding to . Therefore, the numerical control search is conducted on /̄long that is divided into equal steps. Using , we constrain all spectral components of the pulse to lie within a bandwidth around the carrier frequency. The optimization transfer function , that is employed to enforce the constraints, is defined by Equation (168) in Appendix B.
Because our control optimization was carried out on the optimization control sequence we needed to evaluate for and . With Equation (118) we evaluate , and identify that for all and . Finally, we link with through Equation (170) in Appendix B:
[TABLE]
and
[TABLE]
Given the target function and the partial derivatives above, we search for a control sequence as it was described at the beginning of the section. The resulting waveform is shown in Figure 5 and the pulse characteristics are and . In the case of a stochastic operator , one cannot expect to be able to set the integral in Equation (98) equal to zero, since the high frequency components of the noise always retain their decoherence inducing effect. Nevertheless, for a reasonably low amplitude noise, our sequence in Figure 5 would extend the qubit coherence time by a factor of .
5.5 Broadband Dipolar Pulse
The control sequence presented here was engineered for nanoscale nuclear magnetic resonance experiments [50] and was used to increase the proton spin phase coherence time by a factor of 500 under rather difficult control conditions. The resonant control fields for the experiment had a vast range of dependent maximum values . Furthermore, the the strong dipolar interactions between the proton spins as well as chemical shifts limited the spin coherence time to . Furthermore, it had been determined that the amplitude and phase transfer functions for the electronics – and , respectively – had non-trivial character as it can be seen in Figure 6.
We will define our control problem exactly according to the ensemble control setup abstraction that was laid out in Section 3. We say that we have an ensemble of proton spins labelled by ; here, we consider this ensemble to be a representative ensemble of spins in the sample volume of interest. Each has an associated unique transfer function that determines the control amplitudes as a function of the control sequence . We constrain the maximum amplitude of the control sequence to be equal to one. Our transfer functions reflect the effects of Rabi field distribution and the amplitude and phase transfer functions shown in Figure 6.
In addition to the transfer functions defined by Equation (122), we also use an optimization transfer function in order to limit the range of frequency components in to within a bandwidth of as well as to enforce that the pulse starts and ends with zero amplitude. The control sequence has piecewise constant pulse amplitudes for equal periods of duration , so that the total pulse length is . The low-pass filter that is incorporated into the optimization transfer function is illustrated by the shaded area in Figure 6(a). is defined by Equation (168) in Appendix B. The control sequence is then determined by , where , and determines the number of zero amplitude intervals of duration at the beginning and at the end of the sequence.
It is natural to label the transfer functions according to the maximum Rabi field strengths they yield on the nuclear spins, i.e., according to the value corresponding to for all . The transfer functions could then be written explicitly as
[TABLE]
where the function is given by Equation (154) in Appendix B; and being evaluated by interpolating the experimental transfer function data in Figure 6.
The system generators are given by
[TABLE]
while the system propagators evaluate to for all . The control amplitudes that are specified by a matrix for each equate to
[TABLE]
and
[TABLE]
where is the matrix specifying the optimization waveform . In order to enforce for all , we constrain for . The control sequence that is implemented experimentally is specified by another matrix the elements of which are calculated as
[TABLE]
and
[TABLE]
after a suitable is found.
Our objective is to find a control sequence that would yield as well as
[TABLE]
for all . We assign equal weights to each member of the ensemble and define the following combined target function
[TABLE]
We construct a set of Van Loan generators that decompose into a direct sum of and :
[TABLE]
the corresponding Van Loan propagators of which are given as
[TABLE]
Here, we will slightly abuse our sub-matrix index notation, just as it was done in Section 5.3, and write the target function in Equation (129) as a function of :
[TABLE]
We now evaluate the partial derivatives of Equation (137) with respect to the elements of for each :
[TABLE]
In order to carry out the gradient ascent searches to find that yields , we evaluate partial derivatives of with respect to the elements of :
[TABLE]
and
[TABLE]
Our control searches are conducted in the way it was described as it was described at the beginning of the section. The three individual quantities that appear in Equation (137) are the unitary metric , dipolar metric , and metric , which are defined as
[TABLE]
Our representative set for the Rabi strengths is , , , , , , , , , , , , , . The parameters for the control search were , , and . The control sequence and its figures of merit are shown in Figures 7 and 8. The sequence length of corresponds to Rabi cycles for the spins experiencing .
6 Conclusions
We have developed a general method for computing perturbation integrals arising in numerous effective Hamiltonian control schemes, which are generally of the form:
[TABLE]
where the system propagator , , is generated by piecewise constant control sequences . Our method is based on the observation that these integrals can be written as solutions to first order matrix differential equations, which we call the Van Loan equations, that have the same general form as the Schrödinger equation. Consequently, the same optimization algorithms that have been developed for quantum control may be applied to optimizing these expressions, and importantly, this formulation enables the immediate application of methods that also compute and make use of gradient information (e.g. GRAPE [65] or GOAT [46]).
Gradient evaluation for Dyson terms is crucial for ensuring fast convergence to high accuracy solutions, as gradient information increases the efficiency of the search. Here, we have only used first order derivatives, but as is described in [66], incorporating second order derivative information into standard quantum control problems can ensure fast convergence rates near the optimum, and lead to better performance than first order methods [45]. To reiterate, as we rephrase the effective Hamiltonian control problem into a bilinear control theory problem, all of the higher order algorithms and methods outlined in [45] may be applied. Consequently, we expect the methods presented here to be useful tools for efficiently exploring the space of control sequences, to find solutions satisfying large number of design criteria involving Dyson terms.
As a computational task, this method has the benefits of both implementation efficiency, as well as computationally efficiency. With regards to implementation, given any software package for solving general bilinear control theory problems, the only new ingredient necessary to optimize integrals of the type in Equation (144) is the construction of the relevant Van Loan differential equations.444Indeed, the reason our implementation of GRAPE described in Appendix A uses the commutator method for gradients, as opposed to the exponential method, as in [45], is exactly due to the utilization of an already-existing code base for bilinear control problems. In terms of efficiency, Van Loan [42] found — in the context of constant generators — that the number of operations in this approach is favourable as compared to other methods (requiring either fewer, or a comparable number of operations). In particular, for a given accuracy, the method is much faster than using numerical integral rules. This has been corroborated by numerical experiments in [67], and has been observed in [44], in which the method significantly outperformed numerical integration methods for computing spin relaxation theory expressions, enabling simulation of larger systems.
In this manuscript, we have also demonstrated the use of this method in a variety of scenarios, including the successful application in nanoscale magnetic resonance experiments [50]. We note that the examples presented in this manuscript should not be understood as comprehensive; in order to keep our treatment illustrative and concise, we did not present control searches involving Dyson terms higher than first order. However, such extensions follow naturally from the given examples, and in our experience optimizations with higher order terms have also converged to desired solutions. Moreover, we wish to highlight that these block matrix methods are applicable to terms of even more general form than those appearing in perturbative control schemes, as are captured in Equation (144). Any nested integral expressible as an output to Theorem 1 could be incorporated into control optimization.555As a concrete example, a unitary trajectory forms a continuous-time unitary -design if and only if the integral takes a specific value. This integral may be computed, and hence optimized, using the methods presented here.
The methods developed here could, in principle, be used in any setting requiring robust coherent control of quantum systems, whether it be quantum computation, sensing, or spectroscopy. The main requirement for successful implementation is an accurate model of the system generators, and a precise knowledge of the control amplitudes seen by the quantum system, i.e. sufficiently good characterization of the experimental transfer function. Quantum control problems that could be addressed with these tools in the future include, but are not limited to: minimizing the effect of cross-talk in the case of simultaneous multiple qubit control, reducing the effect of counter-rotating terms in the cases where the Rabi strength approaches the qubit level spacing (i.e. Bloch-Siegert type effects), and preventing leakage to higher levels. Another potential application of this framework is the inclusion of non-Hamiltonian Lindblad terms for cross-polarization problems.
Acknowledgements
This work was undertaken thanks in part to funding from the the Canada First Research Excellence Fund (CFREF). Further support was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Excellence Research Chairs (CERC) program (215284), the Canadian Institute for Advanced Research (CIFAR), the province of Ontario and Industry Canada. H. H. and D.P. would like to thank Ian Hincks, Raffi Budakian and William Rose for their insights and fruitful discussions. F. Z. would like to acknowledge Guilu Long for his insights and inspiring discussions.
Appendix A GRAPE Algorithm
In this appendix, we describe the Gradient Ascent Pulse Engineering (GRAPE) algorithm introduced by Khaneja et al [65] for numerically solving bilinear control theory problems for piecewise constant control amplitudes. A control problem over a time interval is bilinear if the control amplitudes determine the system evolution through a first order matrix differential equation:
[TABLE]
where for all and . It is clear that the Schrödinger equation is a special instance of bilinear control with being anti-Hermitian matrices, as are Van Loan equations.
Throughout this manuscript, we only deal with piecewise constant control amplitudes , for which, we split the interval into subintervals with respective durations such that , for all , and . As we said in Section 3 we store the piecewise constant values of in a real valued matrix , the elements of which are determined by Equation (25). As such,
[TABLE]
where the product symbol denotes a sequential matrix multiplication.
A bilinear control theory problem for piecewise constant control amplitudes is stated as: find control amplitudes , or equivalently the corresponding matrix , that adhere to certain problem specific constraints and yields a with some desired properties. Given the desired properties for , we can always write down a target function that is an analytic function and takes the value if and only if its argument has the properties that we want from . Having defined such a target function, the problem of finding a that yields becomes a multivariable optimization problem.
Because is an analytic function of , knowing and its partial derivatives for and , is enough for determining both and . The key insight by Khaneja et al [65] was to point out that the computational cost of the simultaneous evaluation of and is not much more than evaluating alone. As such, it is advantageous to use optimization algorithms that make use of the gradient information, e.g. the conjugate gradient algorithm.
Here, we proceed by evaluating a partial derivative
[TABLE]
with
[TABLE]
The integral term in Equation (148) can be evaluated either through block matrix techniques introduced in Section 4 with Equation (29) or by noticing that
[TABLE]
and approximating the integral with a finite sum of commutators. During the optimization of the examples presented in this manuscript, we always approximate the integral with a finite sum of 15 commutators since the speed of the optimization was not our main concern in this work. We note that a very comprehensive analysis of the performance of various algorithms for bilinear control optimization is given in [68], which also investigates the use of block matrix methods for evaluating the integral expression in Equation (149).
Our scheme for control searches is much the same as the one by Khaneja et al [65]. We always define target functions that are analytic functions of , i.e., , . The partial derivatives of with respect to are then evaluated in terms of using the chain rule. We use an off-the-shelf gradient ascent optimizer to maximize starting from a set of initial control sequences until the algorithm yields a value sufficiently close to one.
Appendix B Transfer Functions
In this appendix, we describe the matrix methods used for performing the control searches for two examples in Section 5. In these cases, the control sequences , , are piecewise constant over intervals of equal length such that . Furthermore, in these cases, the control vectors are of dimension two, i.e., . Consequently, we map the real valued control vectors onto a complex scalar function , with . Given the complex vector representation and the fact that all transfer functions in this work are linear functions that treat and symmetrically, we represent each as an matrix with complex entries such that
[TABLE]
and
[TABLE]
where are the piecewise constant control amplitudes as defined in Section 3 and denotes the th row of the matrix , that specifies the control sequence .
Moreover, for all examples in this manuscript are diagonal in the Fourier domain and fully specified by two real valued scalar functions – the amplitude transfer function and the phase transfer function . For the example in Section 5.5 the were constructed from and given in Figure 6(a) and 6(b), respectively. In order to construct each , we first calculate the discrete Fourier transform matrix, which is a unitary transformation , with its elements given as
[TABLE]
for . We then construct a diagonal matrix , the diagonal elements of which are given as , with
[TABLE]
for . Here, denotes the remainder from diving by . Finally, we evaluate
[TABLE]
for each . The control amplitudes are then evaluated via Equation (150) and (151). The elements of the Jacobian , that are also necessary for the control searches, evaluate to
[TABLE]
for all .
B.1 Optimization Transfer Functions
In Section 3, we argue that to implement certain constraints on the control sequence , one can use an optimization transfer function such that , where and specify the piecewise constant control sequences and , respectively. is constructed in such a way that its output functions always adhere to the necessary constraints. In this subsection, we demonstrate explicitly how to construct that introduces periods of zero pulse amplitudes at the beginning and at the end of the control sequence and limits the bandwidth of in the Fourier domain. This optimization transfer function is used for two control searches in Section 5.
First, we construct an optimization transfer function that ensures that the control sequence has equal periods of zero amplitude at the beginning and at the end of the sequence. To introduce such periods, we define an optimization transfer function . Just like above, we map both and onto complex vectors and , respectively. is the number of piecewise constant elements of , where is the number of zero amplitude elements that are introduced at the beginning and at the end of . The action of is then implicitly defined as
[TABLE]
for any ; denotes a zero vector of length . It is easy to see that a , which has the above property, can be constructed from three blocks:
[TABLE]
where is a rectangular matrix with all its entries being zero. The elements of are then given by
[TABLE]
and
[TABLE]
Of course, the construction of generalizes easily for introducing an arbitrary number of zero amplitude periods of an arbitrary duration into .
Now, to limit the bandwidth of frequency components, we construct an optimization transfer function that acts as a low-pass filter. Here, a low-pass filter should be understood simply as some amplitude transfer function in Equation (154), that takes non-zero values only over some limited range centred around . In this manuscript, we used a particular functional form
[TABLE]
that has smooth frequency cut-offs at in order to prevent introducing long lasting ripples to the pulse waveform . Accordingly, we implement as
[TABLE]
where is given by Equation (154).
The optimization transfer function that we used for the two examples in Section 5 combined the action of and and is calculated as
[TABLE]
The elements of are consequently determined as
[TABLE]
for all , whereas the elements of the Jacobian evaluate to
[TABLE]
for all and .
Appendix C Symbolic Computation Methods
In this appendix, we describe the design and usage of code in Mathematica for symbolically simplifying the structure of the matrices of Theorem 1, with the main goal being to verify Conjecture 1. More generally, automated simplification of the integral expressions arising in Theorem 1 may be useful, as many of the generators naturally constructed in the quantum control context are sparse with repetitive entries, and as such the integral expressions will be amenable to significant simplifications.
The code described here is implemented using Patterns, Pattern Matching, and Replacements in Mathematica. In particular, we use symbolic expressions in Mathematica to represent integrals arising from Theorem 1. Expressions that can be simplified are identified using Patterns, which are used in Mathematica to identify expressions with a particular structure. If a Pattern matches an expression, then we apply Replacements, which transform one expression into another. For example, if one of the matrices in the symbolic expression for an integral is [math], we use a Pattern to identify this, and a Replacement to replace the expression with [math]. See [69] for an introduction to these concepts in the Mathematica programming language.
This appendix is organized as follows:
- •
In Section C.1, we outline a representation for the integrals in Theorem 1 using symbolic expressions in Mathematica.
- •
In Section C.2, we define replacement rules to be used for automatically simplifying the integral expressions.
- •
In Section C.3, we describe the implementation of a function that implements Theorem 1; i.e., given a list of the matrices which are the input to Theorem 1, the function constructs and simplifies the matrices, which are the outputs of Theorem 1.
- •
Finally, in Section C.4, we apply the functionality built in the preceding sections to verify Conjecture 1 for .
The Mathematica notebook containing the code developed in this appendix may be found in the online repository [70].
C.1 Reserved Symbols and Primitive Expressions
The first step is specifying a representation for expressions appearing in Theorem 1 in Mathematica. We reserve the symbol “” to represent integration variables, and the symbol to represent the identity matrix, which can be produced in Mathematica by typing Esc d s 1 Esc.
There are two types of expressions which can appear in simplifications of Theorem 1 that we consider. In this appendix, we call these primitive expressions, and use the term to refer to both the mathematical objects, as well as their representation in the code. The first kind of primitive expression is the nested integral (note the additional appearance of the pre-factor ):
[TABLE]
which is represented in the code using the head “Int”:
[TABLE]
The second primitive is simply an expression not appearing in an integral, for example , which is represented in the code using the head ”Ex”:
[TABLE]
We remark that the programmatic representation of nested integrals is slightly redundant for representing expressions resulting from Theorem 1, as primitives from this theorem will always have , but there is no harm in this redundancy.
A final technical detail is a special representation for time-independent scalar multiples of matrices, i.e. , for and a matrix valued function. These are represented in the code in the following way:
[TABLE]
which eases the discrimination between symbols representing (potentially time-dependent) matrices and time-independent scalars.
C.2 Simplification of Primitives via Replacement Rules
The simplification of primitive expressions is carried out by defining replacement rules, then applying those rules to primitives using the internal functions of Mathematica. The replacement rules we specify correspond to basic simplifications of the nested integrals arising in Theorem 1. The rules are broken into groups.
- •
Zeroes
The detection of zero matrices:
[TABLE]
- •
Linearity of integration
Factoring scalars:
[TABLE]
and linear combinations:
[TABLE]
- •
Performing integrals
In some cases, an integral can be explicitly performed. The case that we handle is when and for some index in Equation (171). As an explicit example, it holds that
[TABLE]
Handling all possible versions of this simplification with replacement rules must be broken into four cases based on the order of the nested integral, as well as where the simplification appears in the nest. Some care must also be taken with powers of ; the pattern will detect powers of when , but not the expressions “” and “” as, symbolically, they are not powers of . Thus, one must create replacement rules for these cases separately.
- –
Case 1: A first order integral that can be performed.
[TABLE]
- –
Case 2: An integral of order where the last integral can be performed.
[TABLE]
- –
Case 3: An integral of order where the first integral can be performed.
[TABLE]
- –
Case 4: An integral of order where an “internal” integral can be performed.
[TABLE]
C.3 Computing General Expressions from Theorem 1
Computing general expressions arising from Theorem 1 consists of programmatically implementing the mapping
[TABLE]
where the right-hand side is the time-ordered exponential of the left-hand side. The function TOExponential implements this mapping. As
[TABLE]
(i.e. depends only on ) rather than specifying , the user specifies the symbols for the matrices as an input. That is, TOExponential takes in two lists of symbols:
[TABLE]
where the first list of symbols represents the diagonal blocks of the time-ordered exponential, and the second list represents the off-diagonal blocks of the matrix to be time-ordered exponentiated. The output is a list of symbols
[TABLE]
representing the upper triangle of matrices in the right-hand side of Equation (183). In the notation of Theorem 1, we have that the expression for is contained in .
The Mathematica notebook has several example applications of this code. One basic use-case is to show that:
[TABLE]
for and , where . In this case, the inputs to TOExponential are the lists
[TABLE]
and the output is the list
[TABLE]
By interpreting the structure of this output according to the data representation for primitives in Appendix C.1, we see that this output is correct.
C.3.1 Implementation of TOExponential
The implementation of TOExponential requires two pieces: the symbolic construction of expressions for the matrices given in Theorem 1, and the simplification of these expressions via the replacement rules of Appendix C.2. Due to the recursive structure of the matrices in Theorem 1, it is both conceptually natural and computationally more efficient to perform this computation in a recursive fashion. The recursion proceeds by computing successive off-diagonals of the matrices (i.e. successive values of for all ). As a reminder, the recursion relation for the matrices for is:
[TABLE]
Roughly, the computation goes as follows:
Base case: Initialize the case using Equation (190) and apply replacement rules to simplify any expressions. 2. 2.
Recursive step: Construct expressions for using Equation (190) and the already computed matrices for . Apply replacement rules to simplify resulting expressions.
The most basic recursive construction step is to produce a single term in the sum in Equation (190). That is, given (a linear combination of primitives) and two symbols and , we must programmatically produce the mapping
[TABLE]
This mapping is performed by two functions both named RecursiveRule, where the second handles the case when the symbol is given as some scalar multiple . The replacement rules are
[TABLE]
and
[TABLE]
and for the version which handles scalar multiples they are
[TABLE]
and
[TABLE]
A full recursion step (for computing the next off-diagonal of elements from the previous) is implemented by RecursiveConstruct, which simply uses the RecursiveRule functions to construct the whole of the right-hand-side of Equation (190).
Lastly, the implementation of TOExponential is to do the initialization step of constructing and simplifying the matrices, then recursively calling RecursiveConstruct to populate the output consisting of all matrices.
C.4 Bivariate Polynomials and Conjecture 1
Here, we describe code for working with and analyzing the generators for the purposes of verifying Conjecture 1. A description of is in Section 4.3.2, in the lead up to the statement of Conjecture 1.
The first step is simply to produce the generator and its time-ordered exponential. The functions BPODGenerator and BPGenerator both serve the function of specifying , with the first producing the upper-off-diagonal pieces of the generator, and the second including the diagonal. The function BPTOExponential simply applies the function TOExponential of Appendix C.3 to the generator given in BPGenerator and BPODGenerator.
The problems of verifying Conjecture 1, and providing explicit linear combinations of the blocks of for computing integrals of the form
[TABLE]
for polynomials of degree at most ), are fundamentally problems of linear algebra, and hence it is necessary to represent the upper right blocks of in a way that they may be analyzed using the linear algebra functions in Mathematica.
C.4.1 Linear Algebraic Representation
Due to the linearity of integration, for a polynomial of degree , the expression in Equation (196) may be viewed as member of the vector space of expressions
[TABLE]
where, again
[TABLE]
and each is considered to be linearly independent from all others. Denoting this vector space of expressions as , we call the the elementary basis for this vector space.
We may represent as column vectors by defining a linear mapping
[TABLE]
which acts on the standard basis as
[TABLE]
where is the column vector with a in the position and a [math] everywhere else (i.e. sends the standard basis of to the standard basis of ). Note that ordering of the images of the standard basis elements of according to this mapping is the lexicographic ordering of the basis elements according to the powers of and appearing in the integral.
The function BPVectorRep constructs a list of replacement rules that corresponding to this mapping. The function BPTopRightBlockVectors produces a matrix whose column vectors correspond to the top-right blocks of under the above mapping, with the blocks ordered in terms of the lexicographic ordering of their indices (which is automatically implemented by the Flatten function in Mathematica).
C.4.2 Testing Conjecture 1
With the terminology of the previous section, the content of Conjecture 1 is simply that the top right blocks of the time-ordered exponential of are a basis for the formal vector space . This is equivalent to the columns of the matrix output by BPTopRightBlockVectors being a basis for . The function TestClaim1 checks if this is the case by computing the rank of the matrix , with the claim being true if and only if the rank is . Note that before checking this, the code checks if any primitive expressions are present (i.e. whether the symbols “Ex” or “Int” remain, which would mean that there are expressions in the top right blocks appearing that are not expected). If any such expressions are found, the function outputs the value indicating that Conjecture 1 has failed catastrophically. If no such expressions are found, the rank is checked, with an output of meaning Conjecture 1 is verified, and [math] meaning it is invalidated. For example, the matrix in the case is
[TABLE]
which is clearly full rank.
The function TestClaim1Range applies TestClaim1 to a range of values of and . We have tested this claim for and have found it to be true in all cases. While this does not constitute a proof, it does provide confidence that it is true in general, and in practical settings one can simply check its validity for a specific and of interest.
C.4.3 Computing Integrals with Bivariate Polynomials
Given an arbitrary polynomial we wish to write the integral
[TABLE]
as a linear combination of the top-right blocks of . Assuming the truth of Conjecture 1, the top-right blocks of this time-ordered exponential are a basis for , which we will call the exponential basis. Let be the mapping which takes an element of and returns the column vector of coefficients corresponding to expanding it according to the exponential basis (in lexicographic ordering).
The matrix output from BPTopRightBlockVectors then simply represents a change of basis from the image of the exponential basis under to the image of the standard basis under . That is, for every , it holds that
[TABLE]
Hence, the vector of coefficients of representing its expansion in the exponential basis is
[TABLE]
As an example, in the case, for an arbitrary corresponding to the polynomial
[TABLE]
we have
[TABLE]
and hence
[TABLE]
Denoting the blocks of by , we may write the desired integral explicitly as a linear combination of matrices by taking the following dot product
[TABLE]
where the vector of matrices is just the top blocks of in lexicographic ordering.
The function PolynomialBlockDecomposition outputs the right hand side of Equation (224) for an arbitrary and by performing the above computation in generality. An example when is included in the code.
Appendix D Proof of Theorem 1
Recall, the goal is to prove that the matrices, defined implicitly by the equation
[TABLE]
satisfy
[TABLE]
and
[TABLE]
for all and . We also show that alternatively, for , they can be given recursively:
[TABLE]
The approach of the proof is to first show that Equations (226) and (228) hold through an application of the uniqueness of solutions to differential equations. Afterwards, we show that the expression in Equation (227) satisfies the same recursion relation as is given in Equation (228), and therefore Equation (227) is also correct.
To prove that Equations (226) and (228) hold, we first use the definition of the time-ordered exponential to obtain the explicit form of the differential equation that the satisfy. Differentiating both sides of Equation (LABEL:equation:appendix_exponential), we see that:
[TABLE]
Hence, for and , the matrices satisfy the differential equation
[TABLE]
with initial conditions
[TABLE]
where we are using a dot to denote differentiation with respect to . From here, the validity of Equation (226) follows from the case of Equations (230) and (231). To establish the recursion relation in Equation (228), we differentiate it:
[TABLE]
where in the last equality we have applied the recursion relation. Hence, the expressions in the recursion relation satisfy the same differential equations as the case of Equations (230) and (231), and hence they are correct under the assumption that this differential equation has a unique solution.
To show that Equation (227) is also correct, we show that it satisfies the same recursion relations given by Equation (228). As the first term in Equation (227) is immediately equal to the term in the sum of Equation (228), and the remaining terms in both are only non-zero when , it is sufficient to establish that
[TABLE]
when . First, expand the left-hand side of the above equation:
[TABLE]
where we have broken the sum into the and case. We may rewrite the first summation as:
[TABLE]
where in the second equality we have applied the recursion relation for Int in Equation (35). We may similarly rewrite the terms in Equation (234) using a series of steps given below, with a description of each step given after the chain of equalities:
[TABLE]
The operation in each equality given by:
- (1)
Break up the inner sum, with the new index having the correspondence . 2. (2)
Swap the order of summation over and . 3. (3)
Apply the recursion relation for Int given in Equation (35). 4. (4)
Move the summation for and past all terms that have no dependence on these indices. 5. (5)
Change the limits for summation over to start at . 6. (6)
Substitute the expression in the brackets using Equation (227). 7. (7)
Increase the upper limit of summation over to include , which does not change the sum as .
To complete the proof, add the new forms for the and terms in Equation (234), given respectively in Equations (235) and (236), to conclude that Equation (233) holds.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] L. M. K. Vandersypen and I. L. Chuang. NMR techniques for quantum control and computation. Reviews of Modern Physics , 76:1037, 2005.
- 2[2] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver. A quantum engineer’s guide to superconducting qubits. Applied Physics Reviews , 6:021318, 2019.
- 3[3] C. D. Bruzewicz, J. Chiaverini, R. Mc Connell, and J. M. Sage. Trapped-ion quantum computing: Progress and challenges. Applied Physics Reviews , 6:021314, 2019.
- 4[4] L. Viola and S. Lloyd. Dynamical suppression of decoherence in two-state quantum systems. Physical Review A , 58:2733, 1998.
- 5[5] L. Viola, S. Lloyd, and E. Knill. Universal Control of Decoupled Quantum Systems. Physical Review Letters , 83:4888, 1999.
- 6[6] L. Viola and E. Knill. Robust Dynamical Decoupling of Quantum Systems with Bounded Controls. Physical Review Letters , 90:037901, 2003.
- 7[7] E. Lucero, J. Kelly, R. C. Bialczak, M. Lenander, M. Mariantoni, M. Neeley, A. D. O’Connell, D. Sank, H. Wang, M. Weides, J. Wenner, T. Yamamoto, A. N. Cleland, and J. M. Martinis. Reduced phase error through optimized control of a superconducting qubit. Physical Review A , 82:042339, 2010.
- 8[8] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm. Simple Pulses for Elimination of Leakage in Weakly Nonlinear Qubits. Physical Review Letters , 103:110501, 2009.