Quantum robust control theory for Hamiltonian and control field uncertainty

Quantum robust control—which can employ fast leading order approximations, slower but more accurate asymptotic methods, or a combination thereof for quantification of robustness—enables control of moments of quantum observables and gates in the presence of Hamiltonian uncertainty or field noise. In this paper, we present a generalized quantum robust control theory that extends the previously described theory of quantum robust control in several important ways. We present robust control theory for control of any moment of arbitrary quantum control objectives, introducing moment-generating functions and transfer functions for quantum robust control that generalize the tools of frequency domain response theory to quantum systems, and extend the Pontryagin maximum principle for quantum control to control optimization in the presence of noise in the manipulated amplitudes or phases used to shape the control field. To provide guidelines as to the types of quantum control systems and control objectives for which asymptotic robustness analysis is important for accuracy, we introduce methods for assessing the Lie algebraic depth of quantum control systems, and illustrate through examples drawn from quantum information processing how such accurate methods for quantification of robustness to noise and uncertainty are more important for control strategies that exploit higher order quantum pathways. In addition, we define the relationship between leading order Taylor expansions and asymptotic estimates for quantum control moments in the presence of Hamiltonian uncertainty and field noise, and apply such leading order approximations to significant pathways analysis and dimensionality reduction of asymptotic quantum robust control calculations, describing numerical methods for implementation of these calculations.


