Exploring the complexity of quantum control optimization trajectories

Abstract

The control of quantum system dynamics is generally performed by seeking a suitable applied field. The physical objective as a functional of the field forms the quantum control landscape, whose topology, under certain conditions, has been shown to contain no critical point suboptimal traps, thereby enabling effective searches for fields that give the global maximum of the objective. This paper addresses the structure of the landscape as a complement to topological critical point features. Recent work showed that landscape structure is highly favorable for optimization of state-to-state transition probabilities, in that gradient-based control trajectories to the global maximum value are nearly straight paths. The landscape structure is codified in the metric R ≥ 1.0, defined as the ratio of the length of the control trajectory to the Euclidean distance between the initial and optimal controls. A value of R = 1 would indicate an exactly straight trajectory to the optimal observable value. This paper extends the state-to-state transition probability results to the quantum ensemble and unitary transformation control landscapes. Again, nearly straight trajectories predominate, and we demonstrate that R can take values approaching 1.0 with high precision. However, the interplay of optimization trajectories with critical saddle submanifolds is found to influence landscape structure. A fundamental relationship necessary for perfectly straight gradient-based control trajectories is derived, wherein the gradient on the quantum control landscape must be an eigenfunction of the Hessian. This relation is an indicator of landscape structure and may provide a means to identify physical conditions when control trajectories can achieve perfect linearity. The collective favorable landscape topology and structure provide a foundation to understand why optimal quantum control can be readily achieved.

---

1 Introduction

Quantum optimal control theory (OCT) provides a basis to explore applied field tuning of quantum system dynamics.^1–5 Recent successes include achieving control of electrons in quantum dots through time-varying gate voltages,⁶ controlling electronic correlations in molecules,⁷ manipulating the dynamics of quantum many-body systems,⁸ cooling ultracold atomic gases,⁹ eliciting quantum revivals,¹⁰ and generating adiabatic time evolution.¹¹ Optimal control experiments (OCE) have found growing success even for complex systems, especially through the use of closed-loop learning algorithms¹² combined with femtosecond laser pulse shaping technology.¹³ For example, these advances have enabled vibrational control of population inversion in Bose–Einstein condensates,¹⁴ strong-field ionization of silver atoms,¹⁵ and coherent energy transfer in light-harvesting complexes.¹⁶ Underlying these achievements is the question of why desirable quantum controls can evidently be found with only modest search effort over the vast space of possible control fields. An important factor in addressing this question is the favorable control landscape topology, which results upon satisfaction of certain physical assumptions.¹⁷ The present paper considers a second contributing factor based on the landscape structure (i.e., non-topological features) reflected in the behavior of the evolving controls during the optimization process.

The dynamics of a closed quantum system interacting with an applied field is described by the time-dependent Schrödinger equation

where H(t) is the time-dependent Hamiltonian and U(t,0) is the corresponding evolution operator (propagator) with 0 ≤ t ≤ T. Here T is the time at which the objective is optimized. Treating the interaction within the dipole approximation, the time-dependent Hamiltonian is given by

H(t) = H₀ − μE(t), (2)

where H₀ is the field-free Hamiltonian, μ is the dipole moment operator and E(t) is the time-dependent control field.

Through eqn (1), a control field can be mapped to a value of a cost functional J[E(t)], which may be either the expectation value J_O of an observable O or the cost functional J_W for creating a target unitary transformation W (see eqn (3) and (4), respectively). The functional relationship between J[E(t)] and the control field E(t) can be cast in terms of a quantum control landscape, where each point J[E(t)] on the landscape corresponds to a particular control field E(t), 0 ≤ t ≤ T. Finding successive control fields that eventually lead to an optimal value of the cost functional entails climbing the landscape. We will refer to ‘climbing’ as moving towards the target value of J, regardless of whether the cost functional is to be maximized or minimized. Using a local search algorithm to climb (e.g., the gradient algorithm utilized in this work) results in smooth paths up the landscape J[E(s,t)] that can be parameterized by the continuous variable s ≥ 0. Every such path up the landscape is associated with a trajectory E(s,t) through control space. Thus, starting with an initial control field E(0,t) that produces a low value of the cost functional, the goal is to find the final control field E(s_max,t) that yields a high (optimal) value of the cost functional at s = s_max.

Many recent studies have considered the topology of quantum control landscapes.^17–21 The control landscape topology is specified by the character of the critical points where , ∀t∈ [0,T]. A critical point at a suboptimal value of J where the Hessian is negative semi-definite would indicate a (second order) trap capable of halting an effort attempting to reach the top of the landscape. The landscape topology has been analyzed for closed quantum systems with N states under the assumptions of: (i) controllability, such that any unitary evolution matrix U(T,0) can be generated by at least one control field at some sufficiently large time T, (ii) surjectivity, whereby the set of functions , ∀i, j are linearly independent over t∈ [0,T], and (iii) full accessibility, allowing for free access to any control field. Upon satisfaction of these assumptions, the quantum control landscape can be shown to contain only trap-free intermediate value critical saddle submanifolds (when they exist for a particular application)^17,18,21 along with the absolute minimum and absolute maximum submanifolds. Thus, the landscape topology is very favorable to performing OCEs or OCT simulations, which is an important factor in explaining the evident ease of locating optimal control fields.

Beyond issues of topology, the structure of the quantum control landscape away from the critical points could greatly influence the nature of a control trajectory winding its way to an optimal field. If landscape structure is highly gnarled, this circumstance could call for complex control trajectories, although the trap free topology assures that the path (e.g., guided by a gradient ascent algorithm) will lead to an optimal solution. The encountered landscape structure upon climbing may be examined by quantifying the linearity of the control trajectories.²² A predilection towards nearly straight trajectories is indicative of very simple control landscape structural features. Prior experimental work aiming to maximize second harmonic generation found that all of the sampled control trajectories were strikingly close to being straight.²⁶ Recent theoretical work has also supported the presence of simple landscape structure by examining the linearity of control trajectories in simulations over the state-to-state transition probability landscape.^22,27,28

In this work, we extend the state-to-state transition probability quantum control landscape structure analysis to consider two additional landscapes:

J_O = Tr(ρ(T)O) = Tr(U(T,0)ρ(0)U^†(T,0)O), (3)

and

J_W = ∥W − U(T,0)∥² = 2N − 2Re{Tr(W^†U(T,0))} (4)

