Second order phase dispersion by optimized rotation pulses

We show that the duration of broadband universal control pulses can be halved by choosing control targets with a quadratic function of phase dispersion. This class of control pulses perform a broadband universal rotation around an axis, in the Bloch sphere representation of two-level systems, given by this phase dispersion function. We present an effective optimal control method to avoid the problem of convergence to local extrema traps. DOI: 10.1103/PhysRevResearch.2.033157


I. INTRODUCTION
The critically important milestones on the road-map of quantum technology development are expected to include quantum optimal control [1,2].The main result of this communication is to reduce the time and energy required to control a quantum system.As an application of optimal control theory, quantum optimal control can be tasked with finding control pulses, usually in the form of electromagnetic radiation, to perform a desired operation on an arbitrary quantum state.Already this method has been applied successfully to a wide range of experiments, including magnetic resonance spectroscopy [3,4], error-correction for quantum computing [5], quantum information registers [6], robust atom interferometry [7,8], and high-resolution medical imaging satisfying legal irradiation constraints [9].
The growing number of reported applications of quantum optimal control are expanded by the concurrent development of optimal control algorithms.There are a number of approaches to optimal control, such as Lagrangian methods [10], time optimal control [11], annealing quantum systems to a desired effective Hamiltonian [4], sophisticated gradientfree searches [12], rapidly converging Krylov-Newton methods [13], and optimal control using analytic controls [14].This development is essential as desired solutions push to the limit of what is physically possible and control problems become computationally and numerically challenging [15].
Quantum optimal control problems can be solved numerically by using a piecewise-constant control pulse approximation [16][17][18][19], and the trajectory gradient can be calculated and followed by using the gradient ascent pulse engineering * david.goodwin@chem.ox.ac.ukPublished by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
(GRAPE) method [3].The initial GRAPE method showed linear convergence to a desired optimum destination and was later extended to give superlinear convergence [20], then quadratic convergence [21,22].This communication introduces an ordered optimal control method, based on the GRAPE method, to avoid local convergence traps encountered in the optimal control problem set out below.This method is termed morphic-GRAPE by the authors.
In addition, and vital to effective experimental applications, optimal control has advanced theoretical studies calculating universal gates for quantum computing [23][24][25][26][27].A subset of universal quantum gates include control targets designed as an effective rotation in the Bloch sphere representation of a two-level system, with a defined rotation axis and rotation angle, for any initial basis-state vector [6,28].The developments reported here will be concerned with this universal rotation class of control pulses.
One of the drawbacks associated with universal rotation solutions is the long duration of control pulses when compared with the easier control problem of optimizing state-to-state problems [1,28].The optimal control method presented here will show that the previously required pulse duration can be halved by defining target propagators as a function of phase dispersion.