Introduction
Quantum robust control techniques are concerned with the control of moments of quantum states, gates and observables with respect to sources of uncertainty relevant to the manipulation of atomic and molecular dynamics [1][2][3][4][5], in addition to moments due to quantum noise and measurement uncertainty that are optimized by quantum optimal control [10]. Such approaches are necessary foundations for model-based quantum control of molecular dynamics [6][7][8][9], or adaptive feedback approaches that combine model-based strategies with learning control [11,12]. Inspired by corresponding concepts in the domain of classical control engineering, quantum control robustness is described in terms of the effect of disturbances in the control field [13,14] and limited knowledge of system Hamiltonian parameters [15,16] on the different moments of the quantum control objective [17][18][19]. The concept of robustness is essential in ensuring control objectives are satisfied in the presence of such noise and uncertainty [20][21][22].
While evolutionary algorithms have been successfully implemented in laboratory quantum control to obtain optimal control profiles [23][24][25], this approach of experimental quantum learning control has primarily been applied to the problem of population transfer between eigenstates [10], whereas applications to control of mixed states are more challenging to achieve without employing model-based robust control. In addition to its applications in the control of quantum observables, robust control of quantum systems has important applications in quantum information processing. One of the important problems in this area is optimal protection of quantum systems against decoherence. However, control field noise and Hamiltonian parameter uncertainty can result in equally significant deviations from the target quantum gate. Applications of quantum optimal control theory (QOCT) to quantum computation include optimal implementation of quantum gates in closed systems [26][27][28][29] and in open systems sustaining environmental decoherence [30][31][32]. In particular, model-based QOCT techniques have been applied in designing radiofrequency pulses for implementation of quantum gates in systems of coupled nuclear spins in nuclear magnetic resonance experiments [33]. For example, shaped pulses obtained from QOCT have been applied to generate single-qubit gates in trapped ion systems with improved robustness to control field noise [34].
Several robust quantum control theories and algorithms, inspired by classical counterparts, have been studied [3,4,[35][36][37][38][39], including robust control via leading order Taylor series [4,37], either for the control of moments or in minimax/worst-case formulations [18]. However, these methods are limited in their accuracy unless the magnitude of uncertainty is small. In addition, whereas in classical control engineering these methods are typically applied with real-time feedback, the latter is currently not possible for many quantum systems. Another important tool in the domain of engineering control for assessing the response of a system to input noise is the method of transfer functions [40]. Transfer functions are used to approximate the moments of the system state vector components or an observable function of the state, when the system is subject to input noise of a specified frequency domain spectrum. While this approach of frequency domain response is commonly applied to linearized control systems, it has not been applied to bilinear systems such as quantum control systems.
In prior work, we introduced the asymptotic theory of quantum robust control, which enables accurate control of moments of quantum observables and gates in the presence of Hamiltonian uncertainty or field noise. In our initial work on quantum robust control [17], we provided an asymptotic theory for computation and control of moments of bilinear quantum systems in the presence of common types of Hamiltonian and control field uncertainty, without relying on linearization or related leading order Taylor expansions. The robustness of the quantum dynamics was analyzed with respect to implementation errors in the classical input variables (in a semiclassical picture of controlled quantum dynamics) and parameter uncertainty in the quantum Hamiltonian. Specifically, the method calculates the effect of noise in the control field and uncertainty in the system's dipole moment on the fidelity of control. In addition, different quantum pathways involved in the controlled dynamics are delineated such that moments of quantum interferences between different order pathways and their contributions to the transition amplitude and probability can be quantified for increased accuracy.
More recent work from the authors [12], which applied these robustness analysis techniques to control optimization in the presence of noise and uncertainty, introduced the Pontryagin maximum principle (PMP) for the optimality conditions for quantum control objective moments in the presence of noise and uncertainty, and also presented evolutionary algorithm-based quantum robust control design based on the quantum robust control PMP. The proposed control strategies were classified into the control paradigms of model-based and model-free robust control, which have the potential for integration through adaptive feedback techniques.
In this paper, we present a generalized quantum robust control theory that extends the previously described theory in several important ways. We present robust control theory for control of any moment of arbitrary quantum control objectives, introducing moment-generating functions and transfer functions for quantum robust control that generalize the engineering tools of frequency domain response to quantum systems, and extend the PMP for quantum robust control to control optimization in the presence of (discrete) field noise. We establish criteria in terms of a quantum control system's dynamical Lie algebra that are relevant to the robustness of control strategies. In addition, we examine the relationship between leading order Taylor expansions for several quantum control moments and asymptotic estimates of these moments, and apply such leading order approximations to significant pathways analysis and dimensionality reduction of asymptotic quantum robust control calculations, introducing bounds on number of significant pathways that need to be considered in such calculations as well as numerical methods for these calculations. We also extend quantum control robustness analysis to include phase encoding, which is required for asymptotic robustness calculations in the presence of noise in the manipulated phases of control fields.
The paper is organized as follows. In section 2, we review encoding schemes for quantum control robustness analysis and introduce phase encoding. In section 3, we extend the asymptotic quantum robustness analysis theory to moments of arbitrary quantum control objectives. In section 4, we introduce moment-generating functions and transfer functions for quantum robust control that can accommodate arbitrary distributions of control field noise and Hamiltonian parameter uncertainty. In section 5, we provide examples from the control of quantum gates in molecular systems and establish criteria based on a control's system's dynamical Lie algebra that play an important role in determining the accuracy of robust control techniques. In section 6, we introduce the PMP for quantum robust control in presence of noise in the manipulated amplitudes and phases of the control field. In section 7, we introduce leading order approximations for quantum robust control and summarize properties of the Hessian of the quantum control landscape that are relevant to robust control. We also present methods for significant pathway identification and dimensionality reduction that leverage properties of the quantum control landscape. Finally, in the appendix we conclude by providing bounds on the number of significant pathways relevant to quantum control moment calculations, which is relevant for numerical implementations of the methodology.

Robustness of controlled quantum dynamics
Asymptotic quantum control robustness analysis constitutes an approach that can provide accurate estimates of the first and second moments of quantum control objectives (and higher moments if desired) of δJ (the deviation in the control objective functional J due to Hamiltonian uncertainty or field noise) suitable for use in either distributional or worst-case robustness criteria for controlled quantum dynamics, along with bounds on the accuracy of those estimates, this approach is more accurate than methods for moment calculations based on leading order Taylor expansions. In prior work [17], we introduced expressions for moments of quantum transition amplitudes and probabilities in terms of quantum pathways corresponding to two types of uncertainty that can affect the predicted dynamics: control field amplitude pathways, and Hamiltonian (dipole parameter) pathways. The equations for field amplitude pathways are summarized in appendix A. Here, we review dipole pathways and introduce control field phase pathways, in order to complete the specification of all major types of pathways relevant to robust quantum control.
For each type of pathway, there is a corresponding encoding scheme that enables the contributions of the pathways to the quantum control objective moments to be determined. Thus through time-independent Hamiltonian and control field encoding of quantum pathways, and knowledge of the moments of Hamiltonian and control field parameters, it is possible to compute the moments of quantum control objectives to arbitrary degrees of accuracy.
Recall that dipole pathways are defined as follows: where the shorthand notations A(ω k ) = A k and ω ji = (E j −E i ) have been used, and where the sum (1,...,lm−1) is over all 1 l i N, i = 1, . . . , m − 1 such that frequency ±ω pq corresponding to dipole parameters μ pq , μ qp appears in the multiple integral α pq times. Phase pathways are defined as those quantum pathways that differ only in the way in which control field phases enter the pathway expressions. Like amplitude pathways, they are used to calculate the effect of noise in a manipulated input variable on the moments of a quantum control objective. For phase pathways, . . , k m ) denotes the number of times mode k appears in the multiple integral, and m min = K k=1 |α k |. Note that a particular combination of phases in a phase pathway does not uniquely specify the pathway order, unlike amplitude and dipole pathways. Normalizing quantum control pathways allows them to be used in moment calculations, as described in [17]. Here, we have: In order to obtain the contributions of phase pathways to a quantum control objective function value, a modulation scheme is required. For phase encoding, the modulation scheme is as follows: where s denotes a timelike modulation parameter, and the transition amplitude consists of: Deconvolution of the transition amplitude term yields mth order phase pathways in a way identical to the amplitude counterpart. Note that in the case of phase modulation, the encoding is Hermitian.

Generating functions and transfer functions for quantum control moment calculations
In our previous work on asymptotic quantum control robustness analysis [17], expressions for the first and second moments of the transition amplitude and transition probability for state-to-state population transfer were derived. However, it is possible to express the moments of any quantum control objective function (including quantum observable or gate control fidelity) more generally in terms of a type of moment-generating function, through two Fourier transforms. The characteristic or moment-generating function provides moments of a distribution in terms of derivatives of its Fourier transform. Defining and assuming amplitude noise that is uncorrelated between modes, the generating function for the expectation of the nth power of F(U(T)) can be written with the expectation expressed in terms of the generating function as where Γ = 0, . . . , γ f denotes the set of encoding frequencies corresponding to the significant pathways.
The quantity with A a spectral amplitude is an order γ transfer function for the specified moment of F. When multiplied by a corresponding input noise moment obtained from the generating function of the noise input, provides the resulting contribution of that order pathway to the moment of F. Note that by changing the γ to run over a subset of Γ, we can extract expectations of specified interferences (in the case n = 1) or any other subset of terms. Similarly, the PMP first-order conditions for optimality for moments (see below) can be expressed in this form as well.
In general below, the notation J(ε(t)) is often used for the quantum control objective expressed as a functional of the control field (whereas F is a function of the final time dynamical propagator).

