Time-Minimal Control of Dissipative Two-level Quantum Systems: the Generic Case

# Time-Minimal Control of Dissipative Two-level Quantum Systems: the Generic Case

Bernard Bonnard,  Monique Chyba  and Dominique Sugny B. Bonnard is with the Institut de Mathématiques de Bourgogne, UMR CNRS 5584, 9 Avenue Alain Savary, BP 47 870 F-21078 DIJON Cedex FRANCE (bernard.bonnard@u-bourgogne.fr).M. Chyba is with the University of Hawaii, Department of Mathematics, Honolulu, HI 96822, USA (mchyba@math.hawaii.edu).D. Sugny is with the Institut Carnot de Bourgogne, UMR 5209 CNRS-Université de Bourgogne, 9 Av. A. Savary, BP 47 870, F-21078 DIJON Cedex, FRANCE (dominique.sugny@u-bourgogne.fr).
###### Abstract

The objective of this article is to complete preliminary results from [4, 13] concerning the time-minimal control of dissipative two-level quantum systems whose dynamics is governed by Lindblad equations. The extremal system is described by a 3D-Hamiltonian 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 time-minimal control analysis of two-level dissipative quantum systems whose dynamics is governed by Lindblad equation. More generally, according to [10], the dynamics of a finite-dimensional quantum system in contact with a dissipative environment is described by the evolution of the density matrix given by the equation

 i∂ρ∂t=[H0+H1,ρ]+iL(ρ), (1)

