A general framework for quantum splines

Quantum splines are curves in a Hilbert space or, equivalently, in the corresponding Hilbert projective space, which generalize the notion of Riemannian cubic splines to the quantum domain. In this paper, we present a generalization of this concept to general density matrices with a Hamiltonian approach and using a geometrical formulation of quantum mechanics. Our main goal is to formulate an optimal control problem for a nonlinear system on the dual of the Lie algebra of the unitary group which corresponds to the variational problem of quantum splines. The corresponding Hamiltonian equations and interpolation conditions are derived. The results are illustrated with some examples and the corresponding quantum splines are computed with the implementation of a suitable iterative algorithm


Introduction
Quantum control theory, and in particular quantum optimal control, is drawing considerable attention from quite different communities in different areas of physics, chemistry, applied mathematics and quantum information (see [23] for a nice introduction to the area with several practical physical and chemical applications). This subject represents an essential ingredient to use the new quantum technologies which are being created in these fields, with a wide range of applications (see, for instance, [5,6,16,18,22,24,25,26]). An important problem is the design of efficient algorithms to generate accurately specified quantum states under certain conditions (in minimal time, with a minimal energy investment, etc). These algorithms usually consist in defining a suitable time dependent Hamiltonian which drives the system to the desired quantum states under the required conditions.
In a recent paper [9], Brody and co-workers introduced the notion of quantum spline as the generalization to the quantum domain of the familiar notion of geometric splines (see [12,13,15,21]). A quantum spline is the solution of the following quantum optimization problem. Let us consider a set of states {|ψ 0 , · · · , |ψ N } in a certain finite dimensional Hilbert space H and a corresponding set of times {t 0 , · · · , t N }. The problem is to find a time-dependent Hamiltonian H(t) defining an unitary evolution which, at the given times t k , passes arbitrarily close to the given states |ψ k , and such that the total change in the Hamiltonian along the curve is minimized. We can write it as a variational problem asking to minimize the cost functional where |ψ(t) = U (t)|ψ 0 represents the trajectory in H defined by the unitary evolution U (t), the scalar product · | · represents the trace of the product and is just a positive coefficient which allows us to change the relative importance of the total change of the Hamiltonian of the system and how close we are from the point we are supposed to target (D represents the corresponding distance).
Some comments are in order: • As the time dependence of the Hamiltonian in quantum control problems is often implemented as time-dependent magnetic fields which interact with the spin system, we can think at the minimization problem as the equilibrium between the accuracy at targetting the points at the designated times and the total energy required to change the magnetic fields during the process. • The (square) distances D 2 (|ψ(t k ) , |ψ k ) are considered on the projective space PH, with respect to the canonical Fubini-Study metric. From the technical point of view, the functional (1) combines thus quantities defined on two different manifolds: the trace-norm of skew-Hermitian operators and the Fubini-Study metric on the projective space. • The solution of the problem above, even for very simple examples such as a two-level quantum system, requires of numerical algorithms. In [9] the authors provide an algorithm for the discretization of their variational equations which is used in the construction of the solution of the example they study.
The aim of this paper is to introduce an alternative formulation of the problem of quantum splines and a generalization of the notion which enlarges the range of potential applications. The main idea is to re-formulate the problem in terms of density matrices and to adapt the conditions on the dynamics to it. Thus, the notion of quantum spline for pure states introduced in [9] is included, while an analogous notion for mixed states does also make sense: New formulation of the interpolation problem: Let us consider a set of density matrices {ρ 0 , · · · , ρ N } on a certain finite dimensional Hilbert space H and a corresponding set of times {t 0 , · · · , t N }. Find a timedependent Hamiltonian H(t) defining an unitary evolution which, at the given times t k , passes as close as possible to the given points ρ k , and such that the change in the Hamiltonian is optimal in the sense that the cost functional is minimized, where the trajectory of the system represented by ρ(t) is a solution of von Neumann equation: and d represents the distance on the space of Hermitian operators This formulation of the problem includes the previous one for the case of pure quantum states, but it allows to consider it also for general quantum states. Furthermore, it relaxes the continuity condition on the derivative of the Hamiltonian H(t) at the intermediate points, as it is often assumed in (quantum) control based on pulses. From a physical point of view, the lack of continuity does not represent a serious condition on the system, and as we are going to see it allows for a much simpler analysis of the problem.
Regarding the generalization to density matrices, notice that all real systems are in contact with some type of environment. In that case, we know that the system is very seldom in a pure quantum state and therefore to extend the notion of quantum spline to mixed states seems as a natural step forward. For instance, we could consider a system which initially is in equilibrium at a finite temperature T (see [7]) with a bath. If the evolution associated to the quantum controls is much faster than the dynamics produced by the interaction with the bath, we could approximate the evolution by a unitary transformation on the set of mixed quantum states. This is the type of problem that we will be considering in this paper. Another interesting generalization would be to consider that the controlled evolution is of the same order as the interaction with the bath and replace von Neumann equation by a more general master equation, such as Lindblad-Kossakowski equation. The problem would be formally analogous to our "New formulation" above, replacing equation (4) by Lindblad-Kossakowski equation. Nonetheless, Lindblad-Kossakowski equation would put some constraints in the possible values of the set of times {t j } and density matrices {ρ j } (remember that Lindblad dynamics is not always controllable [2,3]). We are also studying this problem and it will be considered in a future paper.
There are important aspects that characterize the new formulation: • From the technical point of view, our formulation is much simpler than the one considered in [9], since it can be done at the level of a linear space (the dual space to the algebra of the unitary group u * (n)), instead of considering the Lie group U (n) and the projective space CP n−1 . • The set of quantum states D for our formulation will be the submanifold of u * (n) defined by the elements which have trace equal to one and are positive definite. D is a nonlinear submanifold of u * (n), but the particular choice of our dynamic system allows us to treat the problem as a problem defined on the linear space u * (n). Indeed, as the dynamics on D is defined by von Neumann equation, we know that the solution ρ(t) must belong to a unitary orbit of the coadjoint action of the unitary group. As it is well known (see [20]), those orbits constitute symplectic submanifolds of u * (n) which define the leaves of the symplectic foliation associated with its Lie-Poisson canonical structure. If we consider an initial condition contained in one of the leaves, the whole solution will be contained in it. This property will allow us to consider the problem defined directly at the linear space u * (n) and forget about the nonlinear constraints. Nonetheless, this also imposes some constraints for the definition of splines in the general case. Indeed, as the dynamics we are considering is always unitary, any generalized spline will be contained in a unitary orbit. Thus, if the target points are not contained in one unitary orbit, the best possible solution will only be able to define the closest unitary orbit to the set of target points. We will discuss this point at the end of the paper. • The geometric formulation of Quantum Mechanics (see [4,8,10,19] ) is a natural framework to formulate the problem in a geometric formalism which exhibits these symplectic aspects. In that framework, von Neumann equation defines in a natural way a Hamiltonian vector field. This property is particularly useful when considering numeric integration of the solutions, since by using symplectic integrators we do not need to take into account constraints such as the trace of the quantum state or its rank and the problem can be treated as a free one. • In [9], a numerical algorithm was introduced in order to provide a method of integration for the splines. In our formalism, even if simpler, a numerical algorithm will also be necessary. We introduce an iterative algorithm to approach the optimal solution of the interpolation problem. We may not reach the absolute optimal solution, but we always obtain a very good approximation to it. As our framework is defined on linear spaces the algorithm turns out to be very simple and the efficiency is very high.
The structure of the paper is as follows. Section 2 presents a summary of the geometric formalism of quantum mechanics, with special emphasis on the geometrical structures associated with the description of the set of Hermitian operators and the set of density states, which are the main ingredients of our construction. Section 3 presents the main aspects of our new formulation, analyzing first the quantum control problem from the point of view of the Pontryagin maximum principle, then the analytical tools used to build the system of differential equations defining the solution and finally discussing an iterative algorithm to build solutions in an efficient way and the numerical methods we use to obtain them. Section 4 illustrates our construction with the simplest example for the interpolation between a set of pure states of a two level system and Section 5 presents a more sophisticated example of the interpolation of mixed states of a three level quantum systems. Finally, Section 6 summarizes the main results of the paper and discuss the possible generalization of our method to the case of open quantum systems.