where ∥·∥ denotes the Frobenius norm and N is the number of states in the quantum system. The quantum ensemble control landscape J_O allows for controlling arbitrary observables and consideration of systems which are initially in mixed states ρ(0), and the unitary transformation landscape J_W is relevant to the prospect of creating logic gates for quantum information processing.²⁹ Both landscapes can contain critical saddle submanifolds,^18,21 which are not present in the simpler state-to-state transition probability landscape. Thus, we aim to assess the influence of these saddle submanifolds on landscape structure. Specifically, we will show that encountering saddle submanifolds at intermediate heights on the landscape can drive control trajectories away from straight paths. This phenomenon is sketched in Fig. 1. The trajectory on the right takes a rather straight path up the landscape, while the trajectory on the left is initially attracted to the saddle point, thus altering its direct nature. Nevertheless, we will show in numerical simulations that the trajectories leading to optimal controls are generally still close to being straight, despite the presence of saddles. Lastly, we will derive a mathematical criterion showing that the gradient of the cost functional must be an eigenfunction of the Hessian in order to achieve a perfectly straight gradient-based control trajectory.

Fig. 1 Illustrative sketch of a quantum control landscape (blue surface) and its underlying control space (colored contour map), also displayed as the projection onto the plane spanned by two control variables (in practice many more are used). The vertical axis represents a cost functional J subject to maximization. The two optimization paths shown, in green, differ in the nature of their R values. The path on the right, proceeding from E₁(0,t) to E₁(s_max,t), climbs in a direct manner to the top of the control landscape. Consequently, its projection into control space, in black, has a path length very close to the Euclidean distance along the straight line between the projection's endpoints, shown as a magenta dashed line. Thus, this path's R value is low. On the other hand, the optimization path on the left, from E₂(0,t) to E₂(s_max,t), is attracted away from proceeding directly to a maximum by the saddle point at the center of the landscape. As a result, its projection into control space is much longer than the straight line between its endpoints, and so this path possesses a higher R value. The paths do not strictly terminate at the top and bottom of the landscape, as optimizations generally begin and end at values of the cost functional which approximate J^min and J^max, respectively, to acceptable precision.

The remainder of the paper is organized as follows. The linearity measure R of a control trajectory is defined in Section 2. Section 3 gives the gradient-based search algorithm used in this work to identify control fields that optimize arbitrary observables or achieve unitary transformations, respectively, through J_O or J_W. The presence of landscape saddle submanifolds is examined, and a metric is outlined that gives the distance from any point on a trajectory to a particular saddle submanifold. In Section 4, a derivation is presented specifying the Hessian-gradient eigen-relation for straight control trajectories. Section 5 provides numerical results on the statistical behavior of R over a large ensemble of landscape trajectories for J_O and J_W. We also illustrate satisfaction of the Hessian-gradient eigen-relation for nearly straight control trajectories. Concluding remarks are furnished in Section 6.

2 The control trajectory linearity measure R

The goal of the present work is to explore quantum control landscape structure reflected in the nature of the control trajectories followed during optimizations. Fig. 1 indicates that a landscape trajectory has a projected image in the underlying control space. As we use a myopic gradient algorithm to climb the landscape, the degree of gnarled character encountered on the climb is manifested in the trajectory through control space. To this end, we construct a ratio R delineating the linearity of the control trajectory. Specifically, R is the ratio of the path length of the control trajectory between the initial and final fields, given byto the Euclidean distance between the initial and final fields in control space, given bywhere s parameterizes the control trajectory. R is expressed as

The lower bound for the ratio is R = 1, which holds for a trajectory that takes a straight path between E(0,t) and E(s_max,t). A trajectory with R = 1 corresponds to the simplest possible path a gradient optimization can take to climb the quantum control landscape J[E(t)], and trajectories approaching this limit indicate a lack of encounters with complex structural features. High values of R mean that a control trajectory meanders through control space along a curve whose total length is much greater than the distance between its two endpoints, signifying that the optimization path encountered rugged structural features on the landscape.

In a recent work, we found values of R − 1 as low as ∼10⁻⁴ for a five level state-to-state transition probability landscape.²² By examining the distances between the endpoints of the control trajectories, we also observed that trajectories with low and high R values are distributed almost identically throughout control space. Furthermore, we found that it was possible to characterize the complexity of low R trajectories through a ‘straight shot’ climbing procedure. This procedure operates by calculating the gradient of the objective (e.g., the state-to-state transition probability) at the initial control field, and then proceeds to evolve the field along that specified direction until a local maximum is encountered. For suitably low R trajectories, nearly optimal control fields may be found in this way by merely continuing to march in the initially identified direction. This technique is useful for characterizing the straight nature of a trajectory, but it is still short of a constructive control algorithm without knowledge of a suitable initial field.

3 Negotiating the quantum control landscape

The quantum systems considered here have N-level Hamiltonians H in eqn (2) consisting of the field-free Hamiltonian H₀ as an N × N diagonal matrix, and μ as the dipole moment given by an N × N real symmetric matrix. We will utilize the gradient algorithm to prescribe the path of steepest ascent, thereby forming a natural “rule” to climb the landscape and consequently to navigate the underlying control space. The D-MORPH procedure will be used to implement the gradient search trajectory.^30,31 The process of climbing the landscape may be confounded by the presence of saddle submanifolds located at certain regions of specific intermediate heights on the landscape. The cartoon in Fig. 1 reflects this issue, and one goal of this work is to examine the degree to which saddles distort control trajectories.

We parameterize a trajectory by a continuous variable s ≥ 0, such that the control field becomes E(s,t). As s is increased, the change in the cost functional can be written using the chain rule as

Maximization of J stipulates that , which is assured by settingLikewise, when minimizing J, we require , so in that case we set . Minimization applies to J_W, and J_O may be minimized or maximized depending on the nature of the objective operator O. In the present work, the control variables adjusted by the gradient algorithm are the values of the control field E(s,t_i) at a set of discrete time points t_i, i = 1,2,…. In order to ensure that each simulation for a particular application traverses the same domain up (or down) the quantum control landscape, every trajectory begins at the same initial height on the landscape J^I and ends at the same final height on the landscape J^F. This criterion is met by using the D-MORPH algorithm^30,31 to first guide a climb down (up) the landscape to J^I, if an initial trial field yields a value of J that is above (below) J^I. The resultant field at J^I is taken as E(0,t) for the purpose of determining R. The final height on the landscape J^F is then achieved by E(s = s_max,t), identified from application of D-MORPH starting with the initial field E(0,t).