where is the field-free 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:

 ˙ρnn = −i[H0+H1,ρ]nn−∑k≠nγknρnn+∑k≠nγnkρkk ˙ρkn = −i[H0+H1,ρ]kn−Γknρkn, k≠n (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

 H1=−μxEx−μyEy,

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

 i∂∂t⎛⎜ ⎜ ⎜ ⎜⎝ρ11ρ12ρ21ρ22⎞⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝−iγ12−E∗Eiγ21−E−ω−iΓ0EE∗0ω−iΓ−E∗iγ12E∗−E−iγ21⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜⎝ρ11ρ12ρ21ρ22⎞⎟ ⎟ ⎟ ⎟⎠, (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

 i∂∂t⎛⎜ ⎜ ⎜ ⎜⎝ρ11ρ12ρ21ρ22⎞⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜⎝−iγ12−u∗uiγ21−u−iΓ0uu∗0−iΓ−u∗iγ12u∗−u−iγ21⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜⎝ρ11ρ12ρ21ρ22⎞⎟ ⎟ ⎟ ⎟⎠. (4)

The interaction representation means that we have transformed the mixed-state 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 time-minimal transfer problem from a state to a state . Hence, we have to analyze a time-minimal control problem for a bilinear system of the form:

 ˙q=F0(q)+2∑i=1uiFi(q), |u|≤1,

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 time-minimal synthesis of [13] corresponding to the case where is real but introducing more general tools to handle the problem. It corresponds to a time-minimal control problem of a two-dimensional 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 semi-direct Lie group. This classification is physically relevant to analyze the 3D-case because, using the symmetry of revolution of the problem, it gives the time-minimal 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 two-sphere 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 2D-single 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 2D-case

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 single-input system:

 ˙y = −Γy−u1z ˙z = γ−−γ+z+u1y,|u1|≤1.

It gives the time-optimal 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 time-minimal control problem, this amounts to exchange the trajectories and corresponding respectively to and . Also, according to this property, the time-minimal 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:

 ˙q=F(q)+uG(q)

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:

 G=(−zy), F=(−Γyγ−−γ+z), [G,F]=((γ+−Γ)z−γ−(γ+−Γ)y).

Hence, is given by:

 y[2(Γ−γ+)z+γ−]=0

and if , the singular set is defined by the two lines and . The collinear set is defined by:

 γ+z2+Γy2−γ−z=0

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

 C1: y=γ−1+γ+Γ, z=−Γy.

The collinear set shrinks into a point when . Computing the intersection of with the singular line , one gets:

 Γy2=γ2−4(γ+−Γ)2(γ+−2Γ).

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 Legendre-Clebsch condition, it is splitted into fast and slow directions. In the 2D-case, it is tested by Lie brackets configurations and can be computed by introducing:

 D=det(G,[[G,F],G]), D′′=det(G,F).

The trajectory is time-optimal in the so-called hyperbolic case and time-maximal in the so-called elliptic case . Computing Lie brackets of length 3, we have:

 [[G,F],G]=(2(γ+−Γ)yγ−−2(γ+−Γ)z)

and

 [[G,F],F]=((γ+−Γ)(γ−−γ+z)+Γ[(γ+−Γ)z−γ−](γ+−Γ)2y).

Hence:

 D′′=∣∣∣−z−Γyyγ−−γ+z∣∣∣, D=∣∣∣−z2(γ+−Γ)yyγ−2(γ+−Γ)z∣∣∣.

For the singular direction , we get:

 DD′′=2z2(γ+−Γ)γ+(z−γ−2(γ+−Γ))(z−γ−γ+).

Near the origin, the sign is always positive if . If , the sign is given by . For the singular direction , we have:

 DD′′=y22(γ+−Γ)[γ2−(γ+−2Γ)−4Γy2(γ+−Γ)2].

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 time-optimal 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:

 y=ρcosϕ, z=ρsinϕ,

one gets:

 ˙ρ = γ−sinϕ+ρ[−Γ+(Γ−γ+)sin2ϕ] ˙ϕ = γ−cosϕρ+(Γ−γ+)sin(2ϕ)2+u1.

If we use the coordinates and , the system becomes:

 ˙x = −2Γx+γ−z+z2(Γ−γ+) ˙z = γ−−γ+z+u1√2x−z2.

Making a feedback transformation of the form where is a function of and , we can consider the system:

 ˙x = −2Γx+γ−z+z2(Γ−γ+) ˙z = γ−−γ+z+u1.

If we set where , we obtain the system:

 ˙x = γ2−4(γ+−Γ)−2Γx+(Γ−γ+)z2 ˙z = γ−(γ+−2Γ)2(γ+−Γ)−γ+z+u1.

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:

 ˙px = 2Γpx ˙pz = −2zpx(Γ−γ+)+pzγ−,

and can be easily integrated to compute the time-minimal 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:

 us=−D′D=γ−(γ+−2Γ)2y(Γ−γ+),

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:

 ˙p=−p(∂F∂q+u∂G∂q), u=±1.

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:

 ˙v=(∂F∂q+u∂G∂q)v, u=±1,

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:

 p(0)G(q(0))=p(t)G(q(t))=0.

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 ad-formulae 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 semi-direct Lie group identified to the set of matrices of :

 (10gv), g∈GL(2,R), v∈R2,

acting on the subspace of vectors in : .

Lie brackets computations are defined as follows. We set:

 F(q)=Aq+a, G(q)=Bq,

and are identified to , in the Lie algebra . The Lie brackets computations on the semi-direct product are defined by:

 [(A′,a′),(B′,b′)]=([A′,B′],A′b′−B′a′).

We now compute . The first step consists in determining which amounts to compute . We write where is the center

 R(1001).

We choose the following basis of :

 B=(0−110), C=(0110) and D=(100−1).

The matrix is decomposed into:

 A=(−Γ00−γ+)=(λ00λ)+(s00−s)

and hence and . In the basis (, , ), is represented by the matrix:

 ⎛⎜⎝0−2s0−2s02ε0−2ε0⎞⎟⎠.

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

 B=αv0+βv1+γv2,

with:

 α=εε2−s2, β=−λ2s2(λ2−λ1)(ε2−s2), γ=−λ1s2(λ1−λ2)(ε2−s2).

Hence one gets:

To test the collinearity at , we compute

where the determinant is equal to

 (z20−y20)(αs+2ε(βe−λ1t+γe−λ2t))+2y0z0(λ1βe−λ1t+λ2γe−λ2t).

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:

 diag(0,(0θ−θ0)).

Hence, we have in this basis:

We decompose in the same basis:

 B=αv0+βv1+γv2,

where

 α=εε2−s2, β=−s2(ε2−s2), γ=0.

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:

 [(A′,a′),(B′,b′)]=([A′,B′],A′b′−B′a′).

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 time-minimal control problem in the full control case, replacing the control domain by .

### 2.6 The time-minimal synthesis

We use [3] as general reference on time-minimal synthesis, see also [7]. The initial condition is fixed to and we consider the problem of constructing the time-minimal 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 fish-shaped accessibility set represented on Fig. 3. This set is not closed since the arc starting from is not in .

We next list the micro-local 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 bang-bang 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 non-admissible singular direction. Every extremal curve is bang-bang 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 time-optimal syntheses in the single-input 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 bi-input case

### 3.1 Geometric analysis

The system is written in short in Cartesian coordinates as follows:

 ˙q=F0(q)+u1F1(q)+u2F2(q), |u|≤1.

Introducing the Hamiltonians , the pseudo-Hamiltonian associated to the time-optimal control problem is:

 H=H0+2∑i=1uiHi+p0,

where . The time-optimal control is given outside the switching surface : , by , , with the corresponding true Hamiltonian:

 Hr=H0+√H21+H22,

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:

 x=ρsinϕcosθ, y=ρsinϕsinθ, z=ρcosϕ

and a feedback transformation:

 v1=u1cosθ+u2sinθ, v2=−u1sinθ+u2cosθ.

We obtain the system:

 ˙ρ = γ−cosϕ−ρ(γ+cos2ϕ+Γsin2ϕ) (6a) ˙ϕ = −γ−sinϕρ+sin(2ϕ)2(γ+−Γ)+v2 (6b) ˙θ = −cotϕv1. (6c)

Hence, one deduces that the true Hamiltonian is:

 Hr=[γ−cosϕ−ρ(γ+cos2ϕ+Γsin2ϕ)]pρ+pϕ[−γ−sinϕρ+sin(2ϕ)2(γ+−Γ)]+√p2ϕ+p2θcot2ϕ.

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:

 Hr=−(γ+cos2ϕ+Γsin2ϕ)pr+sin(2ϕ)(γ+−Γ)pϕ+√p2ϕ+p2θcot2ϕ. (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 2D-system, we generate the extremal synthesis for the 3D-system.

More generally, we have for a system on the two-sphere 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 two-sphere 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 :

 ˙ρ = γ−cosϕ−ρ(γ+cos2ϕ+Γsin2ϕ) ˙ϕ = −γ−sinϕρ+sin(2ϕ)2(γ+−Γ)+pϕQ (8) ˙θ = pθcot2ϕQ,

and

 ˙pρ = (γ+cos2ϕ+Γsin2ϕ)pρ−γ−sinϕρ2pϕ ˙pϕ = [γ−sinϕ+ρ(Γ−γ+)sin(2ϕ)]pρ−[−1ρcosϕγ−+(γ+−Γ)cos(2ϕ)]pϕ+p2θcosϕQsin3ϕ ˙pθ = 0,

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 time-optimal solution of the 2D-system, 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 2D-system. 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:

 ˙ρ = γ−ψ−ρ[Γ+(γ+−