The geometrical description of quantum mechanics
Let us briefly review some very general aspects of the geometrical formulation of quantum mechanics, focusing on the ingredients which will be used later in the paper. For more details, we address the interested reader to [4,8,10,19] and references therein. We use Einstein summation over repeated indexes.
2.1. The geometrical structures of u * (n). We hereafter assume H to be an n-dimensional complex Hilbert space. In this case the C * -algebra of operators corresponds to A = End(C n ) and the involution is the usual adjoint operation for complex endomorphisms. In this case, when considering an orthonormal basis, self-adjoint operators are represented by Hermitian matrices.
It is important now to emphasize the following facts: • We can define isomorphisms between u(n), u * (n) and the vector space of Hermitian matrices iu(n) by mappings These isomorphisms define also Lie brackets on each space [·, ·] u , [·, ·] u * and [·, ·] iu . We have a special interest in the identification of the elements of u * (n) with the Hermitian matrices, which will be used later. • On the dual space u * (n) we have a natural linear Poisson structure, the Lie-Poisson structure. The Poisson tensor Λ is defined, on the set of linear functions on u * (n), as follows From this definition on linear functions, we can construct the Poisson bracket on general functions f, g ∈ C ∞ (u * (n)) by extending equation (7) requiring bilinearity, i.e., as the action of the bi-differential operator which on linear functions behaves as equation (7). This defines the canonical Lie-Poisson tensor of the dual of the unitary algebra (see [1]). • The Hamiltonian vector field X iH = −Λ(iH, ·) is the infinitesimal generator of the coadjoint action of U (n) on u * (n) and, therefore, the systeṁ is a symplectic Hamiltonian system on the coadjoint orbit where the initial condition ξ 0 is contained (again we refer the reader to [1] for details) .