Calculation of the gradient requires utilization of the relevant cost functional, either J_O or J_W. Sections 3.1 and 3.2 summarize these details, further include the characterization of the critical submanifolds of the landscapes, and introduce a method to determine the distance of a point on the landscape to any particular critical submanifold. In all of the simulations, the assumptions underlying the favorable landscape topology are taken as satisfied;^17,18,21 consistent with these assumptions, no landscape trajectories prematurely terminated before reaching J^F in our simulations.

3.1 Quantum ensemble control landscape

The quantum ensemble control landscape J_O[E(s,t)] for optimizing an arbitrary observable with a mixed initial state is given by eqn (3). The gradient of the objective may be written as²¹where we have made use of the relation in the second equality. Climbing the landscape is affected by utilization of eqn (10) in eqn (9). This process could possibly be hindered if the trajectory comes near a saddle submanifold where is small at an intermediate, suboptimal value of J_O. A representative sketch is given in Fig. 1, where the optimization path on the left draws near a saddle submanifold, located at the center of the blue landscape surface. In order to assess the impact of encountering saddle manifolds upon the value of R for a control trajectory, we first briefly summarize the nature of the critical submanifolds of the J_O landscape.

At a critical point, the first order variation satisfies . Thus, examining the last line of eqn (10), we see that if is of full rank (i.e., the surjectivity assumption (ii) from Section 1 is satisfied), we must have [ρ(T),O] = 0 at a critical point.²⁴ From this relation, it is possible to obtain the form of the critical matrices U_o, and from that information, the distances D to a particular critical submanifold from any point on the landscape; Appendix A contains further details. The distance metric can be used to judge if a control trajectory has strayed close to the saddle submanifolds, and therefore D allows us to assess whether proximity to these manifolds affects the R value of a trajectory.

3.2 Unitary transformation control landscape

In quantum information science, implementing basic logic gate operations necessitates the control of unitary transformations.²⁹ We therefore consider the optimization of the entire propagator U(T,0), seeking to attain W, which leads to the unitary transformation cost functional given by eqn (4). J_W has a minimum value of 0, corresponding to U = W, and a maximum of 4N, corresponding to U = −W, where N is the number of states in the quantum system. Thus, we would like to minimize J_W.

The gradient for the unitary transformation landscape is²

Optimizations may be carried out on the unitary transformation landscape by combining eqn (9) and (11). Trajectories climbing the J_W landscape may pass by saddle submanifolds. In order to assess the effect of these manifolds on R, we briefly review the nature of the critical points of the unitary transformation landscape.

From the last line of eqn (11), we see that at a critical point, where , if is of full rank, we must have W^†U(T,0) = U^†(T,0)W. This implies that W^†U_o is Hermitian. W^†U_o is also a unitary matrix, so at a critical point, the eigenvalues of W^†U_o must be either 1 or −1. As a result, the critical values of J_W in eqn (4) are 0,4,…,4N, implying that there exist N + 1 critical submanifolds, each of which correspond to a different number of 1's and −1's as the eigenvalues of W^†U_o.^18,19 A discussion of the distance metric D between any point on the landscape and a particular critical submanifold is contained in Appendix B. As explained in Section 3.1, evaluating D along a landscape climb to reach W can reveal the influence of encountering saddles on the path to an optimal control solution.

4 ‘Straight shot’ eigen-relation

Section 3 laid out how to climb the quantum ensemble and unitary transformation control landscapes using a gradient algorithm, and summarized the topology of the critical points of these landscapes. Here, we present a relation which must be satisfied for ‘straight shot’ control trajectories, utilizing the gradient algorithm for any quantum control landscape.

Consider a control trajectory that is perfectly straight, i.e., R = 1. Any such ‘straight shot’ trajectory can be written²²

E(s,t) = E(0,t) + ρ(s)ΔE(t), (12)

where 0 ≤ ρ(s) ≤ 1 is a monotonically increasing function of s as a result of the gradient algorithm assuring that dJ/ds ≥ 0 (i.e., without loss of generality we seek to maximize J), and ΔE(t) = E(s_max,t) − E(0,t) depends only on time. Utilizing eqn (9), the gradient function along this trajectory iswhere .

According to eqn (9), the instantaneous tangent of a control trajectory, , “points” in the same direction in control space as the gradient of the quantum control landscape, , at location s. Furthermore, along a straight control trajectory, the gradient function must point in the same direction at every location on the trajectory; thus, the gradient functions and at two different locations s and s′ on a straight control trajectory are related by a scale factor. Additionally, for straight trajectories, the difference between any two fields on the trajectory, E(s,t) and E(s′,t) is always proportional to the same function ΔE(t). When using a gradient-based climbing algorithm, the function ΔE(t) is proportional to the gradient itself, as shown by eqn (13). These facts allow us to conclude that an eigen-relation exists between the Hessian and the gradient on a straight control trajectory.^32,33

To make these considerations more specific, we expand the gradient of the cost functional J[E(s + ds,t)] to first order in ds along an arbitrary control trajectory at location s

Here, is the Hessian. Now assume that R = 1 for the control trajectory E(s,t). If we bring onto the left hand side, divide by ds, and take the limit, we havewhere eqn (13) was used for the second equality. Using eqn (13), we can rewrite the result as . On the right hand side of eqn (14), if we use eqn (9), we obtainor more compactly, with the time dependence implicitly understood,where g(s) is the gradient of the cost functional, , and the scalar product denotes the integral operation in eqn (16). As the action of the Hessian on the gradient merely rescales the gradient, this equation states that the gradient is an eigenfunction of the Hessian along a straight control trajectory, with an associated eigenvalue . We will illustrate this relationship in Section 5.3. In addition, one may show that repeating the above procedure by expanding the Hessian about the point s and making use of eqn (12) yields a higher order eigen-relation, stating that the gradient is an eigenfunction of the third-order Hessian. Iterating this process yields a hierarchy of eigen-relations between the gradient and all higher-order derivatives of J; however, all of these relationships contain no more information than eqn (17), since the entire hierarchy holds upon satisfaction of eqn (9) and (12).

Near a critical point, the eigen-relation is always true, even for trajectories with R > 1. To see this relation, expand the gradient of J[E(s₀ + ds,t)] to first order about the location s₀, which is infinitesimally close to a critical point at E(s_c,t)

where ds = s_c − s₀. Utilizing the gradient rule in eqn (9) inside the integral and noting that at a critical point, we have

