TimeMinimal Control of Dissipative Twolevel Quantum Systems: the Generic Case
Abstract
The objective of this article is to complete preliminary results from [4, 13] concerning the timeminimal control of dissipative twolevel quantum systems whose dynamics is governed by Lindblad equations. The extremal system is described by a 3DHamiltonian depending upon three parameters. We combine geometric techniques with numerical simulations to deduce the optimal solutions.
Keywords. Time optimal control, conjugate and cut
loci, quantum control
1 Introduction
In this article, we consider the timeminimal control analysis of twolevel dissipative quantum systems whose dynamics is governed by Lindblad equation. More generally, according to [10], the dynamics of a finitedimensional quantum system in contact with a dissipative environment is described by the evolution of the density matrix given by the equation
(1) 
where is the fieldfree Hamiltonian of the system, represents the interaction with the control field and the dissipative part of the equation; is the commutator of the operators and defined by . Equation (1) is written in units such that . The components of the density matrix satisfy the following equations:
(2) 
where and for a level quantum system. The parameters describe the population relaxation from state to state whereas is the dephasing rate of the transition from state to state . Note that not every positive parameter or is acceptable from a physical point of view since the density matrix must satisfy particular properties (trace conservation, hermitian operator and complete positivity [1, 9, 10]).
Particularizing now to the case , we assume that is of the form
where the operators and are proportional to the Pauli matrices and in the eigenbasis of . The electric field is the superposition of two linearly polarized fields and and we assume that these two fields are in resonance with the Bohr frequency . In the RWA approximation, the time evolution of satisfies the following form of the Lindblad equation
(3) 
where is equal to and is the complex Rabi frequency of the laser field (the real and imaginary parts of are the amplitudes of the real fields and up to a multiplicative constant). In Eq. (3), is the difference of energy between the ground and excited states of the system. In the interaction representation, Eq. (3) becomes
(4) 
The interaction representation means that we have transformed the mixedstate with the unitary transformation . Since , the density matrix can be represented by the vector where , and and belongs to the Bloch ball . Equation (4) takes the form:
(5) 
is the set of parameters such that and and they satisfy the following inequations derived from the Lindblad equation [12], the Bloch ball being invariant. The control is where and are two real functions. We can write the control field where and up to a rescaling of the time and dissipative parameters we can assume that .
We consider the timeminimal transfer problem from a state to a state . Hence, we have to analyze a timeminimal control problem for a bilinear system of the form:
where the drift term depends upon three parameters. This problem is a very difficult problem whose analysis requires advanced mathematical tools from geometric control theory and numerical simulations.
Such analysis is motivated by physical reasons. It is a fundamental model in quantum control and more general problems can be handled by coupling such systems. Numerous optimal control results exist in the conservative case e.g. [6], but only partial ones for this problem: a pioneering work [13] assuming real and a second one [4] for complex but restricted to .
The objective of this article is double. First of all, we introduce all the proper geometric tools to analyze the problem using Pontryagin maximum principle. Secondly, our aim is to make a complete study for every generic parameter in , combining mathematical reasonings and numerical simulations based on shooting techniques and including computations of conjugate points to test optimality. Based on the Cotcot code [2], they can be used in practice to compute the true optimal control, once the physical parameters are identified.
The organization of this article is the following. In section 2, we complete the classification of the timeminimal synthesis of [13] corresponding to the case where is real but introducing more general tools to handle the problem. It corresponds to a timeminimal control problem of a twodimensional bilinear system in the single input case. The optimal synthesis for a fixed initial point can be constructed by gluing together local optimal syntheses. We can also make estimates of switching points by lifting the system on a semidirect Lie group. This classification is physically relevant to analyze the 3Dcase because, using the symmetry of revolution of the problem, it gives the timeminimal synthesis for initial points . This geometric property is explained in section 3. Moreover, using spherical coordinates the system can be viewed as a system on a twosphere of revolution coupled with the evolution of the distance to the origin, which represents the purity of the system. According to the maximum principle, smooth extremals are solutions of the Hamiltonian vector field where , , and the control components are given by , . Non smooth extremals can be constructed by connecting smooth subarcs of the switching surface : . A contribution of this article in section 3 is to classify the possible connections. We proved that every non smooth extremal is either a solution of the 2Dsingle input system, assuming real, or occurs when meeting the equatorial plane of the Bloch ball. In the second case, the switching can be handled numerically using an integrator with an adaptative step. In the same section, we combine analytical and numerical analysis to determine the extremals and compute conjugate points. This completes the analysis from [4] in the integrable case. The physical interpretation is presented as a conclusion.
2 The 2Dcase
Following [13], a first step in the analysis is to consider the following reduced system. Assuming real, the coordinate is not controllable and we can consider the planar singleinput system:
It gives the timeoptimal analysis of the control problem when the initial state is a pure state on the axis, that is . We proceed as follows to make the analysis.
2.1 Symmetry group
A discrete symmetry group is associated to reflections with respect to the different axes. More precisely, if then one gets the same system changing into and into . If then one gets the same system changing into . In particular, concerning the timeminimal control problem, this amounts to exchange the trajectories and corresponding respectively to and . Also, according to this property, the timeminimal synthesis is symmetric with respect to the axis. This is connected to the symmetry of revolution of the whole system which is explained later.
2.2 The feedback classification
A preliminary step in our analysis is to consider the feedback classification problem [3]. The system is written in a more compact form as follows:
where and are affine vector fields. To make the feedback classification, we relax the control bound . The geometric invariants are related to the sets:

where are located the singular trajectories.

corresponding to the set of points where and are collinear.
They are obtained by straightforward computations of Lie brackets:
Hence, is given by:
and if , the singular set is defined by the two lines and . The collinear set is defined by:
and is a closed curve formed by the union of two arcs of parabola, containing and , with , which is the equilibrium state of the free motion. More generally, contains the equilibrium points for the dynamics with constant control , since . In particular, the equilibrium point when is given by
The collinear set shrinks into a point when . Computing the intersection of with the singular line , one gets:
If , since , one deduces that the intersection is empty except in the case where . An important consequence is to simplify the classification of the optimal syntheses. We represent on Fig. 1 the sets and for a situation with and .
Another feedback invariant is the optimality status of singular trajectories. Using the generalized LegendreClebsch condition, it is splitted into fast and slow directions. In the 2Dcase, it is tested by Lie brackets configurations and can be computed by introducing:
The trajectory is timeoptimal in the socalled hyperbolic case and timemaximal in the socalled elliptic case . Computing Lie brackets of length 3, we have:
and
Hence:
For the singular direction , we get:
Near the origin, the sign is always positive if . If , the sign is given by . For the singular direction , we have:
Hence, near the origin, the optimality is given by the sign of . Moreover, since , one gets that if and if .
2.3 Computation of a normal form
In order to analyze the timeoptimal control, we compute a global normal form up to a polar singularity for the action of the feedback group. A first step is to linearize the vector field since it is connected to the evaluation of the switching times. One further normalization is to straight the horizontal singular line: . Using polar coordinates:
one gets:
If we use the coordinates and , the system becomes:
Making a feedback transformation of the form where is a function of and , we can consider the system:
If we set where , we obtain the system:
In this simplified model where the control is rescaled by the positive parameter , we keep most of the information about the initial system. In particular, all the feedback invariants are preserved: the collinear set corresponds to and the singular set is identified to with its optimality status, while we have wiped out the singular line . Due to the feedback transformation, we lose however the singularities of the vector fields and and also the saturation phenomenon of the singular control.
For the simplified model, the adjoint system takes the form:
and can be easily integrated to compute the timeminimal synthesis.
2.4 The saturation phenomenon
One interesting property which is not captured by the normal form is when the singular control is saturating along the horizontal singular line . Introducing , we get on the singular line: . The singular control is given by:
and saturation occurs when . Observe that if then the singular control is never admissible when .
2.5 The switching function
For the true system, the switching function is more intricate but can be still analyzed using geometric and numerical techniques. The switching function is given by and switching occurs when . By construction, is tangent to the circle , hence is rotating when we follow an arc curve or . The dynamics of is given by the adjoint equation:
The corresponding dynamics is linear and can be either oscillating if the eigenvalues are complex or non oscillating if they are real.
An equivalent but more geometric test is the use of the standard function introduced in [7] and defined as follows. Let be the tangent vector solution of the variational equation:
whose dynamics is similar to the one of the adjoint vector. Let be two consecutive switching times on an arc or . One can set and . By definition, we have:
We denote by the solution of the variational equation such that and where this equation is integrated backwards from time to time 0. By construction and we deduce that at time 0, is orthogonal to and to . Therefore, and are collinear; is defined as the angle between and measured counterclockwise. One deduces that switching occurs when . In the analytic case, can be computed using Lie brackets. Indeed, for , , we have by definition
and in the analytic case, the adformulae gives :
Here, to make the computation explicit, we take advantage of the fact that we can lift our bilinear system into an invariant system onto the semidirect Lie group identified to the set of matrices of :
acting on the subspace of vectors in : .
Lie brackets computations are defined as follows. We set:
and are identified to , in the Lie algebra . The Lie brackets computations on the semidirect product are defined by:
We now compute . The first step consists in determining which amounts to compute . We write where is the center
We choose the following basis of :
The matrix is decomposed into:
and hence and . In the basis (, , ), is represented by the matrix:
The characteristic polynomial is
and the
eigenvalues are and , ; and
are distinct and real if and we note
, ;
and are distinct and imaginary if
and we note
, . To
compute , we must distinct two
cases.
Real case: In the basis , , , the eigenvectors
corresponding to are respectively:
,
and . Therefore, in this
eigenvector basis, is the
diagonal matrix: . To compute , we use the decomposition
with:
Hence one gets:
To test the collinearity at , we compute
where the determinant is equal to
Using the fact that this last expression has at most two zeros,
one being for , it is straightforward to check that for
, this expression only vanishes at . This proves
the result numerically checked in [13] that there is
at most one switching. This analysis can be generalized to
any initial condition.
Imaginary case: In this case, we note
the eigenvalue associated to the eigenvector
. We consider the real part
and the imaginary part
. In the basis ,
, , takes the normal
form:
Hence, we have in this basis:
We decompose in the same basis:
where
Hence, we get:
Computing we obtain:
which vanishes for two values of in if
. These two values coincide for . Hence
the switchings occur periodically with a period of
which confirms the numerical simulations of
[13].
Case : The Lie bracket is given by:
We note the canonical basis and . We have:
More generally, one gets:
where is given by the recurrence relation:
The computation is intricate but simplifies if
since in this case for
. Numerical simulations have to be used to compute the switchings sequence.
Generalization: This technique can be generalized to the
timeminimal control problem in the full control case, replacing
the control domain by .
2.6 The timeminimal synthesis
We use [3] as general reference on timeminimal synthesis, see also [7]. The initial condition is fixed to and we consider the problem of constructing the timeminimal synthesis from this initial point. This amounts to compute two objects:

The switching locus of optimal trajectories which is deduced from the switching locus of extremal trajectories.

The cut locus which is formed by the set of points where a minimizer ceases to be optimal.
In order to achieve this task, we must glue together local time minimal syntheses which are classified. To be more precise, we recover the case (d) of [13], the gluing being indicated on Fig. 2 on which we have represented the local extremal classifications of [3] which are crucial to deduce the optimal syntheses.
In this case, the cut locus is a segment of the axis and its birth is located at the initial point which is a consequence of the elliptic situation. The switching locus is the union of a fast singular trajectory, corresponding to an hyperbolic point and a curve corresponding to parabolic points (for the terminology see [3]). Note also the importance of the tangential point where arcs and are tangent leading to the fishshaped accessibility set represented on Fig. 3. This set is not closed since the arc starting from is not in .
We next list the microlocal situations we need to construct the synthesis.

Ordinary switching points: The local synthesis is given by or . The two cases are distinguished using for instance the clock form with and which is also useful to get more global results.

Hyperbolic case: Existence of a fast singular admissible trajectory. The optimal synthesis is of the form where is a singular arc (see Fig. 4).

Elliptic case: Existence of a slow admissible singular trajectory. An optimal arc is bangbang with at most one switching. Not every extremal trajectory is optimal and we have birth of a cut locus (see Fig. 5).

Parabolic point: It corresponds to a nonadmissible singular direction. Every extremal curve is bangbang with at most two switchings. In our case, the initial point is fixed and the switching locus starts with the intersection of with the singular line (see Fig. 6).