General noise and uncertainty distributions
The expressions for E [J], var J require expressions for higher moments of noisy/uncertain manipulated input and system parameters like A k , μ ij . Higher moments can be provided in closed form for any uncertainty distribution for which there exists an analytical Fourier transform. To obtain these, we apply the characteristic (moment-generating) function M(y) of the probability distribution function p(x) of the manipulated input or system parameter: If M(y) is available in closed form, the moments of x are obtained via (using the notation · to denote expectation for simplicity of notation) .
Consider the case of Gaussian uncertainty/noise in each of the parameters which is i.i.d: i.e., Then In this case, the moments are given by: For example, Applying the binomial expansion, we obtain any moment (x −x) m . For any odd moment, we find and hence (x −x) m = 0. If the parameter distributions are independent but not identically distributed, we use the following decomposition of the total characteristic function M y : where y = y 1 , . . . , y K T . When the uncertainty/noise in the parameters are not independent (i.i.d.), a multivariate Gaussian is used, in which case we have: For example for this multivariate case withx = 0, we have according to Wick's/Isserliv's theorem [41]:

Examples: quantum gate control in molecular systems
In [12], we presented numerical examples of the implementation of asymptotic robustness analysis and multiobjective Pareto optimization of the first and second moments of population transfer in atomic systems. Here, we illustrate how to extract quantum pathways for optimal quantum gate control in molecular systems of various types and Hilbert space dimensions, for use in robustness analysis and robust optimization for such applications. We consider the control of two types of molecular states that have been proposed for use in quantum information processing, namely rotational and vibrational states.
1. Rotational: one of the initial proposals towards physical quantum information processing, was to use the spin of electrons as qubits-where spin-up, |1 , can be mapped to 1 and spin-down, |0 , to 0. Spin of an electron is intrinsic angular momentum, and the central point is to use the spin angular momentum vector as qubit. This prompts the use of rotational angular momentum vector of the entire molecule towards quantum computation. This may be particularly interesting because, rotational states have been extensively studied using spectroscopic methods and are also easy to be manipulated in lab with the existing technology within their relatively long dephasing times.
Use of rotational modes of a molecular system as qubits and performing quantum information processing has been proposed [42][43][44], wherein the authors make use of multi-target optimal control theory. A typical rotational system consists of a diatomic molecule with a reduced mass m and separated by a distance r. In the absence of external interactions such a molecule with N energy levels, will have the following form for the internal Hamiltonian in the rotational energy eigenbasis: where j is the eigenvalue of the total angular momentum,Ĵ. The z-componentĴ z with eigenvalue m takes 2j + 1 values. The qubits can be represented by any eigenstate |j, m . For example for j = 1, we have m = −1, 0, 1 and we can map |1, 0 as |0 and |1, 1 as |1 .
In order to perform a quantum gate, we need to interact the system with a control field designed for the specific gate. Such a molecular system interacts with external electromagnetic field, ε(t) due to inherent dipole moment, which formally is represented by an N × N symmetric matrix μ. The symmetricity of this matrix is imposed due to the selection rules, where the transitions are only allowed to the neighboring states implying Δj = ±1 and leading to the matrix elements j|μ|j ∝ δ j,j±1 . Thus the interaction Hamiltonian is where the minus sign implies that, due to the torque created by the field, the dipole tends to align with the interacting electric field.
2. Vibrational: Tesch et al [45][46][47] have proposed the use of vibrational modes of the molecule towards quantum information processing. In this case, the vibronic excitations, which correspond to normal modes, are used to define the qubit. Each such mode will have several excitations, and pairs of such can be used to represent |0 and |1 . The implementation of logic gates is performed by interacting the molecule with pre-designed laser pulses, which are calculated using the OCT. Although the formal method of finding the optimal pulse shapes will be the same, the physical implementation will be done using different frequencies of the electromagnetic spectrum. Whereas in the case of rotational states, microwave frequencies may be used, for vibrational states, usually the IR frequency is used. The idea of using harmonic oscillator eigenstates for quantum computation was described in [52]. However the harmonic oscillator potential is superseded by what is called the Morse potential, which is a more realistic description of molecular interactions [48]. The Morse potential V(r − r e ) = D e 1 − e −β(r−r e ) 2 , an excellent approximation of the ground electronic state behavior of a vibrating diatomic molecule, is employed. D e is the dissociation energy, ω 0 is the characteristic angular frequency of the oscillator, and β = ω o m 2D e . Its eigenfunctions are where (K + 1/2) 2 = 2mD e / 2 β 2 , z = (2K + 1)e −β(r−r e ) , and the L (α) ν are the generalized Laguerre polynomials. H 0 is again in its diagonal basis with elements corresponding to the system energy levels: Unlike the rigid rotor model, the Morse potential exhibits no strict selection rule; transitions Δν = ±2, ±3, . . . are only weakly forbidden. The simplest method to compute the dipole matrix elements is direct integration over the Morse eigenfunctions: ν|μ(r)|ν = Ψ ν μ(r)Ψ ν dr. We use parameters for the diatomic molecule CO in both the rigid rotor and Morse oscillator model systems above, for the implementation of quantum gates. Similar results were obtained for HCl in the gas phase. In order to render the rotational and vibrational systems comparable, with the salient differences being the spacing of energy levels in H 0 and the coupling structure of the dipole moment operators, we scaled the μ operators of the two systems such that they have equal and constant off-diagonal elements for |j − i| = 1 (where i, j index the Hamiltonian matrix elements; corresponds to allowed transitions in both systems), set all μ elements for a given |j − i| to the value of the largest element on the Hilbert subspace of states chosen for quantum computation, and scaled the dipole operator of the rotational system such that its nonzero elements matched the values of the corresponding elements of the vibrational dipole operator. Then, the same constant of proportionality was applied to scale the rotational H 0 operator as well. (This resulted in the dynamical timescale of the rotational model decreasing to the low ps scale corresponding to infrared frequencies.) This enables analysis of scale-invariant features of the controlled dynamics for the two structures of molecular Hamiltonians, as will be discussed below, with the primary differences between the systems being the ratios of H 0 eigenvalues and the decay rates of the dipole operator elements with |j − i|, which were all maintained at their physical values. Similar model systems were used in prior work on molecular quantum gate control [26].
Quantum gates are then implemented on input states on a subspace of the Hilbert space of bound molecular states. In our examples, we truncate the Hilbert spaces for the rotational and vibrational states to N = 4 or N = 8. Prior literature (e.g., [45][46][47]) has done the same, since for specified field fluences and evolution times, the effective dynamical dimension of the systems does not span the entire Hilbert space of bound states, to a good approximation. This is because the dipole couplings are not sufficient to transfer population between states beyond a relatively small difference in energy under these conditions (which also results in relatively limited quantum decoherence since the dynamical propagator on the employed subspace is roughly unitary). Moreover, for the vibrational systems, the larger energy level spacing compared to rotational systems results in a smaller effective dynamical dimension and an even more limited decoherence. Due to the assumption of a dilute gas phase ultrafast dynamics, the primary source of decoherence is nonunitary dynamics on the subspace used for information processing. Decoherence is not modeled in the present work due to the truncation of the Hilbert space for illustrative purposes.
The associated optimal control problems for quantum information processing in such systems can be solved using so-called homotopy tracking algorithms [49,50]. 3 These algorithms follow a specified track F v of objective function values, where s denotes algorithmic time, toward the global minimum of each objective, while locally minimizing an auxiliary Lagrange cost (functional of the field) at each step. The following differential equation specifies the evolution of each control field ε in continuous algorithmic time: Here a(v, t) denotes the functional derivative (gradient) of the Hilbert-Schmidt distance with respect to the control ε(t), b(v) = T 0 a 2 (v, t)dt is the norm square of the gradient, and f(v, t) is a 'free' function that arises due to the fact that the control problem is underdetermined in the absence of a Lagrange cost [49]. The choice of control optimization algorithm affects the fluence of the resulting optimal gate control field; following the shortest possible (geodesic) path from U 0 to W identifies higher fluence fields. Setting dF v dv = 0 after convergence to W allows exploration of the level set of control fields that drive the system to W at time T. The choice of free function f(s, t) = −ε(s, t) seeks to reduce fluence while traversing the level set. For gate fidelity problems we use F(U(T)) = W − U(T) 2 , the Frobenius distance between the final-time dynamical propagator and the target gate [26]. Figures 1 and 2 present the optimal control fields and corresponding power spectra obtained from control optimization for maximization of Hadamard gate fidelity in two-qubit rotational and vibrational molecular model systems. Similar systems were studied in the context of deterministic quantum gate control in [26]. For the single-qubit Hadamard gate we have 3 Dynamic optimization (optimal control) problems typically require specialized algorithms since the dynamical constraint is a differential equation that must be satisfied for each feasible control; homotopy tracking algorithms are ideal for multiobjective control problems [51]. Optimal field (left) and power spectrum (right) for driving four-level (two-qubit) system (1) to the Hadamard gate on qubit 1 (W 2Q,1 ), with total evolution time T = 1 ps. Equation (3) with was used to solve for the optimal field. The Frobenius distance between controlled and target gates is less than 10 −5 . Hamiltonians were scaled as described in the text for scale-invariant comparison to the vibrational system. We consider control of the following Hadamard gates on several multi-qubit systems: where the first subscript refers to the number of qubits in the control system and the second denotes the qubit on which the Hadamard gate is implemented. Figure 1 shows the optimal control field and power spectrum for the Hadamard gate W 2Q,1 on the first qubit of the two-qubit rotational system, whereas figure 2 shows the results for the Hadamard gate W 2Q,1 on the first qubit of the two-qubit vibrational system. Figures 3 and 4 present the moduli of the extracted Dyson series orders from the optimally controlled dynamics for the rotational and vibrational systems, for particular matrix elements of the unitary propagator. Figure 3 shows the moduli of extracted Dyson series orders for optimal gate control of the Hadamard gates W 2Q,1 and W 2Q,2 on the first and second qubits, respectively, for the rotational two-qubit system. Figure 4 shows the results for optimal gate control of the Hadamard gates W 2Q,1 and W 3Q,1 on the first qubit for the vibrational two-qubit and three-qubit system, respectively. The complex amplitudes of such quantum pathways can be used to calculate the moments of the unitary propagator or gate fidelity in the presence of field noise with arbitrary distributions, as described above and as illustrated numerically in [12] for uncorrelated noise.