II. BROADBAND OPTIMAL CONTROL TO EFFECTIVE PROPAGATORS
An ensemble of noninteracting two-level systems should include frequency dispersion terms from local environmental conditions.Uniform manipulation of the ensemble is a difficult practical task but by describing the ensemble as a bilinear control problem [16,18], numerical optimization can be used to find solutions that are difficult to find analytically.In formulating this bilinear control problem, the controllable parts of the Hamiltonian are the pulses, and the uncontrollable part is the local frequency dispersion term.Restricting control pulses to phase modulation, ϕ(t ), with constant amplitude A, the time-dependent Hamiltonian for this ensemble of K noninteracting two-level systems can be written as where the angular frequency ω k is variously termed the resonant frequency offset, chemical shift, or detuning and describes the local frequency dispersion within the ensemble.Practically, this ensemble of resonant frequency offsets form a discrete grid of noninteracting two-level systems spread over a relevant bandwidth [29,30].σ (k) x,y,z are operators of the kth two-level system, operating in a Hilbert space and related to the Pauli matrices σx,y,z by where E (k) is a K × K single-entry matrix with the diagonal element E (k) kk = 1.The GRAPE (gradient ascent pulse engineering) method of optimal control [3] proceeds to describe the time-dependent control pulses as piecewise constant over a small time interval t [18,19].This approximation allows for a numerical solution of Eq. ( 1) through time-ordered propagation, with the Liouvillian Ln being an adjoint representation of a Hamiltonian: L = 1 ⊗ Ĥ − Ĥ † ⊗ 1 [20].At a time increment n, the effective propagator U n evolves the system forward in time over the interval [0, t n ], and the effective propagator of the adjoint control problem, V n , evolves backwards from a desired target R over the interval [T, t n ].The effective propagator U N is the solution to the ensemble Hamiltonian (1) over the total pulse duration T .The effective propagator of the adjoint problem, V n , is needed for a fidelity gradient calculation, outlined below.The Hilbert-Schmidt inner product of the desired effective propagator R and the effective propagator over the control duration, U N , gives a measure to numerically maximize and is termed the fidelity [3,4,6,28]: This fidelity measure F is normalized by the dimension of the Hilbert space, d, to give sensible and predictable bounds F ∈ [−1, +1].The desired target R is interpreted as a rotation in the Bloch sphere representation in all that follows.Given an initial guess for control pulses, optimal control pulses can be found by stepping in the direction of the fidelity gradient vector over control manifold, until a maximum fidelity is found.The distance to step is set by an appropriate Newton-type optimization method [31].The gradient-following GRAPE method [3] requires a fidelity gradient vector ∇F and, for phase modulated control, this is constructed from the elements The two partial derivatives on the right-hand side are directional propagator derivatives, and can be interpreted at the derivative of the propagator, Pn , in the direction of the operators σx and σy .The fidelity gradient is calculated by using the same fidelity definition of Eq. ( 4), with the directional propagator derivatives sandwiched by the effective propagator and its adjoint in Eq. ( 3).
For each time slice, the time propagator in Eq. ( 3) and each of the directional derivatives in Eq. ( 5) can be calculated analytically with one exponentiation of a block-triangular matrix [32,33], in each direction σ ∈ { σx , σy }.Considering the large and sparse Hamiltonian formulation of Eq. ( 1), the main computational cost is 2N exponential operations, which can be performed in parallel and can include a efficient propagator recycling scheme [22].
Pulses are linearly scalable with the desired bandwidth [1], with the transform from one pulse p 1 , with duration T 1 , amplitude A 1 , and covering a bandwidth 1 , to a second pulse p 2 , with duration T 2 , amplitude A 2 , and covering a bandwidth 2 , having the same effect if where the subscripts refer to the according property of pulse p 1 and p 2 .The dimensionless bandwidth factor b and scaling factor s are defined here to characterize the pulse, rather than the system-specific quantities of bandwidth, pulse amplitude, and pulse duration.This is common in the discipline of magnetic resonance, and the rotation angle β is included in the scaling factor to normalize with this use in magnetic resonance.In this definition, bandwidths are measured in Hertz, pulse duration in seconds, pulse amplitudes in radians per second, and rotation angles in radians.In the investigation that follows the scaling factor is chosen to be s = 2β πb which, in effect, gives the same pulse amplitude over all bandwidth factors for a defined system.
The choice of constant amplitude, phase modulated, control pulses is one of inference: both chirped pulses [34][35][36][37][38] and universal rotation pulses [28] are close to constant amplitude, and it is expected that pulses produced in this work will have a similar form.
The piecewise-constant time slices used in calculations are a very good approximation in nuclear magnetic resonance (NRM) spectroscopy.In cases where this piecewise-constant approximation is not valid, such as applications to electron paramagnetic resonance (EPR) spectroscopy, a feedback control method [39,40] or a transfer-matrix method [41][42][43] can be used to calibrate pulses due to these hardware-specific cavity effects.Further to this, optimal control can be designed to include a robustness to pulse amplitude miscalibration [29].