The set of density matrices D(H).
The physical states which are used in the Schrödinger formalism correspond to the points of the projective space PH. As it is well known, we can also represent these states by using the projectors on onedimensional subspaces of the Hilbert space. Unitary evolution, associated with the Schrödinger or the von Neumann equation, will define trajectories on this set. We will denote by D 1 (H) the corresponding set of projectors. Nonetheless, D 1 (H) is not enough to represent all the possible physical states of a system, since arbitrary convex combinations of rank-one-projectors also define admissible physical states. Indeed, if our system is not isolated and is surrounded by an environment (as it is the case of all real systems), or just a set of other systems, the representing state will not be, in general, pure. Therefore, we must enlarge the set of states to consider: Definition 2.1. The set of density states D(H) of the system corresponds to the subset of u * (n) obtained by convex combinations of rank-one-projectors, that is, Equivalently, we can consider the following definition: an element ρ ∈ u * (n) is a density operator if and only if An important property of the manifold D(H) is its internal structure, in particular in what regards the restriction of the coadjoint orbits of the unitary group on u * (n). A very important aspect of that structure is the nature of D(H) as a stratified manifold, the strata being defined by the rank of the states (see [?] for details). Indeed, as all density states are self-adjoint and positive definite, they all are diagonalizable and have non-negative eigenvalues. The number of non-vanishing eigenvalues equals the rank of the operator (or the matrix, if we choose a basis). Thus, all states having the same rank belong to the same stratum. As the coadjoint action of the unitary group is known to preserve the spectrum of the density operator, it is obvious that the orbit will stay on the same stratum as its initial condition. Equivalently, we can claim that the Hamiltonian vector field defined by equation (8) will be tangent to the different strata.
Thus, if we consider the Hamiltonian vector field X iH (defined by equation (8)) restricted to D(H), we obtain the von Neumann equation for the density matrix We know that this evolution corresponds to a unitary transformation. Thus, the Hamiltonian vector fields of Λ are the infinitesimal generators of the coadjoint action of the unitary group on u * (n), and the corresponding orbits are known to be symplectic submanifolds. For instance, the set of pure states D 1 (H) defines a symplectic orbit which corresponds to a symplectic leaf of the canonical Poisson foliation. Therefore, the solution of von Neumann's equation on D 1 (H) is always tangent to it and defines a symplectic transformation. This allows us to use the formulation of the problem on u * (n) without using any constraints to specify the submanifold of pure states, since we know that any trajectory with an initial condition on D 1 (H) will remain on the submanifold for all times.
For the case of a mixed state, the situation is analogous, but in this case the symplectic orbit of the coadjoint action will not cover all the stratum. Each stratum will contain several coadjoint orbits of different dimensions and therefore the formulation of the interpolation problem must take into account this property. Indeed, the exact quantum spline for mixed states (i.e., a curve joining all points) makes sense only if the set of points belong to the same coadjoint orbit. As the target points are not reached exactly, the problem still makes sense for a neighborhood of the orbit. In the next Section we will see how it is possible to define an iterative algorithm which drives the evolution to the closest point of the orbit with respect to the target point.
3. Dynamical interpolation problem on u * (n) Consider also a m-dimensional vector subspace h of u * (n), m ≤ n 2 which will represent the space of admissible controls for our interpolation problem.
3.1. The control problems chain. In general, let E be a real vector space of dimension n and U a subspace of E (in our case, E will be the vector space u * (n) × u * (n)). Given a cost function L : E × U → R and the functions G j : E → R, consider the optimal control problem P concerning the minimum of the functional where is a tunable positive parameter, and with u : [0, T ] → U and x : [0, T ] → E being curves, subject to the initial condition x(0) = x 0 , satisfying the control systeṁ and the conditions: -x is continuous in [0, T ] and smooth in ]t j−1 , t j [, with finite one-sided limits at t j , for all j = 1, . . . , N ; -u is smooth in ]t j−1 , t j [, with finite one-sided limits at t j , for all j = 1, . . . , N . This defines a set of N subproblems P j , j = 1, . . . , N , one for each subinterval, where the initial condition is fixed (by the value of the solution at the end point of the previous interval), but the final condition is free. Thus, on interval ]t j−1 , t j [ we consider the functional and search for the optimal solution x : [t j−1 , t j ] → R subject to the initial conditions x(t j−1 ) = x j−1 , which corresponds to the final conditions of the previous subinterval. Each of these subproblems is a particular case of a well known optimal control problem known as a Bolza problem (see [11]). The solution x(t) is asked to be continuous on the whole interval [0, T ]. The solution of each Bolza problem can be obtained straightforwardly by using the Pontryagin maximum principle, the only difficulty being the boundary condition. Although it is a well known method, let us present, for completeness, a sketch of the variational construction needed to prove Pontryagin's result. Let us take, for instance, the first subinterval [t 0 , t 1 ]. Consider the Hamiltonian function h : E × E * × U → R associated with the optimal control problem, defined by For each x ∈ E and p ∈ E * , the partial functional derivative of a function h with respect to x at a is an element of E * denoted by (δh/δx)(a) and the partial functional derivative with respect to p at a is an element of E denoted by (δh/δp)(a). Now, let us enlarge the functional J by the dynamical constraints, using the costate trajectory p : [t 0 , t 1 ] → E, in the following way: Thus, we know that the optimal solution will satisfy the equationṡ Therefore, it is now simple to prove the following result: Theorem 1. If u is the optimal control resulting of the optimal control problems P j , j = 1, . . . , N and x is the associated optimal state trajectory, then there exists a piecewise-smooth costate trajectory p : [0, T ] :→ E such thaṫ

Proof
The problem is defined piecewisely, and hence we will proceed in the same way. For each j, consider the problem on the subinterval [t j−1 , t j ], determine the optimal solution on it, and consider then the value of the optimal solution to define the boundary condition of the next subinterval in order to ensure continuity of the trajectory. On the first subinterval, [t 0 , t 1 ], we have seen that the solution satisfies the condition. Consider now the optimal solution x(t) and the corresponding value x(t 1 ), and use it as the initial condition for the Bolza problem in the subinterval [t 1 , t 2 ]. We will then define an optimal solution for the second subinterval which will define the initial condition for the third one. On each case we produce a smooth costate variable p, defined on each subinterval. If we repeat the procedure on all subintervals, the result follows and the global costate variable will be, by construction, piecewise-smooth.
Observe that system (17)(18) and the time-invariance of the Hamiltonian h guarantee that h is an integral of motion.
3.2. The optimal control problem P for quantum splines. Let us apply now the result above to our formulation of the interpolation problem. As we had formulated it on u * (n), we consider a cost function L : where ξ, η = 1 2 Tr(ξη) for ξ, η ∈ u * (n). The functions G j : u * (n) → R are given by The dynamical interpolation problem is the optimal control problem concerning the minimum of the functional satisfying the control systemρ where the curves verify the same regularity conditions of the previous section. The Hamiltonian function to apply the Pontryagin maximum principle becomes now the function h : u * (n) × u * (n) × u(n) × u(n) × h → R given by The approach considered above, obtained by application of the theorem 1, implies the following result: Theorem 2. If u is the optimal control resulting of the optimal control problems P j , j = 1, . . . , N , and (ρ, H) is the associated optimal state trajectory, then there exists a piecewise-continuous optimal costate trajectory (Γ, Π) that satisfies the systeṁ the initial conditions (23) and the interpolating conditions In the equations above ad and ad * refer to the adjoint and co-adjoint actions of the Lie bracket [·, ·] u * and · specifies the isomorphism : u * → u given by the inner product on u * (n).
Notice that eliminating the controls in (26) we obtain the systeṁ which is Hamiltonian with respect to the functionh : Notice that · specifies the inverse of the isomorphism .
Observe that the Hamiltonian system (28) can be written as followṡ with interpolation conditions Consider now the problem P with u taking values in all space, that is, the space of admissible controls corresponds to h = u * (n) and thus m = n 2 . From Hamiltonian equations (30), we differentiate twice the equationḢ = Π and eliminate the costates Γ and Π, to obtain the equationsḦ = −ad ρ Γ and ...
Moreover, interpolation conditions can be interpreted in terms of ρ and H as follows: ρ is of class C 1 , H is of class C 0 and where K j are the matrices defined by Notice that any solution of equations (35-36) satisfies the second equation of (33), i.e., any of them are cubic polynomials and represent extremals of the continuous part of functional J defined in equation (22).
Let {σ l | l = 1, . . . , n 2 } be an orthonormal basis of u * (n). Denote, with respect to this basis, the coordinates of K j by (k l j ), l = 1, . . . , n 2 , and the coordinate functions for ρ and H by (x l ) and (y l ), l = 1, . . . , n 2 , respectively. The system (33) is then written in the following waẏ for l = 1, . . . , n 2 and t ∈ [t j−1 , t j ], with j = 1, . . . , N and where c l rs are the structure constants of the Lie algebra structure of u * (n). (37) is far from trivial because of the presence of the coefficients {k l j }, which require of the solution of ρ(t) to be determined, since

Implicit numeric integration. The integration of equations
These operators couple the equation for ρ(t) and the equation for H(t), which, otherwise, would be independent. In any case, it is not clear whether the optimal solution exist. Therefore, we have designed an iterative algorithm to obtain a numerical solution for equations (33), always improving the distance to the target points. We consider only one interval [t j−1 , t j ], for the sake of simplicity: • We fix the initial value of K j = K 0 j = 0 (the null matrix) and integrate equation (37). This gives us an initial solution ρ 0 (t) which, in general, will be a bad solution of the interpolation problem from ρ j−1 to ρ j .
• From this solution we determine the corresponding value of K 1 j := 1 [ρ j , ρ 0 (t j )] u * . Notice that this expression corresponds to the extremal of the variations of the cost functional G = 1 2 ρ − ρ j 2 , and therefore it defines the suitable steering acceleration of the solution (ρ(t), H(t)) to make ρ(t j ) to move towards ρ j . For a suitable (not too small) value of , the computation of a new solution of equation (37) with this value of K 1 j , produces a solution (ρ 1 (t), H 1 (t)), where ρ 1 (t) is a better solution of the interpolation problem. If the value of is too small, the acceleration of H 1 (t), even in the correct direction, will be too large and the solution ρ 1 (t) might become a worst interpolation solution than ρ 0 (t).
• From solution ρ 1 (t) we recompute K 2 j := K 1 j + 1 [ρ j , ρ 1 (t j )] u * and integrate again equation (37) to define a new solution (ρ 2 (t), H 2 (t)), where ρ 2 (t) will be a better solution of the interpolation problem. We repeat the process to obtain a sequence of solutions {ρ j (t)}, each one representing a better interpolation solution for ρ(t). Notice that any curve of this series is a solution of equation (33) and therefore provides a solution for the interpolation problem. Each iteration allows us to obtain solutions of equations (33) which are closer and closer to the target points. At the same time, the number of iterations of the algorithm also affects the value of the continuous part of the functional J (equation 22): the larger the number of iterations, the larger value for the integral part. Hence adjusting the number of iterations can also be considered a mechanism to obtain solutions of the interpolation problem which assign different weights to the continuous part with respect to the discrete part of the functional.
• On each of these steps the value of 1 [ρ j , ρ i (t j )] u * decreases, since ρ j − ρ i (t j ) 2 is decreasing. Therefore, after a sufficiently large number r of steps, the value of K r j stabilizes within a certain tolerance. The speed of the process depends on the value of , that is fixing the magnitude of K j (i.e. an acceleration for H(t)) and thus fixes the velocity of von Neumann equation for ρ(t). In any case, notice that (ρ i (t), H i (t)) is not a solution of equation (36) since H i (t) is a solution of an equation obtained with respect to a value of K i j which has been obtained from ρ i−1 (t) as and not as which is the form required for the solution (ρ i (t), H i (t)) to be a solution of equations (37) and (36). Hence we can conclude that it is not the optimal solution for our problem. Nonetheless, it is still a solution of equations (33) and therefore an acceptable one for our interpolation problem, even if more optimal solutions may exist. Furthermore, as we will see later in the examples, the method is very efficient and defines very accurate solutions with a small number of iterations.

Numerical integration: unitary methods.
Our construction has produced the system of equations (37) for the solutions of the optimization problem. It is crucial now to notice that these equations represent the flow of the Hamiltonian system (26), the first one (the x-coordinate in our notation) being the coordinate expression of von Neumann's equation, which was proved to be a Hamiltonian vector field (see equation (11)) with respect to the canonical Lie-Poisson tensor on u * (n). This Hamiltonian vector field represents the infinitesimal generator of the coadjoint action of the unitary group U (n) and therefore we know that its integral curves are contained in the corresponding orbit, which is known to be a symplectic submanifold of u * (n) (see [1]). In particular, if we consider an initial condition which is a pure state, the corresponding symplectic orbit will be a curve on the set D 1 (H) and therefore diffeomorphic to a curve in P H. Notice that, being a unitary transformation, the evolution also defines isometries for the scalar product on u * (n). If we are considering the interpolation of a set of points which do not belong to the same unitary orbit, the situation is a bit different. As the evolution we choose (with the vector field associated to von Neumann equation) is always unitary, any solution of equations (37) must define a unitary transformation. Hence, the minimal distance of the trajectory to the target point corresponds to the distance between the unitary orbit containing the initial condition chosen and the target point. Notice that the algorithm introduced in the previous section makes sense also in this case: the sequence of operators K j defined by equation (38) steer the solutions towards the point in the unitary orbit which minimizes the distance with respect to the target point. We will discuss this point in detail in the following section for a particular example.
From both properties we know that, by implementing a numerical method which preserves the symplectic structure and the distance, we will be ensuring that trajectories of equation (37) will remain on the corresponding orbit and therefore that they will define solutions of our interpolation problem. There exists several numerical methods with those properties, being one of the best known the Gauss(-Legendre) Runge-Kutta (Gauss RK) implicit method (see [14] and [17]), that we will be using in the example of the next Section. The transition matrix of the integrator defines thus a canonical transformation (since the integrator is symplectic) which is also an isometry for the metric structure.

Example 1: the interpolation problem for a qubit
As an application of our method, this section presents the optimal control problem on the Lie algebra u * (2), with the function (20), L, defined on u * (2) × u * (2) × u * (2). Our goal is to determine a quantum spline for a set of points {ρ j } ∈ D 1 (H) and a set of times {t j }.
We will use, for simplicity, the identification of u * (2) with the set of Hermitian matrices iu(2) and consider as basis the set of Hermitian Pauli spin matrices: while the elements of the form H ∈ T ρ iu(2) ∼ iu (2) become With respect to these coordinates, the system (37) turns out to bė Note that the definition of K j implies k 1 j = 0, for all j = 1, . . . , N . The dynamics also implies that the curve ρ in the set D 1 (H) of rank-one-projectors can be identified with a curve on the border of a sphere of radius r = 1 2 , that is, x 2 2 +x 2 3 +x 2 4 = 1 4 . Equations (34) imply that the value of u(t) at the final point of each subinterval must vanish and this is imposed as a boundary condition. The only free parameter is the value of H(0), which can be seen to be irrelevant for the algorithm, since the implicit algorithm adapts itself to it. Remember also that H(t) is required to be continuous.
Consider the uniform partition of times t 0 = 0 < t 1 = 1 5 < t 2 = 2 5 < t 3 = 3 5 < t 4 = 4 5 < t 5 = 1 and the following set of rank-one-projectors: The problem consists in finding the optimal control function u(t) and the optimal state trajectory (ρ(t), H(t)) that minimizes the functional and satisfies the dynamical systeṁ subject to the initial conditions where regularity and interpolation conditions are assumed according to the dynamical interpolation problem P . We have implemented our algorithm described in Section 3.3 for different values of the number of iterations and also different values of the parameter . Notice that both options are not independent since a larger value of leads to a faster convergence for the algorithm. The number of iterations of the algorithm also affects the value of the continuous part of the functional J (equation 43): the larger the number of iterations, the larger the value for the integral part. In the table below we present the values obtained for the distance to the target points but also for the continuous part of the functional. Numerical integration was done with Wolfram's Mathematica, by using the implicit unitary integrator included in the "ImplicitRungeKutta" library of the NDSolveUtilities package. It is immediate to verify that the purity of the state does not change at all because of the unitarity of the integrator (hence the points are always on the surface of the sphere) The resulting curve can be found in the figure 1 and it produces a remarkably accurate result with suitable values for the parameter , as it can be seen in the tables for = 0.005 and 5 iterations (left) and for = 0.005 and 50 iterations (right) We verify that the distance to the target points can be made almost zero for a sufficiently large number of iterations. At each intermediate point t k , we also include the values of the continuous part of the functional (J cont = t k t k−1 1 2 u 2 dt) computed over the subinterval and the total value of J.
The computation time on a PC takes less than two seconds per iteration for the example presented. We will see in next Section how the time increases significantly when considering a three-level case.

Example 2: the interpolation problem for a qutrit
As a second application of our method, this section presents the case of a threelevel systems (a qutrit). The formulation is entirely analogous to the previous case but now the problem is formulated on the space D(C 3 ) ⊂ u * (3). As we did in the qubit case we will use the identification of u * (3) with iu(3) and we will represent the density states as 3 × 3 Hermitian matrices. As a basis, we will consider the set of Gell-Mann matrices The structure constants with respect to this basis are while the elements of the form H ∈ T ρ iu(3) ∼ iu (3) become With respect to these coordinates, we consider now the problem of finding the optimal control function u(t) and the optimal state trajectory (ρ(t), H(t)) that minimizes the functional and satisfies the dynamical systeṁ subject to the initial conditions where regularity and interpolation conditions are assumed according to the dynamical interpolation problem P . As we saw above, we have to solve the system (37). Again, equations (34) imply that the value of u(t) at the final point of each subinterval must vanish and this is imposed as a boundary condition. The value of H(0) can also be seen to be irrelevant for the algorithm, as it happens in the case of two levels.
5.1. First example: points contained in a unitary orbit. Let us consider the following set of points, for times t 0 = 0 < t 1 = 1 6 < t 2 = 1 3 < t 3 = 1 2 < t 4 = 2 3 < t 5 = 5 6 < t 6 = 1: which are at a distance of 0.001 of ρ 1 and ρ 2 , respectively. These are the closest possible points since the unitary transformation must preserve the distance between both unitary orbits which is 0.001 at the initial points (ρ 0 and ρ 0b ). We conclude thus that our method provides an efficient numerical solution for the interpolation problem of general (pure or mixed) quantum states and generates a trajectory contained in the closest unitary orbit to them.

Conclusions and outlook
In this paper we have generalized the notion of quantum spline introduced in [9] to the case of general quantum states, pure or mixed. Brody et al. formulated a variational problem on a projective space defined in Section 1. They considered a globally defined continuous curve on the Lie group U (n) and the corresponding set of continuous variations, and study the set of extremals of the functional defined in equation (1), which combines the curve on U (n) and its action on the complex projective space. From this, they obtain a set of differential equations whose solutions correspond to Riemannian cubic polynomials on the subintervals ]t j−1 , t j [ and a set of boundary conditions on the times {t j }. Finding a global solution is a difficult problem and they also provide a numerical algorithm illustrated by a particular example in the case of a two level system.
Our definition reconsiders the problem and writes it on the complete set of quantum states (pure and mixed), modeled as a submanifold D of the dual space u * (n). Considering (as in [9] ) the case of unitary dynamics, allows us to forget about D and define the problem globally on the linear space u * (n). This simplifies the problem in a remarkable way. Furthermore, we consider a different type of optimization problem and choose a local formulation on the different subintervals of the time domain which does not impose the differentiability of the time dependence of the Hamiltonian. This choice, which is common, for instance, in quantum control solutions based on pulses, simplifies further the problem from the mathematical point of view, and transform it in a chain of Bolza-type problems. We consider then a Hamiltonian formulation of the corresponding control problem based on Pontryagin Maximum Principle, which allows us to consider in a natural way symplectic integrators associated to the geometrical formulation of Quantum Mechanics. Despite the much simpler formulation, the resulting system of differential equations is still difficult to solve exactly. We have introduced then an iterative algorithm to determine good solutions of the problem which works very efficiently for pure states and for mixed states, even for problems where the points to be interpolated do not belong to the same unitary orbit. Our algorithm allows us to define the closest unitary orbit in an efficient way.
Another advantage of our formalism is that it admits further generalizations in a simple way. Indeed, we can consider an analogous formulation where unitary dynamics is replaced by more general master equations, as for instance the Lindblad-Kossakowski equation. This Markovian generalization may impose some extra constraints on the set of possible times and points to be interpolated, but the formulation of the problem makes perfect sense. The conclusions of that analysis, which is being done now, will be presented in a future paper.