Lie algebraic depth and quantum control robustness
From figures 1, 2 and 4, respectively, it is apparent that in this example, the scaled rotational system employs fewer field modes and fewer Dyson series orders to reach the target Hadamard gate. As described   (2) using an optimal control field. The modulus of the true propagator element is within 5% of the sum of significant orders. Frequency gamma denotes the Hamiltonian encoding frequency. above, in the context of robust control, the involvement of more Dyson series orders (quantum pathways) in a dynamical control strategy results in larger contributions from higher moments of the noisy/uncertain parameters and increases the importance of employing more accurate robustness analysis techniques. In order to elucidate features of quantum Hamiltonians that affect the accuracy of robust control, we now consider the differences in the quantum pathways required to control such model systems from the perspective of controllability theory.
Controllability of a physical system implies that, there always exists a control, such that the system can be driven from an initial state to any given final state [52]. Operationally this can be checked using the Lie algebra rank condition (LARC), which states that a system is fully controllable if the associated Lie algebra has full rank, which for U(N ) system is N 2 and for SU(N ) systems is N 2 − 1. The rank condition for full controllability (i.e., any unitary matrix can be produced at some time T) is which is called the dynamical Lie algebra L. A consequence of the LARC is that, one can generate the underlying basis of the corresponding Lie algebra as discussed in [52].  [52] indicates that one has to repeat the commutators until there are no more linearly independent basis elements (matrices). To quantify the number of commutators, a quantity called depth was introduced in [26], and is further developed here.
The interaction dipole moment μ is associated with the external field and consequently the pathways in the Dyson series expansion of the unitary propagator which involves the various integrals of the electric field are associated with the order of μ. This implies that the order of μ within these nested commutators reveals the sensibility of the molecular system to the external electromagnetic field. Lower number of μ means the molecule responds easily and in general requires less energy/fewer field modes and quantum pathways to control. On the other hand higher order of μ means the molecule is more inert to the external field and in general requires more energy/more field modes and quantum pathways to control.
The saturating depth of the dynamical Lie algebra for a given order of μ is defined as the lowest value of k in equation (4) required to maximize the rank of this algebra with a series of commutators of this order in μ. This depth for each order of μ may differ considerably for different physical systems.
The algorithm presented in [52] suggests how to find the basis that spans the Lie algebra, but does not reveal the order of μ, i.e. does not allow determination of how many μ s are involved in nested commutators that span a Lie algebra of given rank. Here we present an algorithm that can reveal the basis and its rank for each order of μ s. Algorithm: In the search for optimal solutions, one would require different energy and different numbers of field modes to locate the solutions in different Lie algebra directions. It may also be possible that two different directions may not differ significantly and so to distinguish these directions requires more energy and field modes The Lie algebraic depths of two-and three-qubit model rotational and vibrational molecular quantum control systems are presented in figure 5. The results for Lie algebraic depth are plotted up to the μ order at which the Lie algebra rank becomes full. Figure 5 top shows the results for the model two-qubit rotational (left) and vibrational (right) systems, whereas figure 5 bottom shows the results for the model three-qubit rotational (left) and vibrational (right) systems. It can be observed that while for the lowest order in μ the vibrational systems have higher Lie algebra rank, this rank is not sufficient for gate control, and rotational systems otherwise have higher rank and lower depth for most orders in μ. The difference in Lie algebraic depth between the rotational and vibrational model systems is magnified for larger numbers of qubits (e.g. three vs two). The higher rank/lower depth of vibrational systems at low order in μ is due to the rapid decay in dipole coupling magnitudes with distance |j − i| between eigenstates, which have increasingly small contributions for brackets of higher order in μ (μ elements with j − i = ±1 dominate for these brackets).
In [26], a computational complexity scaling analysis of quantum control optimization was reported for several canonical gate control systems, including molecular systems proposed for quantum information processing. In particular, this work studied rigid rotor and anharmonic oscillator Hamiltonians of Hilbert space dimensions 4, 8, 16, and 32 (the former two also being studied herein). As in the current work, in order to render the dynamics comparable, that study set corresponding dipole matrix elements with |j − i| = 1 for the two classes of systems to constant, identical values. In that work, the scaling of optimization effort for different dipole structures (rotor, oscillator) was compared for the same H 0 . It was observed that systems with more sparse dipole operators scaled more poorly with respect to optimization effort, and it was conjectured that this was due to a higher Lie algebraic depth. Herein we have described how to explicitly compute the Lie algebraic depths for the different classes of systems and found that in general, for systems with different H 0 's, the dipole operator sparseness is in general not sufficient to predict depth. Nonetheless, consistent with [26] we have confirmed herein that the optimal fields for systems of higher Lie algebraic depth in our examples do employ more field modes and quantum pathways. The same principle applies to control of other types of quantum control objectives, such as state control, since LARCs for controllability exist for each [52]. Top left: rank of the dynamical Lie algebra as a function of the Lie algebraic depth of nested Lie brackets, for various numbers of μ operators appearing in the Lie brackets for a four-level (two-qubit) diatomic rigid rotor model. Top right: rank of the dynamical Lie algebra as a function of Lie algebraic depth, for various numbers of μ operators appearing in the Lie brackets for a four-level (two-qubit) diatomic Morse oscillator model. Bottom left: rank of the dynamical Lie algebra as a function of Lie algebraic depth, for various numbers of μ operators appearing in the Lie brackets for an eight-level (three-qubit) diatomic rigid rotor model. Bottom right: rank of the dynamical Lie algebra as a function of Lie algebraic depth, for various numbers of μ operators appearing in the Lie brackets for an eight-level (three-qubit) diatomic Morse oscillator model. In each case, for every order in μ, the depth is plotted until the Lie algebra rank saturates. Once an order in μ is found such that the rank is full, higher orders in μ are not considered. System parameters for each model are described in the text.
In conclusion, given the fact that quantum control systems with higher Lie algebraic depth generally require higher order quantum pathways for high fidelity control, robust optimization rather than deterministic optimal control is important to achieve high fidelity in the presence of noise or uncertainty. Indeed, it is well-known that robustness of the dynamical propagator for linear control systems can be expressed in terms of only a finite number of Lie brackets. Similarly, local controllability of bilinear quantum systems-which is assessed in terms of the quantum control local controllability Gramian matrix that depends only on the gradient of the dynamical propagator or state vector with respect to the control field [26,53]-does not exploit the complete set of Lie brackets available to the full bilinear system. Given these results, in the next section we introduce the PMP for robust control optimization in the presence of field noise.

