Optimal pulse design in quantum control: A unified computational method

Abstract

Compelling applications for quantum control have received particular attention and have motivated seminal studies in wide-ranging areas from coherent spectroscopy and MRI to quantum optics. Designing and implementing time-varying excitations (rf pulses) to manipulate complex dynamics of a large quantum ensemble on the order of Avogadro’s number is a long-standing problem and an indispensable step that enables every application of quantum control (1). For example, magnetic resonance applications often suffer from imperfections such as inhomogeneity in the static magnetic field (B₀ inhomogeneity) and in the applied rf field (rf inhomogeneity). In addition, there is dispersion in the Larmor frequency of spins due to chemical shifts. A good pulse design strategy must be robust to these effects, and such variations need to be considered in the modeling and pulse design stages in order to match theoretical predictions to experimental outcomes. As difficult experiments with more demanding performance specifications have emerged, the complexity of finding optimal pulse sequences has drastically increased. For example, as high-field NMR spectrometers are increasingly more accessible and required, broadband excitation pulses are needed to cover a wide ¹³C chemical-shift range (e.g., up to 40 kHz). In addition, to design excitation and inversion pulses that are practical for a typical NMR spectrometer, methods must accommodate realistic maximum rf power and pulse duration while accomplishing the desired spin evolution. In the majority of cases, the length of an rf pulse is constrained by the fixed delays that dictate a certain coherence transfer. Such limitations and imperfections cause a substantial increase in the complexity of the pulse design problem.

From early work using physical intuition (2, 3) to modern methods like composite pulses (4, 5), an enormous body of pulse sequence design techniques has been proposed over 50 y (6, 7), and the process of innovation is ongoing. Highly customized methods, however, have limited scope, such as the Shinnar–Le Roux algorithm, which is robust to Larmor dispersion but not able to compensate for rf inhomogeneity (8). For relatively simple cases, theoretical methods, such as average Hamiltonian theory, provide intuitive guidelines for constructing pulse sequences (9). Heuristic numerical optimization methods have been used extensively for the design of single and multiple pulses in a pulse sequence (10). However, these approaches have a number of shortcomings, such as slow convergence rates and being easily trapped into local optima. In recent years, there have been attempts to look at pulse design problems from a control theory perspective (11–14). In particular, state-of-the-art methods such as gradient ascent and Krotov algorithms are based on principles of optimal control theory (15, 16) and have been used successfully to design broadband and relaxation-optimized pulses that maximize the performance of quantum systems in the presence of relaxation (17–19). These algorithms, although effective, rely on intensive computations, as for system propagators and gradients, as well as a large number of discretizations in the time domain over which to evolve the system.

To overcome these defects, in this article, we provide a systematic framework for optimal pulse design in quantum control. We first present a general mathematical model for pulse design as an optimal control problem of a continuum of bilinear systems. Employing Lie algebra tools from our prior theoretical work (20), we show the control problem of pulse design can be mapped to a problem of polynomial approximation. This notion guides us to develop a unified computational method for optimal pulse design based on multidimensional pseudospectral (PS) approximations, by which a continuous-time optimal control problem of pulse design can be discretized to a constrained optimization problem using interpolating polynomials. The presented multidimensional PS method is a natural extension of the PS method (21) to consider an ensemble of systems. We present simulation and experimental results of optimal broadband excitation and inversion pulses as applied to protein NMR spectroscopy derived by the proposed PS methods.

Theory

Control Problem of Robust Pulse Design.

The statistical properties of an ensemble of quantum systems are described by its density matrix. Under the Markovian approximation, the evolution of the density matrix of an open quantum system (i.e., in the presence of relaxation that is the dissipation due to random interactions between the system and its environment) is governed by the master equation in the Lindblad form,

[1]

where H(t) is the system Hamiltonian that generates unitary evolution, [ , ] is the matrix commutator, and L(ρ) models relaxation (nonunitary dynamics) (22). The Hamiltonian H(t) can be decomposed into H(t) = H_f + H_c(t), where H_f is the free evolution Hamiltonian and is the control Hamiltonian used to manipulate the system by application of electromagnetic pulses u_i(t) of appropriate shape and duration. A typical problem in quantum control is to steer a system from some initial density operator ρ(0) at t = 0 to, or as close as possible to, a desired target operator C. This is achieved by maximizing the expectation value of C, 〈C〉 = tr{Cρ(t)}, for any time t∈[0,T], in the presence of relaxation. Combining the master equation as in Eq. 1 and the relation yields a system of ordinary differential equations that describe the time evolution of an open quantum system for the desired transfer ρ(0) → C (23),