Saturating case:
A fast singular trajectory is saturating at a point : birth of a switching curve at (see Fig. 7). 
A case:
A fast singular trajectory meets the set and becomes slow.
2.7 Classification of the optimal syntheses
We describe the different timeoptimal syntheses in the singleinput case. Without loss of generality, we restrict the study to the initial point . The classification is done with respect to the relative positions of the feedback invariants and and to the optimal status of singular extremals which are fast or slow according to the values of , and .
For , the set is restricted to the origin and we have two cases according to the sign of . Note that the form of the extremals and starting from depends on the sign of . Two cases for are presented in [13]. We complete this study with the optimal synthesis for and displayed in Fig. 9a.
For , we distinguish four cases according to the signs of and . One case ( and ) is treated in [13]. We consider here three types of optimal synthesis represented in Fig. 9b, 9c and 9d. Note that in a same class of synthesis the reachable set from the initial point depends on the dissipative parameters which can modify the structure of the synthesis. The last case and can be deduced from the case and since the horizontal singular line plays no role in both cases. The synthesis of Fig. 9d is very similar to the one of Fig. 2 except the fact that a part of the horizontal singular line is admissible. The switching locus has been computed numerically using the switching function .
The role of the parameter is clearly illustrated in Figs. 9a and 9c. The case is a degenerate case where the set shrinks into a point. The variation of induces a bifurcation of the control system leading to new structures of the optimal synthesis. For , the set is the union of two branches of parabola. The optimal status of the vertical singular line changes when this line crosses the set in Fig. 9c.
3 The biinput case
3.1 Geometric analysis
The system is written in short in Cartesian coordinates as follows:
Introducing the Hamiltonians , the pseudoHamiltonian associated to the timeoptimal control problem is:
where . The timeoptimal control is given outside the switching surface : , by , , with the corresponding true Hamiltonian:
whose solutions (outside ) are smooth and are called extremals of order 0. More general non smooth extremals can be obtained by connecting such arcs through .
To make the geometric analysis and to highlight the symmetry of revolution, the system is written using the spherical coordinates:
and a feedback transformation:
We obtain the system:
(6a)  
(6b)  
(6c) 
Hence, one deduces that the true Hamiltonian is:
From this, we deduce the following lemma:
Lemma 1.
(i) The angle is a cyclic variable and is a
first integral (symmetry of revolution).
(ii) For , using the coordinate , the
Hamiltonian takes the form:
(7) 
Hence, is an additional cyclic variable and is a first integral. The system is thus Liouville integrable.
As a consequence, we can deduce two properties. First of all, the axis is an axis of revolution and the state is a pole. This means that by making a rotation around () of the extremal synthesis for the 2Dsystem, we generate the extremal synthesis for the 3Dsystem.
More generally, we have for a system on the twosphere of revolution described by Eqs. (6b) and (6c) coupled with the one dimensional system (6a) describing the evolution of the physical variable corresponding to the purity of the system. Moreover, the system is invariant for the transformation which is associated to a reflexional symmetry with respect to the equator for the system (6) restricted to the twosphere of revolution. This property is crucial in the analysis of the integrable case.
If then the situation is more intricate. The extremals solutions of order 0 satisfy the equations which are singular for :
(8)  
and
where .
3.2 Regularity analysis
The smooth extremal curves solutions of are not the only extremals because more complicated behaviors are due to the existence of the switching surface : . Hence, in order to get singularity results, we must analyze the possible connections of two smooth extremals crossing to generate a piecewise smooth extremal. This can also generate complex singularities of the Fuller type, where the switching times accumulate. In our problem, the situation is less complex because of the symmetry of revolution. The aim of this section is to make the singularity analysis of the extremals near .
The structure of optimal trajectories is described by the following proposition.
Proposition 1.
Every optimal trajectory is:

Either an extremal trajectory with contained in a meridian plane and timeoptimal solution of the 2Dsystem, where .

Or subarcs solutions of , where with possible connections in the equator plane for which .
Proof.
The first assertion is clear. If then extremals are such that and up to a rotation around the axis, they correspond to solutions of the 2Dsystem. The switching surface is defined by: . We cannot connect an extremal with to an extremal where since at the connection the adjoint vector has to be continuous. Hence, the only remaining possibility is to connect subarcs of with at a point of leading to the conditions and . ∎
Further work is necessary to analyze the behaviors of such
extremals near .
Normal form: A first step in the analysis is to construct
a normal form. Taking the system in spherical coordinates and
setting , the approximation is: