# Errors in quantum optimal control and strategy for the search of easily implementable control pulses

###### Abstract

We introduce a new approach to assess the error of control problems we aim to optimize. The method offers a strategy to define new control pulses that are not necessarily optimal but still able to yield an error not larger than some fixed a priori threshold, and therefore provide control pulses that might be more amenable for an experimental implementation. The formalism is applied to an exactly solvable model and to the Landau-Zener model, whose optimal control problem is solvable only numerically. The presented method is of importance for applications where a high degree of controllability of the dynamics of closed quantum systems is required.

## 1 Introduction

The theory of optimal control (OC) that has been mathematically formulated in the last century by the seminal works of Pontryagin, Bellman, Kalman, Stratonovich [1, 2, 3] has been instrumental for the achievement of highly reliable electronic devices used to control, for instance, mechanical systems such as airplanes, cars, etc., but also to control chemical reactions or to design ultra-fast laser pulses for manipulating molecules (e.g., to break a certain bond while leaving other bonds intact [4]), and today even to optimize (stochastic) financial analyses [5].

The topic has recently attracted the attention of physicists working in quantum information and computation science, because of the need to engineer accurate protocols. To this aim different numerical techniques have been devised in order to minimize (or maximize) some performance criterion or, alternatively called, objective functional. We mention the most used ones in open-loop quantum control: the Krotov iterative method [1, 6, 7] and the gradient ascent pulse engineering algorithm [8]. Although these are powerful tools for the search of OC pulses they do not provide an assessment of the tolerable error against distortions and do not guarantee to obtain control pulses easily realizable in the laboratory. Such an issue is of paramount importance for quantum information processing (QIP), since the error allowed by fault-tolerant quantum computation ranges between 0.01% to fractions of a percent [9, 10].

A possible (empirical) approach relies on applying an arbitrary distortion to the OC solution and then looking at the error that it produces on the objective functional, or by selecting a region in the Hilbert space that is robust against noise (decoherence free-subspace) [11], whose existence follows from the symmetry properties of the noise. Recently, a more systematic methodology based on the Hessian analysis of the cost functional has been proposed [12] or by using an improved genetic algorithm in the presence of control noise [13].

The aim of this work is to provide an alternative method to evaluate to which extent a given OC scheme can tolerate errors. The method is based on the Hessian approximation of the cost functional, as in Ref. [12], but it is applicable to any system Hamiltonian and Hilbert space dimension. Such arbitrariness is important in several circumstances where either the control pulse does not appear linearly in or where the state to be controlled is an auxiliary state (e.g., the motional state of an atom [14]) and not the quantum bit itself (e.g., an atomic internal state). The control of such states is relevant not only for several QIP implementations, but also for quantum metrological purposes, where the control of large quantum superpositions may increase the sensitivity of precise measurements.

Even though usually it is not possible to know analytically the OC pulse, we underscore that our assumption is that the parameter obtained with some numerical algorithm is very close to the global optimum. More precisely, the error on the cost functional obtained with the numerically found OC pulse has to be much smaller than the error allowed by the process we are interested to optimize. Besides the interest on its own, we believe that our approach might be of importance for experiments, where, typically, optimal pulses are extremely difficult to achieve. To this aim, our method could help to find easily implementable control signals (EICS), while still being able to satisfactorily fulfill the performance criterion we are interested in. Here with “easily implementable control signals” reference is made to pulses that can be utilized to control an experiment at the quantum level. More precisely, since nowadays the experiments are typically controlled by a computer, an obvious requirement for the control signal is that its Fourier spectrum has to match the bandwidth of the transducer or it can not vary faster than the clock frequency of the processor. Besides this, since the computer during the course of the experiment controls some device (e.g., the applied voltage on electrodes or electric current [14, 15]) the control pulse has, for instance, to take into account the bandwidth of those devices. These conditions might be not satisfied by the optimal control pulse obtained with the aforementioned optimization algorithms. Even though, recently, some extensions of those optimization methods in order to include spectral constraints on the control pulses have been made [16, 17], these are not always easy to be handled, especially when the dynamics of a many-body quantum system is concerned.

## 2 The method

Let us briefly review what is the purpose and what are the methodologies adopted commonly in the quantum control research area concerning closed quantum systems and the link between OC theory and analytical mechanics. The usual problem is the engineering of a system Hamiltonian such that the initial state at time of the quantum system under consideration is brought at time to some desired goal state (see also Fig. 1). To this aim there are two (equivalent) techniques for searching optimal pulses at our disposal: the variational method and dynamical programming. In both cases the goal is the minimization of the cost functional

(1) |

over all admissible control pulses and state trajectories
given a certain initial
state .
Here by admissible we mean any for which the Schrödinger equation is
well-defined and has a unique solution given the initial
condition . The first term on the r.h.s. of
Eq. (1) is the terminal functional, e.g., the overlap infidelity . The second term provides
additional constraints on the control pulse ^{1}^{1}1For instance,
the “laser electric field fluence” with , where is an electric field
amplitude.. In the variational
method the additional constraint
is introduced [18], where we set and is a Lagrange multiplier often
referred to as costate, which ensures that the state satisfies the
Schrödinger equation. The search of an extremum of the cost
functional produces a set of equations for the state, costate and the
control pulse.

On the other hand, dynamical programming, based on the Bellman’s optimality principle [19], produces an equation, the so called Hamilton-Jacobi-Bellman equation, for the optimal cost . This equation is expressed in formally the same way as the Hamilton-Jacobi equation of classical mechanics, but with the difference that it is propagated backwards in time. The role of the Hamilton function in classical mechanics is played in control theory by the (quantum) Pontryagin Hamiltonian [20, 21], which is not to be confused with the system Hamiltonian that we control through . From the Hamilton-Jacobi-Bellman equation it is possible to retrieve the equations of motion for the state and the costate, the same as in the variational approach, which look, given the aforementioned Pontryagin Hamiltonian, formally as the Hamilton equations for the phase space variables in analytical mechanics. Beside this, the solution of the Hamilton-Jacobi-Bellman equation, defined through the solution of the latter (the Hamilton boundary-value problem), is given by [21]:

(2) |

which is very similar to the action in classical mechanics. We note, however, that the variational approach, based on the so-called Pontryagin maximum principle [19], yields necessary and sufficient conditions for local minima, whereas dynamical programming produces results that are globally optimal.

Motivated by this analogy between OC theory and analytical mechanics we can view the cost functional as an “action functional”. Indeed, in analogy to the Hamilton’s principle where the actual evolution of a classical system is an extremum of the action functional, which produces the well-known Euler-Lagrange equations of motion, an extremum of produces the OC equations for the state and the costate.

Now, let us describe how our method works. To begin with, we fix the desired cost, , that the system (at least) has to attain. Since the advantage of optimizing quantum dynamics is connected with the possibility of reaching by exploiting the interference of several paths in the space of control parameters , we define a path integral in such a space. To this aim, we introduce the weight with being some suitable measure on the space . The form of this weight resembles the Feynman propagator, which is motivated by the previous discussed analogy between analytical mechanics and OC theory, and by the expression (2) for the optimal cost. The choice of the exponential, however, does not emerge from fundamental physical requirements. We think that any well-behaved function peaked around the optimal control pulse will allow to estimate the robustness of an optimal control problem. This conjecture is based on the observation, see the discussion in the following, that according to our analysis the curvature of the cost functional around the optimal solution is the relevant quantity to analyze.

The numerical factor given in the weight , which we shall refer to “infidelity tolerance”, has the following property: when approaches (from above) the minimum of , then . Hence, is the analogous of in Feynman path integral. Indeed, when we retrieve the OC equations for the state and the costate we mentioned before. Besides this, we assume that exists and it is unique for a given OC problem, regardless of the form of the distortion . The uniqueness is first of all a necessary condition or else our method would be ineffective. There is, however, a more deep reason why we make such an assumption. This is more easily understood in the case where in Eq. (1). In such a scenario the cost functional relies solely upon the state at the final time , that is, the cost functional will depend on some time integral of over (and eventually on other time integrals for its time derivatives over depending on the particular control problem one is interested in). Thus, minimizing the cost functional means to find the right conditions for those integral functionals which depend only on time . Hence, the form of in the interval does not matter, if those functionals of the control pulse and its time derivatives fulfill the right conditions at the final time . In this respect the infidelity tolerance has to be independent from the distortion we utilize, that is, it is unique. This is also what is shown by the analysis carried out on the examples we will discuss later in the paper.

Even though it relies on the particular control problem we have at hand, from the above outlined discussion, at least in the most relevant cases for QIP where in Eq. (1), the set of EICS, , is dense, since what matters is the fulfillment of the right conditions at the final time . For instance, in the first example we are going to consider in the next section, that is, the optimal transport of a particle confined in a moving harmonic trap, Ref. [22] has showed that if is the optimal solution then is optimal as well (here , see also Sec. 3). Thus, there is a group of optimal solutions, parametrized by the continuous variable , which is topologically dense. Similarly, this will occur for , whose set of control pulses forms another (dense) subset in . Each of these possible control signals will have a precise spectrum that has to be within the bandwidth of , namely the largest bandwidth of the elements of . Now, if the EICS needs to have a specific bandwidth, then one has to select from the that have the right bandwidth, namely restrict (continuously) the bandwidth of such that the right subset of becomes , even though such a procedure might produce an empty set.

To be specific let us set, for the sake of simplicity, and consider only one control pulse , that is, only one control pulse is applied to the system we aim to steer during the time interval . We also assume that the state obeys the Schrödinger equation with the time-dependent Hamiltonian operator and initial condition . Contrarily to (1), where the state is an independent variable, hereafter we render explicit the dependence of on the control parameter . Thus, we introduce the reduced cost functional . Then by performing the Taylor expansion of around to second order in we obtain

(3) |

where is some norm in , is the minimum, which, without loss of generality, we will set to zero, is the transposed of the vector , and the time interval is divided in equal parts . The Hessian is a real, symmetric, and positive defined matrix, and therefore it can be diagonalized. Additionally, we assume the boundary conditions , that is, the values and are fixed. We underscore that henceforth we shall work with the discretized system time evolution since in most cases the optimization of a given control problem is performed numerically, and therefore it is discretized, but most importantly, because any experiment is performed with a finite number of control time steps.

Given that, let us define a suitable norm in for by means of , in order to obtain a quantitative appraisal of the tolerated error by a given OC scheme. A suitable choice for such a norm is given by: with . Because of the second order approximation in Eq. (3), the integrals appearing in over the control pulse at different times can be replaced (by making a change of integration variables) with . This norm is indeed an “average” in of the reduced cost functional itself with the exponential function being a “probability distribution”. A similar approach has been introduced by Rabitz [23] where the terminal functional in Eq. (1) is averaged over a distribution function . Within our method can be identified with the exponential function in . We remark that we are not interested in some specific noise model, but rather to identify the portion which enable us to satisfy . Basically, is determined by , which only relies on and , and any kind of noise has to yield pulses such that . Hence, our approach is applicable to any OC problem and cost functional (1). Furthermore, we note that in the next we shall consider only real-valued (e.g., an electric current [14]), but we underscore that the path integral can be easily generalized to complex-valued like the amplitude and the phase of a laser field.

Close to the optimal solution we can approximate with its second order expansion, and therefore the norm becomes

(4) |

where are normalization factors, and are the non-zero eigenvalues of with . The above formula makes good sense, because when also .

In order to apply the above outlined formalism to some concrete example we shall consider hereafter , with . We note that the state is implicitly depending on the whole history of . Assuming that is a differentiable functional of its arguments we have

(5) |

Here we used the fact that . The most difficult part of our method is the computation of the Hessian matrix . To this aim we have two possibilities at our disposal: either we estimate by means of the Broyden-Fletcher- Goldfarb-Shanno (BFGS) formula [24] or we compute the first and second derivatives of the state with respect to the control . Here we choose the latter approach, where we have to solve the following equations

(6) |

which apply to any quantum closed system. Here is the solution of the Schrödinger equation for , and , are the Gateaux derivatives of the system Hamiltonian (defined as: ). These derivatives are always analytically computable, since the dependence of the system Hamiltonian on the control pulse is always known, whereas the analytical dependence of both and on is known only in few cases (e.g., the driven and parametric harmonic oscillator [25]). Because of the latter, we need to solve the equations (6) which allow us to determine the matrix . However, depending on the particular control problem, the BFGS method might be more efficient. Beside this, we note that the case of a linear quadratic regulator (or even Gaussian) control, that is, a system for which the state equation is linear and the performance criterion to be minimized is a quadratic form of the state and eventually also of the control pulse [19], is not contemplated in our scheme. Indeed, the system state is always dependent on the control pulse through the Schrödinger equation of motion, even in the simple scenario where the system Hamiltonian depends linearly on . Indeed, we are interested in the Hessian of the reduced cost functional , whose dependence on might be not trivial through . Hence, the cases in which the Hessian relies only upon the state do not concern .

We also remark that the above outlined formalism concerns state vectors in Hilbert spaces. If we would be interested in the optimization of unitary transformations , then we should perform the Taylor expansion of , where is the ideal unitary we wish to accomplish, and is the dimension of the Hilbert space. Thus, the previuos analysis can be easily generalized to unitary operations.

## 3 Applications of the method

Let us first apply our method to an exactly solvable model. We consider the transport of a particle in a movable one-dimensional harmonic trap potential, for which the OC pulse and the functional dependence of the state on the controller are analytically known [22]. Such control problem is of relevance for the realization of a quantum ion processor. Indeed, optimistic estimates show that transport processes may account for 95% of the operation time of a quantum computation [26]. Beside this, the harmonicity of the ion confinement in segmented Paul traps has been proven both numerically [27] and experimentally [28].

The system Hamiltonian is given by with (we use harmonic oscillator units). We focus our attention on the ground state, that is, when the particle is initially prepared in the lowest vibrational state of the trap. The aim is to transport such a state over the distance in a time such that , where is an unimportant phase factor and is the Gaussian harmonic oscillator ground state wavefunction. Thus, our goal is to find a prescription such that when the OC pulse is perturbed the system has to reach at least the - a priori fixed - value of fidelity .

In Ref. [22] the analytic solution for the time evolved state is provided. This enables us to compute analytically the Gateaux derivatives of the state of the system without the need of solving (6). In Fig. 2(a) we show results for a distortion given by: , with being the OC pulse of Ref. [22] [see Eq. (5) therein]. Such simple distortion modulates the OC pulse at the rate and gives a direct quantitative measure of the distortion degree applied to : the larger is, the larger the infidelity. Thus, Fig. 2(a) shows which is the largest admissible value of for a fixed value of the (reduced) cost functional, whereas in the inset we show the deviation from the linear behaviour of for .

Now we write the cost functional as . In this specific example it turns out, numerically, that the matrix elements of are of the form , with such that . This means that the matrix elements are almost equal to each other. The eigenvalue problem for such a matrix can be well approximated by , that is, only an eigenvalue is non-zero. Thus, by defining we get for the non-vanishing eigenvalue the simple expression . Given these remarks, the norm (4) reduces to .

Finally, in order to determine we proceed in the following way: since the norm is the “average cost functional” and has to be unique, we simply use the inverted formula for a given choice of , and perform the substitution . (For several non-zero eigenvalues of we would have had .) Then, we choose some distortion and by varying the strength of such a distortion we collect the values of versus the numerically exact overlap infidelities (i.e., without second order approximation), which is basically identified with , the fixed error threshold. We tested, however, that different kinds of distortions (e.g., Fourier-like model) produce practically the same curve as the one shown in Fig. 2(b). This is easily understood, since what is relevant for the overlap infidelity terminal functional is the final state . Hence, our method is general and it applies to any kind of distortion model. The reason for choosing the single frequency noise for the results displayed in Fig. 2(a), was only to show the impact of the distortion on the cost functional in a simple and analytical way.

Secondly, we perform a fit of the obtained curve for as a function of the overlap infidelity. For the present example, showed in Fig. 2(b), we found that the following function well represents the data: , with and . Given that, all distortions of the optimal control pulse that satisfy the inequality for a fixed (a priori) value of desired fidelity , will yield an overlap fidelity . From an operational point of view, this means that if we randomly choose a distortion such that , than we will certainly attain a state whose overlap fidelity is greater than (see also Fig. 1). This permits us to find control pulses that might be more appropriate for an experimental implementation. This conclusion follows from the procedure we established for the determination of the infidelity tolerance . We underscore, however, that the threshold depends only on , because relies only on , and that it is the result of an average in through the path integral we defined. This provides a “universal” character to for the given control problem.