[2]

where is the state vector for t∈[0,T], whose elements x_i, i = 1,…,n, are expectation values of the operators participating in the transfer; is a matrix representation related to the operator H_f and the superoperator L; and are related to the operators H_i.

The above Eq. 2 is an idealization assuming that is the same for all members of a quantum ensemble and that each system experiences the same effect of the applied control field u_i(t). In practice, however, (e.g., NMR experiments) there may be variations in the matrices and across the spin population due to differences in chemical environments, relaxation rates, and applied control fields. Practical considerations including power and time limitations, path and terminal constraints as well as those variations in the system Hamiltonian give rise to a new class of optimal control problems of parameterized bilinear systems,

[3]

where the cost functional contains a terminal cost φ, depending on the final state of the system at the terminal time t = T, and a running cost , depending on the time history of the state and control variables; , , is a parameter vector characterizing variations in the system dynamics, and Ω is a compact d-dimensional interval; , are square matrices; and e and g denote endpoint conditions and path constraints, respectively. Note that in more general cases, the system and control Hamiltonian can be time-dependent, and hence and , due to, e.g., random fluctuations.

Pulse design problems modeled as in Eq. 3 are in general analytically intractable. It is then imperative to develop robust numerical methods to solve these problems. Before computing an optimal pulse, understanding whether a pulse exists for a desired transfer is of fundamental importance.

Motivation for RF Pulse Design via Polynomial Approximations.

A representative example in the control of a continuum of dynamical systems as modeled in Eq. 3 is broadband pulse design in NMR spectroscopy and imaging. The evolution of the bulk magnetization of a sample of nuclear spins in a static magnetic field, B₀, follows the Bloch equations (7). In NMR experiments, one needs to simultaneously excite spin populations, on the order of Avogadro’s number, between desired states of interest with dispersion in their Larmor frequencies and in the presence of rf and B₀ inhomogeneity. Mathematically, it is then natural to consider a continuum of Bloch systems parameterized by and ε∈[1 - δ,1 + δ], 0 < δ < 1, modeling Larmor dispersion and rf inhomogeneity, respectively; i.e.,

[4]

where M(t,ω,ε) denotes the magnetization vector at time t of the spin characterized by the parameter vector s = (ω,ε), (B_1x, B_1y) is the respective rf pulse applied in the x and y axis, and Ω_x, Ω_y, and Ω_z are the generators of rotation around the x, y, and z axis, respectively. In the context of NMR and MRI, for example, a control law (B_1x, B_1y) that makes the transfer from a constant function M(0,ω,ε) = (0,0,1) to another constant function M(T,ω,ε) = (1,0,0) corresponds to a broadband pulse with the duration T in the presence of rf inhomogeneity. The existence of rf pulses that achieve a desired transfer is then related to the controllability property of the system as in Eq. 4. This system is said to be controllable if for any given pair of initial and final states there exists at least one control law (B_1x, B_1y) that drives the system between these two states.

The controllability of a family of Bloch systems of Avogadro-scale, mathematically represented as in Eq. 4, can be understood by studying the Lie algebra generated by the set of matrices {ωΩ_z,εΩ_y,εΩ_x}, where ω∈[ω₁,ω₂] and ε∈[1 - δ,1 + δ]. Performing recursive Lie brackets of these matrices, e.g., [ωΩ_z,εΩ_y], [ωΩ_z,[ωΩ_z,[ωΩ_z,εΩ_y]]], etc., yields new generators of the kind {ω^mε^2ℓ+1Ω_x}, m, ℓ = 0,1,2,…. Using these generators, we can produce rotation of the form,

[5]

simply by piecewise constant controls (SI Text). Therefore, producing any desired (ω,ε)-dependent rotation around the x axis, exp{α(ω,ε)Ω_x}, is possible by appropriate choice of the coefficients t_mℓ in Eq. 5 such that and over [ω₁,ω₂] × [1 - δ,1 + δ]. Similar Lie bracket construction of higher orders and the subsequent polynomial approximations can be applied to produce (ω,ε)-dependent rotation around the y axis. The ability of synthesizing these parameter-dependent rotations allows us to generate any three-dimensional rotation, Θ(ω,ε), by Euler angle decomposition Θ(ω,ε) = exp{α(ω,ε)Ω_x} exp{β(ω,ε)Ω_y} exp{γ(ω,ε)Ω_x}. Therefore, there exists a pulse that is able to produce any arbitrary rotation dependent on the Larmor dispersion, ω, and rf inhomogeneity, ε.