III. ROTATION AXES WITH QUADRATIC PHASE DISPERSION
Derived from Landau-Zener-Stückelberg-Majorana theory [37], the ensemble of Eq. ( 1) can be controlled with adiabatic evolution [44][45][46] and is described as a linear frequency sweep over the two-level systems.The phase dispersion resulting from a linear frequency sweep is quadratic and can be used to define the local targets of an optimal control method.Previous work on this theme has optimized state-to-state problems with linear phase dispersion [47][48][49] and quadratic phase dispersion [44,[50][51][52].
Although state-to-state solutions are useful, they are only functionally defined for a specified state preparation: the desired control is only effective for a given initial state of the system.A universal control method, independent of the initial state, attempts to find desired unitary propagators which represent a rotation of the Bloch sphere; rotating all components of the coordinate system about an axis through the origin [6,28].
The effective unitary propagator representing a universal rotation can be formulated for each member of the ensemble, k, for a rotation angle β.The desired axes of rotation, nk , can be defined as where the angle α k is a phase on the transverse plane of the Bloch sphere and is described by a phase dispersion function.This study uses a quadratic phase dispersion similar to that of the linear frequency swept chirped pulses [34][35][36][37][38]: This is the phase that the transverse components of states acquire during the pulsing time [47].The introduction of Q is defined here as the quadratic coefficient of a phase dispersion function, which is closely related to previous work using a linear phase dispersion function [47].The case when Q = 0 has been previously published as broadband universal rotation by optimized pulses (BURBOP) [28].
It should be emphasized that the desired pulses in the communication are an improvement towards universal rotation pulses with reduced control duration, using a phase dispersion function similar to chirped pulses, rather than a generalization of chirped pulses to universal rotation pulses.Furthermore, the function in Eq. ( 10) is a quadratic function of Q but could be formulated, in principle, by any function of Q.This new class of pulses are named second order phase dispersion by optimized rotation (SORDOR) pulses by the authors.

IV. MORPHIC OPTIMAL CONTROL
In the authors' experience, the standard procedure of starting GRAPE from a random guess or an informed guess, e.g., an adiabatic passage pulse, produces a fidelity profile periodically reducing to zero over the bandwidth, particularly when setting Q 0.1 in Eq. (10).The reason for this is that, as the optimization progresses, fidelity sections either side of the zeroed fidelity increase, also increasing the average fidelity, trapping these local fidelity minima into the control manifold [15,22].This can make the process of finding the desired optimal solutions set out in Eq. ( 8) a difficult task.Furthermore, a high average fidelity over the bandwidth is not sufficient for acceptable pulse performance-the pulse should uniformly manipulate the entire bandwidth.
An ordered optimization procedure of morphic-GRAPE is set out below to avoid these traps.This method is effective, yet easy to use, and proceeds by morphing one optimization problem into an incrementally modified optimization problem; using the solution of one as the new starting point of the next.
A grid of coefficients, Q ∈ [0, 1], and bandwidth factors, b ∈ (0, 18], sets a morph area for the results that follow.SORDOR pulses are created with four directional morphs: forward morph, incrementing Q with + Q; backward morph, decrementing Q with − Q; compressed morph, interpolating a control pulse from b to b − b; and expanded morph, interpolating a control pulse from b to b + b.An additional smoothing stage can be used to remove discontinuities in the fidelity profile, using sequences of forward or backward morphs from the largest extrema in the dF dQ profile.The ordered recipe of finding SORDOR pulses is outlined as follows: (1) Start with a universal rotation pulse (Q = 0) and use morphic-GRAPE: ).The smoothing morph starts from the eleven largest extrema, with the first five of these shown in Fig. 1. Figure 1(a) shows the places on the Qcoefficient grid where the smoothing starts, after the forward and backward stages, and Fig. 1(b) shows how the fidelity increases and decreases over these stages.Both are plotted as a function of the number of gradient calculations which shows the computational cost of morphic-GRAPE at b = 18.0, which can be regarded as the seed for the compressed and expanded stages.The forward and backward stages at different values of b show a similar profile to those stages in Fig. 1.
A direct comparison of the computational cost of this morphic-GRAPE method to the previously published universal rotation pulses [28] is not appropriate, because morphic-GRAPE essentially includes that previous work.