Pontryagin maximum principle for quantum robust control
The PMP for quantum control [10,26], used in the development of optimality conditions and control optimization algorithms, can be extended to the control of moments; when all moments are controlled, this amounts to control of the corresponding Fokker-Planck equation for the probability density of the state in the presence of input noise and parameter uncertainty. For linear systems with additive noise, it is possible to obtain an analytical solution for the time evolution of the first and second moments of the state variables (e.g. for the Ornstein-Uhlenbeck process). This is the basis for the linear quadratic (Gaussian) regulator (LQR/LQG) problems in feedback control, where the first and second moments are generally controlled (only relevant moments for Gaussian noise inputs).
The LQG is derived in terms of the Hamilton-Jacobi-Bellman partial differential equations. The LQG formulation is used in the theory of linear quantum filtering and real-time feedback control [54]. For bilinear systems (multiplicative noise) such an analytical solution does not exist. However, it is possible to derive a PMP for quantum control of moments of any objective function in the presence of field or Hamiltonian uncertainty, starting from the robustness analysis theory developed above. This provides the theoretical foundation for model-based robust control of molecular quantum systems. It also has applications to the robust control of other bilinear systems with multiplicative noise.
The state and costate equations for the robust quantum control PMP are partial differential equations. Both state and costate are functions of t and the timelike variable s. A generalized expression for the moments of quantum observables is required in order to derive the PMP (see above). The costate equation and PMP (first-order conditions for optimality) for quantum robust control can be derived analogously to the deterministic quantum control PMP. The costate equation for quantum robust control is a partial differential equation in t, s: T, s)).
Analogously to E[U(T)], the expectation of the quantum control costate at the final time, Φ(T), can be obtained as follows for amplitude noise:

s)∇ U F(U(T, s))U † (t, s)μ I ε(t, s)U(t, s) .
Using the PMP-Hamiltonian function we can obtain first-order conditions for optimality of moments. Expressing the LagrangianJ [26], which incorporates the dynamical constraints, in terms of H and integrating Φ(t), dU(t) dt by parts, we get The first-order variation of this Lagrangian is H provides the equation of motion for the costate that was introduced above. This expression can be evaluated by first writing the s-evolved LagrangianJ and then considering its first-order variation The corresponding first-order conditions (Euler-Lagrange equations) follow from the requirement that δE[J] = 0 for any specified deterministic variation δε, and hence for any deterministic variation δU(t).
In reference [12] we presented the quantum robust control PMP optimality condition in the presence of Hamiltonian parameter uncertainty. For amplitude uncertainty, Here, A k = A (ω k ) denotes a subset of noisy amplitude modes. Note that use of the PMP for amplitude noise implies that the control optimization is executed in the frequency domain through modification any of the amplitudes A(ω), of which the above subset is subject to noise. Appendix B includes numerical methods for the evaluation of frequency domain PMP gradients and provides the expression for the deterministic PMP gradient in the frequency domain, for reference.

Leading order Taylor expansions
As discussed above, quantum control systems of higher saturating Lie algebraic depth are those for which asymptotic robustness analysis is more important for accurate results, as demonstrated above by the significant contributions from higher order quantum pathways (which are associated with higher order moments of noisy/uncertain parameters) in those cases. It is, however, possible to employ asymptotic and leading order Taylor approaches for quantum control robustness analysis together, in order to mitigate computational expense while retaining suitable accuracy.
In this section we consider leading order Taylor approximations for quantum control moments and their relationship to asymptotic robustness analysis methods. We consider robustness of the control performance measure (e.g., transition probability) to variations δθ in the parameters. Assuming the covariance matrix of parameter estimates is available, as in the posterior distribution of δθ modeled as a multivariate normal distribution, i.e., θ ∼ N (θ, Σ), through choice of a confidence level c, we can specify the set of possible realizations of δθ corresponding to that confidence level as: where χ 2 K denotes the chi-square distribution with K degrees of freedom, K denoting the number of noisy or uncertain parameters. The distribution of δθ can be used to estimate the corresponding distribution of the control performance measure J. Note that the choice of confidence level can also be used to set the maximum relevant pathway order in asymptotic robustness analysis, as discussed in the appendix.
Consider as an example the case of dipole operator uncertainty and let J = P ji , the transition probability between states i and j. A lst order Taylor expansion can provide only a normal distribution with variance var J ≈ Tr Σ∇ θ J(∇ θ J) T , where (setting = 1 for simplicity of notation) and X i is the Hermitian matrix obtained by setting The leading term in the expansion for E[δJ] due to parameter uncertainty is of 2nd order: where H θ, θ = d 2 J dθ dθ denotes the Hessian matrix with respect to Hamiltonian parameters. The elements of the Hessian matrix with respect to Hamiltonian parameter uncertainty can be written as follows in the most general case of the control of an arbitrary quantum observable O starting from any initial density matrix ρ 0 : From this expression, the rank of the Hessian at the optima of the quantum control landscape, as well as the eigenvectors corresponding to the downward directions, can be identified. In practice leading order approximations beyond first and second-order approximations [18,37] are almost never used, due to the need for expressions for higher order derivatives of J and the lack of available numerical methods for their evaluation beyond term by term calculation, in contrast to the asymptotic robustness analysis techniques introduced above.
For quantum systems that are not fully controllable [10], parameter uncertainty can present an opportunity to improve E [J]. To see this consider the overlap: where Λ 0, andṼ is an orthogonal matrix. If a quantum system is not fully controllable, it is possible that E[J] > J max nom where J max nom denotes the maximum possible value of J that can be achieved by control field optimization, given the estimated values of the time-independent Hamiltonian parameters. In this case-which can frequently arise in practice due, for example, to laboratory constraints on pulse shaping [12] or the structure of the Hamiltonian-it is possible to search for fields where, given the uncertainty spectrum, the overlap of the Hessian with the directions in parameter space associated with largest positive eigenvalues is maximized. This will be discussed in detail in a separate work.
The above leading order approximations are convenient in problems wherein robustness with respect to Hamiltonian parameter uncertainty is to be quantified and optimized. A first-order approximation to the variance of observable control fidelity due to field noise can be found in closed form: where Acf denotes the autocovariance function of the field noise process and We have the following expression for the leading order correction to E[δJ] due to control field noise: where Hess denotes the Hessian of the objective function with respect to control field variations (as will be discussed below, this can also be represented in the frequency domain, which is especially natural when the noise is in manipulated pulse shaping parameters). It is useful to consider the Hessian nullspace at the global maxima of the control landscape, which affects the robustness of optimal control fields to field noise. Hess t, t is finite rank kernel: which has rank 2N − 2, where N is Hilbert space dimension, for state-to-state population transfer [10,58] (μ(t) denotes the time-evolved dipole operator). Thus to second order, field noise at most frequencies does not affect population transfer. (Laser noise can only decrease E [J nom + δJ] since the Hessian is negative semidefinite). Note from the expressions above that leading order methods use only the first or second moments of the noisy or uncertain parameters in computing quantum control objective moments, whereas as shown in section 5.1 systems with higher Lie algebraic depth often exploit higher order quantum pathways-which require consideration of higher order moments of noisy or uncertain parameters-in achieving control fidelity.
Since leading order methods provide fast approximations to asymptotic moment calculations, the expense of asymptotic calculations can be reduced through the use of significant quantum pathways identified through leading order sensitivity analysis methods. Consider the problem of minimizing var(J) in the presence of field noise when either J nom has been maximized (the 'top' of the quantum control landscape), where J nom denotes the nominal value of J for a given control field and estimated value of the time-independent Hamiltonian parameters in the absence of uncertainty/noise. The first-order Taylor approximations for both E[J] and var(J) are both 0 at the optima of the nominal quantum control landscape and thus second order approximations must be applied at those optima. The leading order approximation for E [δJ] at the global extrema of the quantum control landscape [55][56][57][58] can (with the appropriate sign) be interpreted as a mean absolute error with respect to J nom , which is an alternative measure of variability. This is because at the optimum value of J nom , E[δJ] has a (semi)definite sign across all realizations of the noisy control field or the uncertain Hamiltonian parameter values 4 . Since the first-order contribution to E[δJ] is zero at J max nom , the leading order term is (assuming a stationary noise process; not required): where Acf denotes the two-point autocovariance function, acf denotes the autocovariance function for lag τ = t − t, F denotes the Fourier transform, P(ω) denotes the power spectral density of the stationary field noise process δε t , and where the last line is according to the Wiener-Khinchin theorem [37]. P(ω) is convenient for estimation of the autocovariance of the stochastic process δε t when a single long time sample is used for estimation. More generally, when repeated samples over dynamical time intervals [0, T] are used for estimation, acf(τ ) can be estimated directly, or the frequency domain representation of the stochastic process δε ω can be used directly [17] to approximate E[δJ] through the frequency domain autocovariance matrix Acf(ω , ω): This representation enables the identification of significant quantum pathways. As noted above, studies of quantum control landscapes have revealed that Hess is a finite-rank kernel of rank 2N − 2 for state to state population transfer on Hilbert space dimension N [10]. At the 'top' of the quantum control landscape, for var(J) calculation, the noisy/uncertain parameters corresponding to terms in the integral (5) above a threshold magnitude (of which there are at most 2N − 2 in the basis where the Hessian is diagonal) should be included in the asymptotic moment calculation. By defining the interaction Hamiltonian in terms of these noisy spectral perturbations to the nominal control field, it is easy to see that if the nominal field generates a propagator that maximizes J, this asymptotic moment correction in terms of significant pathways provides an estimate of E[δJ] with a consistent negative semidefinite sign, reduced computational expense compared to the full asymptotic calculation, and increased accuracy compared to leading order methods because all significant order products of these spectral parameters are included in the moment calculation (whereas the leading order correction is expressed in terms of the same basis functions but is limited to second order). Therefore at the 'top', transformation to a Fourier basis wherein Hess(ω, ω ) is diagonal, along with expression of the noise/uncertainty distribution on these transformed parameters, enables significant reduction in computational expense of the asymptotic var(J) (or E[δJ]) calculation, through these use of at most 2N − 2 encoding frequencies for population transfer problems (analogously for Hamiltonian parameter uncertainty). This provides a much more accurate estimate for var(J) than leading order techniques, while similarly omitting the contributions of parameters for which the sensitivity is lower. These conditions result in dimensionality reduction in the number of encoded field or Hamiltonian parameters needed for accurate asymptotic moment calculations.

Summary and prospective
We have presented the theory of robust control of quantum dynamics in the presence of noise/uncertainty in classical input variables (a phase-coherent quantized field interacting semiclassically with an ensemble quantum system) or the quantum Hamiltonian. This theory has been presented from the perspectives of both asymptotic and leading order robustness analysis, and methods for integrating these two methodologies have also been described. Robust optimal control theory has also been extended to the optimal control of quantum control objective moments in the presence of field noise. In addition, the engineering tools of frequency domain response theory have been extended to quantum systems.
Whereas we did not consider the problem of incoherence in field photons (where noise in phases implies lack of coherence among photons), which results in incoherent evolution because different individual quantum systems in an ensemble interact with photons with different phases, future work will study robust control in the presence of molecular interactions leading to dissipative (nonunitary) evolution. This choice of focus is due to the fact that molecular interactions leading to decoherence are more physically relevant than incoherence in the time-dependent control input in the ensemble control applications of greatest interest to chemical physics, where the primary applications of interest lie. Recent advances in the theory of open quantum systems should enable robust control of both Markovian and non-Markovian quantum chemical dynamics to be implemented using similar computational methods, due to the ability to cast non-Markovian quantum dynamics in terms of Markovian dynamics [59,60].
The field noise we considered in this paper is described by moments (means, covariances, higher order moments) of a discrete set of frequencies; although there are no restrictions as to the source of this noise, a common scenario is when the discrete set of field mode amplitudes and phases that are used to shape the control field are subject to noise. Future work should also consider the corresponding theory in the case where field noise is described by a continuous spectrum (stochastic noise processes in the time domain).

Data availability statement
No new data were created or analysed in this study.

Appendix A
• Expression for amplitude pathways where the shorthand notations A(ω k ) = A k and ω ji = (E j −E i ) have been used, and where the sum (k 1 ,...,k m) is over all 1 k i K, i = 1, . . . , m such that mode k appears in the multiple integral α k times.
• Amplitude pathways in terms of encoding/decoding Amplitude encoding: where γ k is the modulating frequency specific to the amplitude power α k associated with a particular pathway ( α). Using the modulation, the Schrödinger equation can be propagated in the time variable t and dummy variable s, for which the resulting encoded mth order transition amplitude U m ji (T, s) is comprised of the following terms: The encoded total transition amplitude contains the following terms: Deconvolution of the total transition amplitude leads to where γ ∈ [0, . . . , α 1 γ 1 + · · · + α k γ k ]. This suggests that all amplitude pathways of different orders can be extracted through deconvolution of the encoded transition amplitude if all γ s associated with each pathway is uniquely known, i.e. U ji (T, γ = α 1 γ 1 + · · · + α K γ K ) → U m ji (T, α). We can thus concisely define amplitude pathways α.
• Dipole pathways in terms of encoding/decoding. Dipole encoding would reveal the contribution of the dipole moments in the transition amplitude. Here, each of the dipole matrix elements is encoded with a Fourier function: The encoded and propagated unitary propagator consists of the different order dipole pathways with the encoded total transition amplitude:

Appendix B
In addition to the identification and encoding of significant pathways based on leading order sensitivity analysis, the computational expense of asymptotic robustness analysis can be reduced by (a) defining the interaction picture in terms of the deviations of field and Hamiltonian parameters from their nominal values, and encoding only these deviations; (b) dividing the dynamical time T into subdomains of durations for which the fast Fourier transform (FFT) can extract contributions of all relevant pathways, and encoding the dynamics separately on each of these time domains. For both (a) and (b), an expression for the bounds on the pathway order that will contribute to a quantum control objective moment is required. We derive these bounds below. Future work will describe how these bounds can be used in numerical implementation of asymptotic robustness analysis using methods (a) and (b). We note that these bounds are less relevant for systems of low saturating Lie algebraic depth for the control objective of interest.
The upper bound of pathway calculation is . . . where | j|μ|i | < d, and i, j ∈ [1, N]. Correspondingly, the upper bound of the Dyson term is given as: The upper bound on the error associated with the robustness calculations can in turn be determined: where m eps denotes the smallest m such that eps, eps denoting the smallest floating point number that can be represented on the computer • Compute m max and bound on error for specified • Note, it is possible to derive an analytical bound on the error, but it is less accurate and not necessary. • Compute time-domain integrals efficiently through FFT of μ(t) (N(N + 1)/2 FFTs of complex functions − N diagonal elements are real-valued)-Fourier transform provides both 1st and 2nd integrals above via μ(ω) at all frequencies ω • For a given J, only need to compute FFT of one scalar function of time.