As a second example we consider the Landau-Zener model which has been proven useful to describe the tunneling of Bose-Einstein condensates in accelerated optical lattices [29] and the dynamics of a quench-induced phase transition in the quantum Ising model [30]. The model is described by the following system Hamiltonian: , where are Pauli matrices. The goal is to bring the system from the ground state of the Hamiltonian to the ground state of the Hamiltonian through the avoided level crossing. As objective functional we consider , that is, we also control the phase of the state. Even though the analytical dependence of the system Hamiltonian is known, and therefore the Gateaux derivative (), for such a control problem both the and the dependence of on are not analytically known. Concerning the latter, this implies that is not possible to compute analytically the Gateaux derivatives and , and therefore the Hessian . Hence, we need to solve (6). Given that, we seek an by using the Krotov iterative method [1, 31]. We then proceed on, as in the former example, by determining the infidelity tolerance , for which we get a similar fit with , , and . The criterion to be satisfied is then again given by: , but with a different numerical value for .

In Fig. 3 it is showed the exact result of the infidelity (red) and the second order approximation (black) for 50 random realizations of for a distortion similar to the one of Fig. 2(a). Instead, the desired upper limit of infidelity , that the system has at least to attain, is represented by the horizontal blue line. Analogously to the former example, only the realizations of that fulfil have an infidelity below 1%. As for the previous example, the single frequency model has been adopted for convenience in order to illustrate the effectiveness of our approach as depicted in Fig. 3. Such a choice, however, does not invalidate our methodology, which we have tested for a Fourier-like distortion model (i.e. a sum over a finite number of harmonics with different amplitudes), similarly to the above outlined determination of the infidelity tolerance.

In general, in order to identify robust control pulses amenable to an experimental implementation one should proceed as follows: 1) determine the optimal control with some numerical optimization algorithm; 2) compute the eigenvalues of the Hessian matrix ; 3) determine the infidelity tolerance as described in the first example; 4) randomly generate distortions such that the condition is fulfilled for a fixed (a priori) value of ; 5) define the new control pulse as ; 6) choose the most suitable control pulse to be employed in the laboratory from the ensemble . We note, however, that such a “recipe” does not guarantee that we are able to obtain with certainty a control pulse more amenable for use in an experiment, but at least it helps to design new ones that are still able to yield a value of the cost functional below the fixed threshold . In other words, the set might be empty, and the elements of are not necessarily close to the optimal (e.g., see Ref. [22]). In this respect, such a procedure might help the quest of both robust and experimentally feasible pulses for controlling different quantum phenomena.

## 4 Conclusions

In conclusion we have presented a method to assess the error of solutions to OC problems. The method might be a helpful tool for experimentalists in order to design experiments robust against source of noise. For instance, to estimate how much the schemes for realizing quantum gates are robust against imperfections of the optimal pulse shape. Compared to other methods, such as the one of Ref. [13], our technique does not need the simulation of a large ensemble of samples in order to minimize on average both the mean and the variance of the objective functional. Instead, once the Hessian matrix and the infidelity tolerance are known, which only rely on the optimal (known) control pulses, one has simply to randomly generate the distortion and perform the matrix multiplication , which is a less demanding computational task than the application of a genetic algorithm, as the one proposed in Ref. [13]. Besides this, our method can be applied to both known error models and to evaluate unknown errors such as random telegraphic noise. For the future, we plan to extend our formalism to dissipative quantum systems (e.g., systems governed by the Born-Markov master equation).

## Acknowledgments

We acknowledge financial support by the IP-SOLID and the National Research Foundation and Ministry of Education Singapore (R.F.), IP-AQUTE and PICC (T.C.), SFB/TRR21 (A.N.,T.C.), the Marie Curie program of the European Commission (Proposal No. 236073, OPTIQUOS) within the 7th European Community Framework Programme and the Forschungsbonus of the University of Ulm and of the UUG (A.N.). A. N. acknowledges conversations with K. Urban, J. T. Stockburger, and I. Kuprov.

## References

## References