V. RESULTS
The main results of this work are shown in Fig. 2. The SPINACH optimal control toolbox [20,22,33] is used to simulate the ensemble of two-level systems, defined by Eq. ( 1), with minimal modification to enable optimization to the desired unitary propagators in Eq. (8).A linearly spaced resonance offset profile in Eq. ( 1), with K = 1 + 10b , and a piecewise constant approximation of Eq. ( 3), with N = 50b, is optimized with the -BFGS method [20], forming a Hessian approximation from the previous twenty gradients.
It is not only interesting to see how the fidelity increases over the grid of Q coefficients but it is also practically useful because experiments requiring combinations of pulses with different rotation angles β should be matched, i.e., a sequence consisting of pulse p 1 followed by pulse p 2 , should have b p 1 = mb p 2 with Q p 2 = mQ p 1 , where m is an arbitrary multiplier.The maximum fidelities found over the grid of Q coefficients and bandwidth factors, from the forward, backward, compressed, and expanded morphic optimizations, are shown in Figs.2(a) and 2(b).
It should be noted that the morphic optimizations for β = π 2 are difficult, numerically and computationally; the backward, compressed, and expanded stages of the morphic optimizations give little improvement when β = π , and were designed primarily for the case when β = π 2 .Figure 2(c) shows the maximum fidelity at each bandwidth factor compared with the fidelity of an equivalent universal rotation pulse, broadband universal rotations by optimized pulses (BURBOP): the shortest duration universal rotation pulses known to the authors [1,28].The BUR-BOP in this comparison have the same constant amplitude A and scaling factor s as the SORDOR pulses.BUR-BOP pulses start from the same ones as published previously [28], then are further optimized with the same exact gradients [Eq.( 6)] and optimization tolerances.It is clear that SORDOR pulses show significant improvement over BURBOP.

VI. EXAMPLE OF AN APPLICATION IN MAGNETIC RESONANCE SPECTROSCOPY
This section will proceed to apply the phase modulated pulses above, the subject of this article, to the experimentally relevant system of nuclear magnetic resonance spectroscopy.
As an example, a spin-1 2 13 C nucleus in a static magnetic field of B 0 = 14.0954T forms the ensemble of two-level systems in Eq. (1).For the analysis of small molecules, pulses should cover a bandwidth of = 37.5 kHz [53] at this magnetic field due to chemical shielding effects from varying molecular environments.This chemical-shift range is linearly spaced with ω k ∈ [−π , +π ] covering a bandwidth of = 40 kHz, slightly more than the required 37.5 kHz, and represents a 13 C spin ensemble with K = 451 members.
Pulses produced with the method outlined are morphed from one grid point to another in Fig. 2, and this morphing is key to obtaining high fidelities.This indicates that pulses should be related to each other and asks the question of whether pulses close on the grid look like each other.Twelve pulses produced by the expanded stage of the morphic optimal control method are shown in Fig. 3, six of each have β = π 2 in Fig. 3(a) and β = π in Fig. 3(c), at (200, 250, 300, 350, 400, 450) μs.
Interestingly, SORDOR π 2 and SORDOR π both are of very similar constitution.As shown in Figs.3(a) and 3(c), pulse shapes in all cases are symmetric in time, while no symmetry constraints are imposed during the optimization procedure.In essence, a parabolic time course of the pulse phase is observed, as is reminiscent of linearly frequency swept chirped pulses [34][35][36]54,55], but with regular spike-like features that also frequently come with approximate π , 2π , or 4π jumps in phase.The features are even better visualized if the quadratic phase contribution is subtracted [Figs.3(b) and 3(d)].Similar phase behavior has previously been reported for BIP [56] and BIBOP [57,58] β = π inversion pulse shapes, but not for β = π 2 rotation or excitation pulses.
The comparison of the SORDOR pulses to a chirped pulse show an additional phase modulation, where the chirped phase is from the formula The bandwidth of this chirped phase, c , is chosen to give the same pulse amplitude as the SORDOR pulses, A = (2π )10 kHz.Clearly, this does not match the quadratic component of all SORDOR pulses in Figs.computation time because only half the time propagators and their directional derivatives need to be calculated with Eq. ( 6).Two SORDOR pulses are further investigated in combination: a π 2 pulse from Fig. 2(a), and a π pulse from Fig. 2(b), each with Q = 0.78 and b = 18.This corresponds to a pulse duration of T = 450 μs, a time increment of t = 0.5 μs in Eq. ( 3), and the constant pulse amplitude of A = (2π )10 kHz in Eq. ( 1).Both of the SORDOR pulses were obtained from the expanding stage of morphic optimal control method outlined above, these being the best performing combination of pulses produced in Fig. 2.
Illustrative examples of sequences using these two SOR-DOR pulse are shown in Fig. 4. Figures 4(a

VII. CONCLUSIONS
Broadband universal rotation pulses with second order phase dispersion define a new class of pulses, named SORDOR pulses by the authors, with the universality of rotation pulses and an evolution character of adiabatic pulses.The difficulty in producing universal rotation pulses with second order phase dispersion, similar to adiabatic passage pulses, poses an interesting question regarding methods of optimal control.
A simple but effective morphic-GRAPE optimal control method has been designed to avoid local extrema traps and find shorter-duration control pulses when compared with the current best universal rotation (BURBOP) pulses.The trend of this comparison shows that SORDOR pulses are approximately 50% of BURBOP durations: a significant improvement.
The advance in pulse performance, compared with BUR-BOP, derives from a phase dispersion function used to define target unitary propagators.The particular phase dispersion function used in this work consists of a constant term, a quadratic term, and a variable multiplier.The exact form of this function is not derived within this work and is used in this form simply because the function produces highfidelity pulses.Studies using other phase dispersion functions, e.g., of higher or arbitrary order, are left to future investigations.

FIG. 1 .
FIG. 1.An indication of the computational cost of calculating SORDOR pulses at b = 18.0.Panel (a) shows the number of fidelity gradient calculations as the Q coefficient is increased and decreased, with the forward, backward and smoothing stages of this morphic optimal control method.Panel (b) shows the corresponding fidelity as a function of gradient calculations.Horizontal arrows in panel (a) show the point from which the waveform is morphed to start the smoothing stage.Dashed lines in panel (b) show the maximum fidelity obtained.

FIG. 2 .
FIG. 2. Fidelity of (a) SORDOR β = π 2 pulses, and (b) SOR-DOR β = π pulses, over a grid of Q coefficients and bandwidth factors.The dashed lines in panels (a) and (b) are the maximum fidelities at each bandwidth factor.Plot (c) shows these maximum SORDOR fidelities compared with the previous best universal rotation pulses (BURBOP).Dashed lines in panel (c) are reference plots exp(−0.50b)and exp(−0.25b).

FIG. 3 .
FIG. 3. A selection of SORDOR pulses showing the unwrapped phase of the constant amplitude pulse as a function of time, measured from the center of the pulse.(a) β = π 2 pulses at six pulse durations; (b) the same pulses shown with a quadratic phase subtracted.(c) β = π pulses at six different pulse durations; (d) the same pulses shown with a quadratic phase subtracted.The formula for the quadratic phase, shown as a dashed line in panels (b) and (c), is outlined in the main text.
) and 4(b) depict the performance of SORDOR π 2 | x and π | x pulses (subscripts denote the axis about which an ω = 0 rotation is defined) projected onto a Bloch sphere for three different initial states.