Eqn (19) shows that the gradient function is always an eigenfunction of the Hessian in the immediate neighborhood of a critical point, whether or not the global trajectory is a straight shot. When R = 1, however, this behavior extends to the entire control trajectory. From the definition of ds, one can see that if s_c < s₀, so that the gradient trajectory is leaving the critical manifold, the eigenvalue in eqn (19) is positive, and likewise if s_c > s₀, so that the trajectory is arriving at the critical manifold, the eigenvalue is negative. This circumstance will be shown in Section 5.3 for the case of a trajectory moving away from the critical point at the bottom of the landscape where ds < 0 and then transitioning to the neighborhood of the top of the landscape where ds > 0. Thus, along a straight shot from the bottom to the top of the landscape, the Hessian eigenvalue associated with the gradient will change sign (i.e., at least once, and more often if saddle critical points are closely approached) at some point along the trajectory. In addition, the eigen-relation is nominally an independent consideration at each value of s. However, the diffeomorphic nature of the gradient climb over s in eqn (9) implies that a smooth eigen-relation should exist, regardless of which eigenvector (corresponding to the gradient function) initiates the climb. This behavior is affirmed in the simulations in Section 5.3.

We note that the relations examined in this section were derived without making use of quantum mechanics, as the analysis leading to the gradient being an eigenfunction of the Hessian only utilizes the gradient climbing rule, eqn (9), and the straight shot condition eqn (12). The interest in this paper is directed towards applications where J[E(s,t)] arises in a quantum control context.

5 Illustrations

This section provides numerical illustrations of the landscape structure analysis tools outlined in Sections 2–4. Section 5.1 shows that for randomly chosen control trajectories, the distribution of R takes on very modest values, which did not exceed 2.0 in the present examples. The results also demonstrate that control trajectories with either very low or very high values of R are spread in a like fashion throughout control space. Section 5.2 addresses the effect of saddle submanifolds upon optimization paths. Encountering saddle critical submanifolds on the way from the bottom to the top along an optimization path tends to modestly increase the associated R value. Finally, Section 5.3 provides numerical evidence for the validity of the Hessian-gradient eigenrelation in Section 4. Dimensionless units are used throughout the simulations.

5.1 Random sampling of control trajectories

In order to assess the structural complexity of the quantum ensemble and unitary transformation control landscapes, we recorded the R values of landscape trajectories constructed using randomly chosen initial control fields.

5.1.1 Quantum ensemble control landscape. Here we consider several eight-level quantum systems to illustrate the landscape structural findings. As the nature of a landscape's saddle manifolds is determined by the number and degree of degeneracies of the initial density matrix and observable matrix,^21,23,25 we will consider combinations of matrices with varying degrees of degeneracy so as to examine a broad class of landscapes. Increasing numbers of nondegenerate eigenvalues generally correspond to larger numbers of saddle submanifolds. In particular, we are interested in the effects of such submanifolds upon landscape structure. Accordingly, five different quantum ensemble control landscapes are sampled using the ρ and O matrices presented in Table 1, and 1000 optimizations of J_O are performed over each landscape. The Hamiltonian in Table 1 is utilized together with the dipole matrixwhere the signs on the dipole matrix elements were chosen randomly under the constraint of μ remaining symmetric. Once these signs were chosen, they remained fixed for the 1000 optimizations over each individual landscape.

Table 1 Quantum ensemble landscape matrices

Matrix	Diagonal entries
ρ 1	(1,0,0,0,0,0,0,0)
ρ 2
ρ 3
O 1
O 2
O 3
H 0	(−10, −8, −4, 2, 10, 20, 32, 46)

---

The density matrix ρ₁, which corresponds to the case of a pure initial state, represents an almost completely degenerate density matrix. In contrast, ρ₂ contains two eigenvalues which are both four-fold degenerate. The observables used in conjunction with these density matrices O₁ and O₂. As these observable matrices also exhibit varying levels of degeneracy, the pairs ρ₁ and O₁, ρ₁ and O₂, ρ₂ and O₁, ρ₂ and O₂ are used to generate four distinct landscapes. Lastly, we also utilized ρ₃ and O₃ as an example of matrices with fully nondegenerate spectra to generate a fifth landscape.

The initial fields of these optimizations were parameterized as

with the target time T being 10. The amplitudes a_n and phases ϕ_n were chosen randomly from the uniform distributions [0,1] and [0,2π], respectively. The frequencies are ω_n = n and M = 60 to ensure that the initial control field was nearly resonant with every possible transition in the Hamiltonian in Table 1. F is a normalization factor, picked so that the initial fluence of each field is 1. This normalization was chosen to avoid operating with strong initial fields. After the initial choice of the control field, the time interval [0,T] was discretized into 1001 evenly spaced time points, and the value of the field at each of these time points served as the control variables. We define the limits of the optimizations on the landscape as follows. We denote J^max_O as the maximum value that J_O = Tr(ρ(T)O) can attain for a particular pair of ρ(0) and O, and J^min_O denotes the minimum value. When a random field drawn from eqn (21) is used to calculate J_O = Tr(ρ(T)O), if J_O < J_O^I = J^min_O + 0.01[J^max_O − J^min_O], then D-MORPH is used to climb until J_O = J^I_O. Likewise, if J_O > J^I_O, D-MORPH is used to descend the landscape until J_O = J^I_O. Then, the resultant control field is optimized using D-MORPH until J_O = J^F_O =J^max_O − 0.01[J^max_O − J^min_O]. This latter trajectory, which takes the field from a height of J^I_O to a height of J^F_O on the landscape, is used to calculate R. The large window from J^I_O→J^F_O is chosen so that each of the optimizations proceeds over a wide swathe of the landscape and thereby encounters as many structural features as possible. In this regard, we have numerically shown that for the state-to-state transition probability landscape, pushing the fidelity of J closer to either 0 or 1 at the beginning and end of the optimizations, respectively, has a diminutive effect on R.²²

Histograms of the resulting R values from 1000 optimizations over the five different landscapes are presented in Fig. 2 for the cases (a) ρ₁ and O₁, (b) ρ₁ and O₂, (c) ρ₂ and O₁, (d) ρ₂ and O₂, and (e) ρ₃ and O₃. The distributions are skewed towards the right, indicating that the R values tend to accumulate near 1.0. The smallest values are found in the histogram of Fig. 2(a) which has a mean of R = 1.22, followed by Fig. 2(b) and (c), which have mean values of R = 1.31 and R = 1.33, respectively. The highest values of R are found in the histograms of Fig. 2(d) and (e), which have mean values of R = 1.65 and R = 1.53, respectively. The distinctions between the encountered structural features are modest over all of the cases examined. In an absolute sense, all of the trajectories produced R < 2.0 indicating only very mild landscape structural features (considering that there is no a priori reason for R to be close to 1). Prior experience with state-to-state transition probability landscapes suggests that ‘straight shot’ trajectories are reserved for cases with R < 1.05, and a search is generally needed to find a proper E(0,t) to initiate such a trajectory.²²