- [1] V. F. Krotov, Global methods in optimal control theory, Monographs and Textbooks in Pure and Applied Mathematics 195 (Marcel Dekker, New York, 1996)
- [2] A. C. Doherty, S. Habib, K. Jacobs, H. Mabuchi, and S. M. Tan, Phys. Rev. A 62, 012105 (2000)
- [3] P. S. Maybeck, Stochastic Models, Estimation, and Control, Mathematics in Science and Engineering, Volume 141-1 (Navtech Book & Software Store, Arlington VA, 1994)
- [4] T. Brixner, G. Krampert, T. Pfeifer, R. Selle, G. Gerber, M. Wollenhaupt, O. Graefe, C. Horn, D. Liese, and T. Baumert, Phys. Rev. Lett. 92, 208301 (2004)
- [5] B. Øksendal, Stochastic Differential Equations (Springer-Verlag, Heidelberg, 2000)
- [6] S. E. Sklarz and D. J. Tannor, Phys. Rev. A 66, 053619 (2002)
- [7] D. Reich, M. Ndong, and C. P. Koch, arXiv:1008.5126v1
- [8] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, J. Magn. Reson. 172, 296 (2005)
- [9] A. M. Steane, Phys. Rev. A 68, 042322 (2003)
- [10] E. Knill Phys. Rev. A 71, 042322 (2007)
- [11] P. Zanardi and M. Rasetti, Phys. Rev. Lett. 79, 3306 (1997)
- [12] T.-S. Ho, J. Dominy, and H. Rabitz, Phys. Rev. A 79, 013422 (2009)
- [13] F. Shuang, and H. Rabitz, J. Chem. Phys. 121, 9270 (2004)
- [14] P. Treutlein, T. W. Hänsch, J. Reichel, A. Negretti, M. A. Cirone, and T. Calarco, Phys. Rev. A 74, 022312 (2006)
- [15] P. Böhi, M. F. Riedel, J. Hoffrogge,J. Reichel, T. W. Hänsch, and P. Treutlein, Nat. Phys. 5, 592 (2009)
- [16] C. Gollub, M. Kowalewski, and R. de Vivie-Riedle, Phys. Rev. Lett. 101, 073002 (2008)
- [17] M. Lapert, R. Tehini, G. Turinici, and D. Sugny, Phys. Rev. A 79, 063411 (2009)
- [18] A. P. Peirce, M. A. Dahleh, and H. Rabitz, Phys. Rev. A 37, 4950 (1988)
- [19] D. E. Kirk, Optimal control theory (Dover Publications, Inc. Mineola, New York, 2004)
- [20] V. P. Belavkin, A. Negretti, and K. Mølmer, Phys. Rev. A 79, 022123 (2009)
- [21] G. Gough, V. P. Belavkin, and O. G. Smolyanov, J. Opt. B: Quantum Semiclass. Opt. 7, S237 (2005)
- [22] M. Murphy, L. Jiang, N. Khaneja, and T. Calarco, Phys. Rev. A 79, 020301 (2009)
- [23] H. Rabitz, Phys. Rev. A 66, 063405 (2002)
- [24] A. Mordecai, Nonlinear Programming: Analysis and Methods (Dover Publishing, New York, 2003)
- [25] D. C. Khandekar, and S. V. Lawande, J. Math. Phys. 20, 1870 (1979)
- [26] G. Huber, T. Deuschle, W. Schnitzler, R. Reichle, K. Singer, and F. Schmidt-Kaler, New J. Phys. 10, 013004 (2008)
- [27] K. Singer, U. Poschinger, M. Murphy, P. Ivanov, F. Ziesel, T. Calarco, and F. Schmidt-Kaler, Rev. Mod. Phys. 82, 2609 (2010)
- [28] G. Huber, F. Ziesel, U. Poschinger, K. Singer, and F. Schmidt-Kaler, Appl. Phys. B: Lasers Opt. 100, 725 (2010)
- [29] A. Zenesini, H. Lignier, G. Tayebirad, J. Radogostowicz, D. Ciampini, R. Mannella, S. Wimberger, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 103, 090403 (2009)
- [30] W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005)
- [31] S. Montangero, T. Calarco, and R. Fazio, Phys. Rev. Lett. 99, 170501 (2007)