Although the procedure of implementing successive Lie bracketing to generate polynomial evolutions as in Eq. 5 introduces a constructive scheme to design an actual rf pulse, it is not efficient. These ideas, however, lead us to develop a numerical method based on PS approximations.

A Unified Computational Method for Pulse Design.

The idea of converting the problem of pulse design into a problem of polynomial approximation suggests that constructing relevant polynomials should play a central role in a unified computational method for control of quantum ensembles. Solving optimal ensemble control problems as in Eq. 3 numerically requires not only discretization in the time domain t∈[0,T], but also sampling in the parameter domain . To fulfill these requirements, we use PS techniques to approximate the continuous state and control functions of the system in Eq. 3 with interpolating polynomials. The PS method employs spectral collocation in which the differential equation describing the state dynamics is enforced at specific nodes. Developed to solve partial differential equations, these methods have been recently adopted to solve optimal control problems (21, 24).

We focus on Legendre PS methods and, without loss of generality, consider the transformed optimal control problem on the time domain t∈[-1,1]. The underlying dependence on the Legendre orthogonal polynomial basis leads to the property of spectral accuracy, which is analogous to the Fourier series expansion for periodic functions (25). The Chebyshev equioscillation theorem states that for a fixed order an interpolating polynomial is the best approximating polynomial to a continuous function on the interval of [-1,1], as evaluated by the uniform norm. The central idea is, therefore, to approximate the continuous state and control functions, x(t) and u(t), by Nth order interpolating polynomials, I_Nx(t) and I_Nu(t), in terms of the Lagrange polynomial basis ℓ_k(t),

[6]

[7]

By using the Legendre–Gauss–Lobatto (LGL) nodes as the collocation nodes, Eqs. 6 and 7 achieve close to optimal interpolation error (25). The LGL nodes, t_j, j = 0,…,N, are defined by the union of the endpoints, {t₀ = -1,t_N = 1}, and the roots of the derivative of the Nth order Legendre polynomial. By the collocation property, ℓ_k(t_j) = δ_kj, of the Lagrange polynomials we have and for j = 0,…,N. Note that the LGL node distribution is nonuniform, and the high density of nodes near the endpoints is one of the key properties of PS discretizations by which the Runge phenomenon is effectively suppressed (21). The derivative of I_Nx(t) at the LGL node t_j is given by

[8]

where D_jk are elements of the constant (N + 1) × (N + 1) differentiation matrix (25). With Eq. 8 and the natural discretization of the dynamics at the LGL nodes, the differential equations describing the system as in Eq. 3 can be transformed to a series of algebraic equality constraints.

In addition, the integral cost functional in the optimal control problem (Eq. 3) can be approximated by the LGL integration rule, as will be presented in Eq. 11. This approximation is exact for , polynomials of degree at most 2N - 1, when is evaluated at the LGL nodes (21).

PS Discretization and Optimal Sampling.

A multidimensional extension can be devised to approximate the state and control functions with variations, x(t, s) and u(t, s), . For example, for , x(t, s) can be approximated by a two-dimensional interpolating polynomial,

[9]

with the derivative at the LGL nodes t_i and s_j in the t and s domain, respectively,

[10]

where . It is straightforward to extend the same concept to the case of higher dimension for , d > 1. Applying the Gauss–Lobatto integration with Eqs. 9 and 10, we finally arrive at the following finite-dimensional constrained minimization problem (for ):

[11]

where and are matrices representing the Hamiltonians with respect to the rth parameterized system; and are the LGL integration weights with respect to the time and parameter domains, respectively; and is the value of the control function u_i at the jth LGL node t_j. Solvers for this type of constrained minimization problem (e.g., KNITRO, SNOPT) and implementation environments (e.g., MATLAB, AMPL) are readily available and straightforward to implement. See SI Text for more details about the PS method and its implementation.

Simulation Results

A fundamental challenge in pulse design for magnetic resonance applications is to develop robust pulses that compensate for the Larmor dispersion and rf inhomogeneity over the sample. We demonstrate our method by the development of broadband π/2 and π pulses in the presence of rf inhomogeneity. The derived pulses have been implemented experimentally, which verifies both the numerical results as well as the PS method as an effective technique for pulse design.

We convert the related optimal control problem, which describes such broadband pulse design objectives, through PS discretization and sampling to yield a constrained nonlinear optimization problem. As presented in Eq. 4, in the presence of Larmor dispersion and rf inhomogeneity, the bulk magnetization of spins can be represented as M(t,ω,ε). Using the PS method, this magnetization can be approximated by an interpolating polynomial of three variables; i.e., M(t,ω,ε) ≈ I_N×Nω×NεM(t,ω,ε), with the derivative at the LGL nodes t_i, ω_j, and ε_k in the t, ω, and ε domain, respectively,

where .

In a π/2 broadband pulse design optimization, we take the objective, or performance, unless otherwise noted, to maximize the average of the x component of the terminal states. Fig. 1 plots the optimal broadband excitation pulse shape for a bandwidth of [-20,20] kHz with limited amplitude (less than 20 kHz) and maximum duration of 100 μs. With relatively coarse discretization (N = 30) and sampling (N_ω = 8) chosen with the LGL nodes, the PS method achieves a performance of 0.98. In the commonly used gradient methods, optimizations are usually performed in 0.5-μs steps over time (corresponding to N = 200) and 0.5 kHz over the resonance offset space (N_ω = 80), which may achieve a comparable performance (13). Increasing the number of discretizations or samples in the PS method, as shown in the convergence section, would yield an increased objective value. In addition to achieving high performance values, the PS method requires no offline computation of gradients or system propagators, which are essential for gradient methods. Such calculations can be time consuming and create a barrier of entry for experimentalists seeking to apply these methods to solve their problems. For a detailed comparison of the PS method with other numerical algorithms see SI Text.

Open in a separate window

Fig. 1.

The calculated shape of the optimal broadband π/2 pulse and the simulated excitation profile as a function of resonance offset. (A) The optimal broadband pulse shape was derived by the multivariate PS method with the respective number of discretizations in time and frequency, N = 30 and N_ω = 8, to cover the bandwidth [-20,20] kHz with limited rf amplitude for all t and maximum duration T = 100 μs. (B) The corresponding excitation profile has an average x value of 0.98.

The PS method, as a flexible framework for solving optimal control problems, can easily be modified to solve other types of objective criteria, including minimum-energy and time-optimal pulses. Fig. 2, for example, shows the pulse shape and the corresponding simulated excitation profile as designed to minimize rf energy, subject to the performance requirement of 0.98, while compensating for the bandwidth ω∈[-20,20] kHz as well as 10% rf inhomogeneity with maximum duration of 100 μs. The pulse achieves an average performance of 0.98 with N = 24, N_ω = 8, and N_ε = 1. An optimal pulse designed with the same parameters, but with an objective to maximize the average of the x component of the terminal states illustrates that the pulse in Fig. 2 achieves the same performance (within 0.3%) with 16% less rf energy (see SI Text). Note the PS method naturally finds continuous and smooth pulses because of the use of polynomial approximations, which are experimentally preferred. The adaptability and generality of the PS method is what makes it an excellent candidate for a universal and unified method for pulse design.

Open in a separate window

Fig. 2.

The calculated shape of the optimal broadband π/2 pulse and the simulated excitation profile as a function of resonance offset and rf-inhomogeneity. (A) The minimum-energy broadband pulse shape was derived by the multivariate PS method with the respective number of discretizations in time, frequency, and rf-inhomogeneity, N = 24, N_ω = 8, and N_ϵ = 1, to cover the bandwidth [-20,20] kHz with 10% rf-inhomogeneity and with limited rf amplitude B₁(t) ≤ 20 kHz for all t and maximum duration T = 100 μs. (B) The corresponding excitation profile is shown with an average x value of 0.98.

It is of fundamental importance to show that the solution achieved from the discretized problem as in Eq. 11 is still a solution to the original continuous problem as in Eq. 3. In particular, as the number of discretizations (or samples) gets large, the objective value of the discretized optimization problem should converge to the objective value of an optimal solution of the original continuous problem. Figs. 3 and ​and44 present the numerical convergence of the performance of a broadband π/2 pulse over increasing N and N_ω, respectively. These figures also illustrate the fast convergence rate and low discretization number characteristic of the PS method.

Open in a separate window

Fig. 3.

Numerical convergence of PS discretization in the time domain. The performance of the optimal broadband π/2 pulse converges to unity as the number of discretizations, N, gets large. Parameter values: maximum rf amplitude = 20 kHz, bandwidth = [-20,20] kHz, duration = 100 μs, and fixed N_ω = 8.

Open in a separate window

Fig. 4.

Numerical convergence of optimal sampling. The performance of the optimal broadband π/2 pulse converges to unity as the number of samples, N_ω, gets large. Parameter values: maximum rf amplitude = 20 kHz, bandwidth = [-20,20] kHz, duration = 100 μs, and fixed N = 30.

Experimental Results

The excitation profile of the optimal broadband π/2 pulse presented in Fig. 1A is shown in Fig. 5, which was recorded with resonance offset ranging from -30 to 30 kHz in steps of 2 kHz. The result demonstrates a uniform excitation over the bandwidth [-20,20] kHz as designed. The optimal pulse derived by the multidimensional PS method has a noticeably higher average excitation profile than the conventional square pulse, especially near the edges of the bandwidth. The phase errors at the edge of the designed excitation range are less than 10° for the optimal pulse and greater than 60° for the hard pulse. Fig. 6A illustrates an optimal broadband inversion pulse designed to cover the bandwidth [-40,40] kHz with its amplitude limited under 20 kHz. The corresponding inversion profile recorded with resonance offset ranging from -60 to 60 kHz in steps of 2 kHz is shown in Fig. 6B and highlights the larger difference between the resonance offset profiles for the hard pulse (black) and optimal pulse (red) for the case of the longer duration and wider bandwidth inversion pulse.

Open in a separate window

Fig. 5.

Experimental excitation profiles of broadband π/2 pulses. The excitation profiles correspond to the optimal broadband pulse (red) in Fig. 1 (designed to cover the bandwidth [-20,20] kHz) and a conventional hard pulse (black).

Open in a separate window

Fig. 6.

Experimental excitation profiles of broadband π pulses. (A) The optimal broadband π pulse shape was derived by the multivariate PS method with the respective number of discretizations in time and frequency, N = 36 and N_ω = 12, to cover the bandwidth [-40,40] kHz with limited rf amplitude B₁(t) ≤ 20 kHz for all t and maximum duration T = 120 μs. (B) The excitation profiles correspond to the optimal broadband pulse (red) and conventional hard pulse (black).

Furthermore, we show that these pulses can be readily incorporated in standard pulse sequences routinely used in solution state NMR spectroscopy. In high-field magnets (> 700 MHz), it is challenging to uniformly excite the entire spectral range with a standard hard pulse. This is particularly difficult when exciting the carbon nuclei in protein NMR spectroscopy where the chemical shifts of interest range from 10 ppm (methyls) to 140 ppm (aromatics). This 130-ppm range on a 900-MHz spectrometer is approximately 30 kHz. The inability of a standard hard pulse to cover this bandwidth is readily observable in a ¹H-¹³C correlated heteronuclear single quantum coherence (HSQC) spectrum. In Fig. 7A we show that the peaks in the downfield region (60–80 ppm) of the ¹H-¹³C correlated HSQC spectra of [¹³C,¹⁵N] GB1 were still observable when the carrier frequency of the optimal inversion pulse was shifted to the up-field (31 kHz). However, when we used an adiabatic ¹³C inversion pulse (Fig. 7B) or a conventional hard π pulse (Fig. 7C) at a resonance offset of 250 ppm (31 kHz), a majority of the resonances either were absent or significantly diminished in intensity. The optimal pulse yields approximately a 20 times sensitivity enhancement over the adiabatic pulse (see SI Text) with a much shorter duration (120 versus 500 μs).

Open in a separate window

Fig. 7.

HSQC spectra comparison. The ¹H-¹³C correlated HSQC spectra of [¹³C,¹⁵N] GB1 sample using optimal inversion pulses (A), standard Bruker adiabatic inversion pulses (Crp60.0.5.20.1) (B), and conventional hard inversions (C). Spectra on the left were recorded when the ¹³C-inversion pulses were on resonance (at 50 ppm), and the spectra on the right were recorded when the ¹³C-inversion pulses were 31 kHz (at -200 ppm) off resonance to the up-field.

Conclusion

Manipulating a large ensemble of dynamical systems is one of the key steps in many compelling applications such as protein NMR spectroscopy, quantum computing, and quantum optics. In this article, we introduced the notion of mapping pulse design to a problem of polynomial approximation by the use of tools from control theory. Accordingly, we developed a unified computational method for optimal pulse design based on PS approximations, which features coarse discretization with spectral accuracy producing smooth solutions. This computational method provides flexibilities in designing optimal pulses for any desired objective and can be applied to consider any initial distribution of system states. Optimal broadband pulses with reduced rf energy and duration derived by this method have been implemented in protein NMR spectroscopy, and the experimental excitation profiles demonstrate significant performance improvements over conventional hard and adiabatic pulses. The method is not restricted to the systems studied here but is universal to closed and open quantum systems and further to broad areas of control of ensembles, such as those arising in systems biology and neuroscience (26).

Materials and Methods

Supplementary Material

Acknowledgments.

Footnotes

References