Fig. 2 The distribution of R values for five different quantum ensemble control landscapes. The cases are generated using: (a) ρ₁ and O₁; (b) ρ₁ and O₂; (c) ρ₂ and O₁; (d) ρ₂ and O₂; and (e) ρ₃ and O₃. For each of these runs, the Hamiltonian H₀ of Table 1 and dipole μ of eqn (20) were used. All of the distributions lean towards the left, showing that R tends to accumulate near 1.0. 1000 optimizations were carried out for each landscape.

A natural question is whether control trajectories that possess low or high R values are concentrated in some particular region(s) of control space. As the initial controls in the statistical sampling considered above were randomly chosen, this assessment directly bears on the investigation of the range of R values encountered. To address this question, Euclidean distances were calculated between (1) all pairs of initial fields, (2) all pairs of final fields, and (3) all initial and final fields from the collection of control trajectories for the landscape associated with ρ₁ and O₁. These distances were calculated using

for control fields E_i and E_j, and the results are displayed in Fig. 3. Importantly, in computing these distances, we first used only the 250 trajectories with the lowest values of R out of the 1000 trajectories that were sampled in generating Fig. 2(a). Then, separately, we used the 250 trajectories with the highest values of R. In particular, trajectories with R < 1.173 were used to generate the histograms in Fig. 3(a), and likewise with R > 1.256 for the histograms in Fig. 3(b); in both cases the histograms show all possible pairs of distances, with 250 × 249/2 ∼ 3.1 × 10⁴ pairs. The distribution of distances between initial fields, in blue, is centered to the left and narrower than the distribution corresponding to the distances between optimal fields, in red. This behavior implies that the optimal fields are scattered more widely throughout control space than the randomly picked initial fields. The histogram in green of the distances between all initial-final field pairs is centered between the other two distributions. The most significant feature in Fig. 3 is that the distributions are very similar, whether generated from trajectories with low values of R (Fig. 3(a)) or from trajectories with high values of R (Fig. 3(b)). This result suggests that trajectories with low and high values of R are distributed in the same way across control space. Thus, straighter, more direct trajectories and those less so are likely to be equidistant from a trajectory originating from a random initial control field. Similar results have been observed in other cases with different ρ and O matrices.

Fig. 3 The distributions of pairwise distances between the initial fields, final fields, and initial-final fields for the landscape generated by ρ₁ and O₁ from Table 1, respectively. The distances are calculated with eqn (22). Plot (a) was produced using the 250 optimizations with the lowest values of R, while plot (b) was produced using the 250 optimizations that had the highest values of R. The similarity between the two plots suggests that low R and high R trajectories are distributed in essentially the same way throughout control space.

5.1.2 Unitary transformation control landscape. Here we consider a four-level quantum system for the target unitary transformationwith the Hamiltonian and dipole matrices beingwhere the random signs on the dipole matrix were again chosen so that it remained symmetric. A total of 1000 initial fields were randomly chosen according to eqn (21), with T = 10 and M set to 20. The initial fields were adjusted as explained in Section 5.1.1 to arrive at a value of J^I_W = 0.99 × 16 = 15.84, and the final achieved value after optimization was J^F_W = 0.01 × 16 = 0.16. Thus, in this case we climb down the landscape as a minimization process. Fig. 4 is a histogram of the R values collected in this fashion. The plot is skewed right with a mean of R = 1.41, indicating that the unitary transformation control landscape contains minimally complex structural features. The quantum ensemble landscapes considered earlier were for eight level systems, while this unitary transformation landscape is for a four level system. Together with previous work which observed that typical R values for generating unitary transformations scales linearly with the number of states N,²⁷ these findings suggest that the unitary transformation landscape appears to be structurally similar to the ensemble control landscape with fully nondegenerate ρ and O matrices.

Fig. 4 The distribution of R values for a unitary transformation control landscape. 1000 optimizations were carried out, using the target unitary transformation W of eqn (23) and Hamiltonian H₀ and dipole μ matrices of eqn (24). Modest values of R, all less than 2.0, characterize this landscape's structure as simple.

Fig. 5 displays the histograms resulting from calculating the pairwise distances between the control fields using eqn (22), as explained in Section 5.1.1. The two sets of control trajectories, each of size 250, were used to produce the histograms in Fig. 5(a) and (b) which satisfy R < 1.34 and R > 1.48, respectively. Significantly, the histograms computed with the lowest values of R and with the highest values of R are almost identical, indicating that trajectories with varying values of R for generating unitary transformations are distributed in the same way across control space, as found in Fig. 3 for the ensemble control landscape. This result was also found for a five level state-to-state transition probability landscape,²² suggesting that similar distributions of fields for high and low R trajectories is a general property of quantum control landscapes.

Fig. 5 The distributions of pairwise distances between initial fields, final fields, and initial-final fields for the unitary transformation control landscape, calculated using eqn (22). Plot (a) was produced using the 250 optimizations with the lowest values of R, while plot (b) was produced using the 250 optimizations that had the highest values of R. The similarity between the sets of distributions in (a) and (b) suggests that low R and high R trajectories, as well as all of the associated fields, are distributed in the same way throughout control space.

5.2 Interactions of optimization trajectories with saddle submanifolds

Although the saddle submanifolds of the quantum ensemble and unitary transformation landscapes are not able to prevent reaching the global optimum, they may still significantly influence the optimization trajectories. In order to examine the extent to which saddle submanifolds affect the climb up the landscape, the metrics set out in Sections 3.1.3 and 3.2.3 (and Appendices A and B, respectively) were used to calculate the distance from each point along a landscape climb trajectory to all critical submanifolds. As examples, the quantum ensemble control landscape was generated by the initial density matrix ρ₂ and observable operator O₁ from Table 1, while the unitary control landscape was generated with W in eqn (23).

5.2.1 Quantum ensemble control landscape. The values of the critical submanifolds of the system on the quantum ensemble control landscape²¹ of ρ₂ and O₁ are given by . Here, J¹_O and J²_O correspond to saddles, although a trajectory passing through these values may not actually come close to the associated submanifold of saddle controls. The normalized distances D to these four critical submanifolds, using eqn (28), are computed during the optimizations with the lowest and highest R values out of the 1000 runs performed, which were R = 1.17 and R = 1.59, respectively. The results are shown in Fig. 6(a) and (b). The magnitude of the gradient for each trajectory, scaled down by a factor of 60, is also plotted to show the rate of climbing, as an optimization typically slows down near a critical submanifold. In Fig. 6(a), the run with a small R value does not closely approach any of the saddle submanifolds. This observation is corroborated by examining the magnitude of the gradient, which falls only slightly in the middle of the climb. The gradient is naturally small near the beginning and end of the climb, where the trajectory leaves (approaches) the critical submanifolds corresponding to the bottom (top) of the landscape. In contrast, Fig. 6(b) shows that the run with a high R value hovers in the vicinity of the two saddle submanifolds for extended periods during the climb, indicated by the distance to the first and second saddle manifolds, D₁ and D₂, respectively, becoming very small during portions of the run. Concurrently, the magnitude of the gradient falls sharply in those two areas, recovering only when D₂ starts to rise, thus confirming that the run passed closely by these two saddle submanifolds. Comparison of Fig. 6(a) and (b) suggests that control trajectories that are ‘diverted’ to pass closely by saddle submanifolds can have higher R values.

Fig. 6 In (a) and (b), the distances D to the critical submanifolds (in arbitrary units) are plotted using eqn (28) along a low and high R value optimization trajectory, respectively, on a quantum ensemble control landscape. It can be seen in (b) that the distance to saddle submanifolds 1 and 2 become very small along the high R trajectory, indicating that the saddles distorted the trajectory by acting as attractors. The norm of the gradient function along the climb, scaled down by a factor of 60 and shown as the thick green curve, also noticeably dips in (b), further confirming that saddle submanifolds were encountered. For a unitary transformation control landscape in (c) and (d), the distances D to the critical submanifolds, using eqn (29), are plotted along a low and high R value optimization trajectory. In (c), no saddle submanifolds are encountered, but in (d) the black and green curves (i.e., for saddles 1 and 3, respectively) dip sharply along with the norm of the gradient (scaled down by a factor of 20), signaling attraction towards these saddle submanifolds and resulting in the corresponding high value of R.

5.2.2 Unitary transformation control landscape. Fig. 6(c) and (d) show the distances to the critical submanifolds while traversing the unitary transformation control landscape using eqn (29); the magnitude of the gradient (scaled down by a factor of 20) is also shown. D₁, D₂, and D₃ denote the normalized distances to the saddles at landscape heights of J¹_W = 12, J²_W = 8, and J³_W = 4, respectively, while D_min and D_max denote the normalized distances to the bottom and top critical submanifolds at J^min_W = 16 and J^max_W = 0, respectively. Fig. 6(c) corresponds to a run with a low R value of 1.15, and Fig. 6(d) corresponds to a run with a high R value of 1.73. The distances to the saddle submanifolds D₁ and D₃ assume lower values during the high R run (as compared to the low R run), and correspondingly the magnitude of the gradient also dips when the distances to D₁ and D₃ are at a minimum. Neither of the trajectories approach the second saddle submanifold, as indicated by the behavior of D₂ on the trajectory. These plots also support the notion that saddle submanifolds are capable of attracting nearby paths of steepest descent away from proceeding in a direct manner to a global minimum.

5.3 Hessian-gradient eigen-relation

The results of Section 4 identify that the gradient of the cost functional should be an eigenfunction of the Hessian along a straight trajectory. Here, we provide evidence to verify the Hessian-gradient eigen-relation for nearly straight shots found in simulations. In order to illustrate this phenomenon as clearly as possible, we consider a simple three level quantum system with Hamiltonian and dipole matrixand we use the cost functional of eqn (3), withcorresponding to maximizing the pure state transition probability P_1→3(T) with T = 10. This landscape contains no saddle submanifolds, and therefore the control trajectories tend to have smaller values of R. We used the stochastic Particle Swarm Optimization (PSO) algorithm^34,35 to locate an initial control E(0,t) taking the form in eqn (21) with M = 20 that produced a nearly straight trajectory with R = 1.002. In order to assess whether a lower bound for R (other than 1.0) might exist, a Derandomized Evolution Strategy (ES)^36,37 was used in a more extensive search for a straight path. Using this algorithm together with an n-point Gauss-Legendre quadrature rule (with n = 1000) to accurately evaluate the integrals in eqn (5) and (6) allowed reaching an even lower value of R = 1.00003. As this result was at the limit of our numerical accuracy, we cannot rule out the prospect of finding a still lower value of R. Below, the eigen-relation is tested with the former conservative near straight trajectory of R = 1.002.

We seek to affirm eqn (17), which states that the gradient function aligns with some eigenfunction of the Hessian. As in eqn (17), we simplify notation with the Hessian and gradient given by and g(s), respectively, and a ‘dot’ · denotes integration over time t ∈ [0,T]. Thus, we calculated , where is the function resulting from the action of the Hessian on the unit gradient function, . In Fig. 7, is plotted versus s as the thick black circles, along with every eigenvalue of the Hessian shown as the smaller green circles. Intriguingly, at every point on the climb, the black circle overlaps almost exactly with one of the green circles, suggesting that the gradient function can be identified with some eigenfunction of the Hessian. The sign of reveals that and g are parallel at the beginning of the climb and antiparallel at the end, which is consistent with the assessment discussed below eqn (19).

Fig. 7 The thick black circles are values of along a nearly straight shot to the top of a state-to-state transition probability landscape. Here is the Hessian and ĝ is the normalized gradient function. The eigenvalues of the Hessian are shown as the smaller green circles, for each point s along the climb of the landscape. The black circles match with a green circle at each point on the climb, supporting the identification of the gradient with an eigenfunction of the Hessian. As the associated eigenvalue of the Hessian changes from positive to negative, the function changes from being parallel to antiparallel to the gradient function.

We will now illustrate the validity of eqn (17), rewritten in the notation here as . This assessment is performed by replotting versus s in Fig. 8 along with ρ′′(s)/ρ′(s), calculated for this control trajectory. According to eqn (13), for a R = 1 straight control trajectory, . Since the trajectory considered here is not exactly straight with R = 1.002, we determined the average value of ρ′(s) by setting , where ΔE(t) = E(s_max,t) − E(0,t) and ρ′′(s) was then obtained by numerically differentiating ρ′(s). Examining Fig. 8 reveals that to good accuracy, providing evidence for the validity of the eigen-relation derived in Section 4.

Fig. 8 A plot of , which is the factor by which the Hessian scales the normalized gradient ĝ. From eqn (17), along a straight control trajectory we should have . Thus, the plot also shows ρ′′(s)/ρ′(s) evaluated along the gradient landscape climb. The good agreement between the two functions signifies the validity of eqn (17) along a straight control trajectory.

6 Conclusions

The exploration of landscape structure first performed for the state-to-state transition probability²² is extended in this paper to more general quantum ensemble and unitary transformation control landscapes. These collective results point towards the existence of simple quantum control landscape structure, which was first encountered experimentally²⁶ as well as hinted at by preliminary numerical studies.^27,28

A special feature of the landscapes considered in this work is the presence of critical submanifolds that possess saddle topology. These submanifolds do not prevent optimization trajectories from reaching the top of the landscape, but they can have an impact on the search effort, reflected in the identified R value. We find that the presence of saddle submanifolds on the landscape has a small but discernible effect on the ease of optimization reflected in the distortion of the control trajectories. Importantly, the appearance of saddle submanifolds did not cause R to take on values in excess of 2.0 for the systems examined here, indicating that these quantum control landscapes, in addition to state-to-state transition probability landscapes,²² possess markedly simple structure.

A condition for the existence of control trajectories with R = 1.0, which would indicate exceedingly simple landscape structure, was also derived. It was found that in order for the gradient algorithm to yield a straight trajectory, the gradient function must not only be separable in s and t, but must also be an eigenfunction of the Hessian. This situation also implies that all higher-order derivatives of the cost functional possess a common eigenfunction. The eigen-relation may form a basis to identify the conditions under which exactly straight control trajectories derived from the gradient algorithm exist. Additional simulations exploring landscape structure on a wider variety of landscapes would also allow for a better understanding of the conditions that lead to low values of R. One direction for further investigation is whether achieving control of unitary transformations, which requires control over every matrix element in the propagator, leads to higher values of R than even for the fully nondegenerate ensemble control landscapes. A simple assessment of this matter was performed for a two level system using an ES algorithm^36,37 to search for straight trajectories. For two-level systems, in the cases of maximizing P_1→2 we find R − 1 ≤ 10⁻⁶, while for unitary transformation control R = 1.0016 was achieved.

Quantum control landscapes are expected to nominally be very complex, as they are defined on extremely high dimensional spaces and represent the generally nonlinear relationship, encoded through Schrödinger's equation, between the unitary propagator and a dynamical set of controls. Thus, a natural expectation is that trajectories traveling over these landscapes would require considerable effort to reach an optimal outcome. However, the topology of quantum control landscapes is very simple when certain reasonable physical assumptions are met, providing one factor to explain the relative ease of performing optimizations.^4,17,38,39 These topological results are complemented by the finding in this work of simple landscape structure. This study, along with ref. 22, supports the notion that simple structure is a universal and robust property of quantum control landscapes. Collectively, the twofold simplicity of landscape features, encapsulated by topology and structure, provides a foundation for more clearly explaining the widespread efficiency of finding optimal controls over quantum phenomena.

Appendix A critical submanifolds of the quantum ensemble control landscape

This Appendix gives a brief overview of the critical submanifolds of the quantum ensemble control landscape, along with a metric which can be used to determine the distance to all critical submanifolds from any given point in control space. A complete characterization of the critical submanifolds is given in ref. 23.

The analysis of Section 3.1 reveals that at a critical point, [O,ρ(T)] = 0. This condition implies that the number of critical submanifolds is in one-to-one correspondence with the number of ways the eigenvalues of ρ = ρ(0) can be permuted such that their overlaps with the eigenvalues of O are distinct. These permutations are induced by critical propagators U(T,0), each of which belongs to a distinct critical submanifold, and can be enumerated using contingency tables.²³ These tables possess n rows and m columns and can be constructed as follows. First, ρ and O are expressed in diagonal form with their eigenvalues in descending order. Then, a particular permutation is applied to the eigenvalues of ρ; such a permutation matrix U_o is associated with the said critical submanifold. The number of times the eigenvalue λ_i of ρ overlaps with the eigenvalue ε_j of O, after the permutation is applied, is the value c_ij of the (i,j) entry in the contingency table. The rows of the contingency table must add up to the sum of the degeneracies of the eigenvalues of O, and the columns of the table must add up to the sum of the degeneracies of the eigenvalues of ρ. This fact can be used to easily construct all possible contingency tables and therefore reveal every critical submanifold along with its associated height on the landscape.

The contingency tables can also be used to calculate a distance to a particular critical submanifold from any point on the landscape. The distance is expressed in terms of the unitary transformation U[E(t)] corresponding to a particular control field E(t). A short exposition of this metric is given here; a full discussion is given in ref. 40. Consider a point on the landscape generated by U_o belonging to some critical submanifold; i.e., U_o is a permutation matrix. The matrix U_o may be divided into blocks {U^ij_o}, starting in the upper left corner, of sizes {o_i × p_j}, where o_i is the degeneracy of the ith eigenvalue of O, and p_j is the degeneracy of the jth eigenvalue of ρ(0), with the eigenvalues arranged in descending order. It can be shown that the singular value decomposition of each block contains only 0s and 1s, and that for the ijth block of size o_i × p_j the number of 1s is equal to c_ij, which is the entry in the ith row and jth column in the contingency table. This fact can be used to specify the distance to each critical submanifold. For each block and its associated entry c_ij, define an N × N diagonal matrix S_ij, which has c_ij values of 1 starting in the upper left corner in the diagonal:

To find the distance of an arbitrary unitary matrix U from the critical submanifold U_o, we divide U into blocks {U^ij} in the manner described above. Each block U^ij of U is written in its singular form Σ_ij, with the singular values on the diagonal sorted in descending order. For each block in U, we define V_ij = S_ij − Σ_ij. The distance from U to the critical submanifold of U_o is then

where p and o are the number of distinct eigenvalues of ρ(0) and O, respectively. It is convenient to normalize this metric.⁴¹ To do so, the unitary matrix U can be assigned a vector consisting of the singular values of all of its blocks {U^ij}. Since the squares of these singular values sum to N, each vector can be thought of as lying in the spherical polygon spanned by the vectors corresponding to the U_o matrices of the critical submanifolds. This formulation implies that the point which is the greatest distance from a critical submanifold is another vertex, which itself corresponds to another critical submanifold. By calculating the pairwise distances between all the critical submanifolds, we can find the maximum distance any point can be from a particular critical submanifold, and then normalize the metric by dividing all distances to that critical submanifold by the maximum distance, thereby giving 0 ≤ D ≤ 1. This normalization is used in the results reported in Fig. 6(a) and (b).

Appendix B critical submanifolds of the unitary transformation control landscape

Here we outline a distance metric for the unitary transformation control landscape; full details are given in ref. 40. Recall, from Section 3.2 that for a unitary transformation U_o corresponding to a critical point, the eigenvalues of W^†U_o are ±1. The distance from an arbitrary unitary matrix U to a submanifold containing the critical point U_o can be specified by noting that the eigenvalues of W^†U take values on the unit circle in the complex plane. The ith eigenvalue of W^†U is λ_i, and the eigenvalues are sorted in ascending order by their real parts. To find the distance to the critical submanifold which corresponds to α eigenvalues of −1 and N − α eigenvalues of 1 for W^†U_o, we compare the real part of the smaller α eigenvalues of W^†U with −1, and the real part of the larger N − α eigenvalues with 1.⁴⁰ The distance between U and the critical submanifold containing U_o then isThe maximum distance to a critical submanifold in this case is 2N, and the denominator of eqn (29) assures that 0 ≤ D ≤ 1.

Acknowledgements

We thank Greg Riviello and Julia Yan for helpful discussions. A.N. acknowledges support from the Program in Plasma Science and Technology at Princeton University. We also acknowledge support from ARO Grant No. W911NF-13-1-0237, ARO-MURI Grant No. W911NF-11-1-2068, DOE Grant No. DE-FG02-02ER15344, and NSF Grant No. CHE-1058644.

References

1. G. G. Balint-Kurti, S. Zou and A. Brown, Advances in Chemical Physics, Wiley, 2008, p. 43.
2. T.-S. Ho, J. Dominy and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 79, 013422.
3. D. Dong and I. R. Peterson, IEE Proc.: Control Theory Appl., 2010, 4, 2651.
4. K. W. Moore and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 012109.
5. K. W. M. Tibbetts, C. Brif, M. D. Grace, A. Donovan, D. L. Hocker, T.-S. Ho, R. B. Wu and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 86, 062309.
6. T. Blasi, M. F. Borunda, E. Räsänen and E. J. Heller, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 241303.
7. M. Nest, M. Ludwig, I. Ulusoy, T. Klamroth and P. Saalfrank, J. Chem. Phys., 2013, 138, 164108.
8. P. Doria, T. Calarco and S. Montangero, Phys. Rev. Lett., 2011, 106, 190501.
9. A. Rahmani, T. Kitagawa, E. Demler and C. Chamon, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 87, 043607.
10. E. Räsänen and E. J. Heller, Eur. Phys. J. B, 2013, 86, 1.
11. C. Brif, M. D. Grace, M. Sarovar and K. C. Young, New J. Phys., 2014, 16, 065013.
12. R. S. Judson and H. Rabitz, Phys. Rev. Lett., 1992, 68, 1500.
13. A. Weiner, Rev. Sci. Instrum., 2000, 71, 1929.
14. R. Bücker, T. Berrada, S. van Frank, J. Schaff, T. Schumm, J. Schmiedmayer, G. Jäger, J. Grond and U. Hohenester, J. Phys. B: At., Mol. Opt. Phys., 2013, 46, 104012.
15. N. X. Truong, P. Hilse, S. Göde, A. Przystawik, T. Döppner, Th. Fennel, Th. Bornath, J. Tiggesbäumer, M. Schlanges, G. Gerber and K. H. Meiwes-Broer, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 81, 013201.
16. R. Hildner, D. Brinks, J. B. Nieder, R. J. Cogdell and N. F. van Hulst, Science, 2013, 340, 1448.
17. R. Chakrabarti and H. Rabitz, Int. Rev. Phys. Chem., 2007, 26, 671.
18. H. Rabitz, M. Hsieh and C. Rosenthal, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 72, 052337.
19. M. Hsieh and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 77, 042306.
20. K. W. Moore, R. Chakrabarti, G. Riviello and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 83, 012326.
21. T.-S. Ho and H. Rabitz, J. Photochem. Photobiol., A, 2006, 180, 226.
22. A. Nanduri, A. Donovan, T.-S. Ho and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 033425.
23. R. Wu, H. Rabitz and M. Hsieh, J. Phys. A: Math. Theor., 2008, 41, 015006.
24. R.-B. Wu, R. Long, J. Dominy, T.-S. Ho and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 86, 013405.
25. M. Hsieh, R. Wu and H. Rabitz, J. Chem. Phys., 2009, 130, 104109.
26. J. Roslund and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 80, 013408.
27. K. W. Moore, M. Hsieh and H. Rabitz, J. Chem. Phys., 2008, 128, 154117.
28. A. Donovan, V. Beltrani and H. Rabitz, Phys. Chem. Chem. Phys., 2011, 13, 7348.
29. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2004.
30. A. Rothman, T.-S. Ho and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 72, 023416.
31. A. Rothman, T.-S. Ho and H. Rabitz, J. Chem. Phys., 2005, 123, 134104.
32. R. Magnanini, Appl. Anal., 1991, 41, 171.
33. D. Drucker and S. A. Williams, Am. Math. Mon., 2003, 110, 869.
34. J. Kennedy and R. Eberhart, Proceedings of the IEEE International Conference on Neural Networks, Perth, WA, IEEE, Piscataway, NJ, 1995, vol. 4, p. 1942.
35. J. Kennedy and R. Eberhart, Swarm Intelligence, Morgan Kaufman, San Francisco, CA, 2001.
36. N. Hansen and A. Ostermeier, Evol. Comput., 2001, 9, 159.
37. R. Ros and N. Hansen, Lect. Notes Comput. Sci. Eng., 2008, 296, 5199.
38. H. Rabitz, M. Hsieh and C. Rosenthal, Science, 2004, 303, 1998.
39. H. Rabitz, T.-S. Ho, M. Hsieh, R. Kosut and M. Demiralp, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 74, 012721.
40. Q. Sun, R. Wu, G. Riviello and H. Rabitz, to be published.
41. G. Riviello and H. Rabitz, to be published.

---