Time-minimum control of quantum purity for 2-level Lindblad equations

TL;DR

This paper investigates the problem of minimizing the time to reach maximum purity in a two-level quantum system governed by Lindblad dynamics, using optimal control techniques including bang-bang and singular controls.

Contribution

It provides a detailed analysis of time-minimum control strategies for quantum purity, including numerical methods and the construction of optimal solutions with singular arcs.

Findings

Multiple local minima exist with unbounded controls.

Optimal controls are bang-bang with bounded controls.

Singular arcs can be part of optimal solutions.

Abstract

We study time-minimum optimal control for a class of quantum two-dimensional dissipative systems whose dynamics are governed by the Lindblad equation and where control inputs acts only in the Hamiltonian. The dynamics of the control system are analyzed as a bi-linear control system on the Bloch ball after a decoupling of such dynamics into intra- and inter-unitary orbits. The (singular) control problem consists of finding a trajectory of the state variables solving a radial equation in the minimum amount of time, starting at the completely mixed state and ending at the state with the maximum achievable purity. The boundary value problem determined by the time-minimum singular optimal control problem is studied numerically. If controls are unbounded, simulations show that multiple local minimal solutions might exist. To find the unique globally minimal solution, we must repeat the algorithm for various initial conditions and find the best solution out of all of the candidates. If controls are bounded, optimal controls are given by bang-bang controls using the Pontryagin minimum principle. Using a switching map we construct optimal solutions consisting of singular arcs. If controls are bounded, the analysis of our model also implies classical analysis done previously for this problem.

Figures

Tables

Table 1. Table 1. Numerical results from time-minimal controls with solutions of various orders.

M	Time	Energy
1	1.9372	7.5830
3	1.9366	8.6873
5	1.9361	1.6368
7	1.9359	1.3765

Table 2. Table 2. Switching times for the example in Figure 3 .

Initial control	Time until first switch	Time between first and second switches
$u=1$	2.1943	0.4685
$u=-1$	1.1532	0.4905

Equations

$$\frac{d}{dt}|\psi (t)⟩=-iH(t)|\psi (t)⟩,$$

$$\frac{d}{dt}\rho =[-iH,\rho ].$$

$$\frac{d}{dt}\rho =[-iH,\rho ]+\sum _{j=1}^N\left(L_j\rho L_j^{†}-\frac{1}{2}\left\{L_j^{†}{L}_j,\rho \right\}\right)$$

$$\frac{d\rho }{dt}=[-iH,\rho ]+\sum _{j=1}^N{L}_j\rho L_j^{†}-\frac{1}{2}\{L_j^{†}{L}_j,\rho \},$$

$$\frac{d\vec{q}}{dt}=\vec{b}+(A-\text{tr}(A))\vec{q}+\vec{u}\times \vec{q},whereA:=\frac{1}{2}\sum _j\overrightarrow{l_j}{\overline{{\vec{l}}_j}}^T+\overline{{\vec{l}}_j}{\vec{l}}_j^T,and\vec{b}:=i\sum _j{\vec{l}}_j\times \overline{{\vec{l}}_j},$$

$$\begin{array}{rl}2r\dot{r}& =2⟨\dot{q},q⟩=2⟨b+u\times q+Bq,q⟩=2⟨b,q⟩+2⟨Bq,q⟩=2r⟨b,\hat{q}⟩+2r^2⟨B\hat{q},\hat{q}⟩\end{array}$$

$$r\frac{dr}{dt}=⟨q,b+Bq⟩:=f(q).$$

$$f_q(r)=⟨q,Bq⟩r^2+⟨q,b⟩r.$$

$$g(q):=-\frac{⟨q,b⟩}{⟨q,Bq⟩}.$$

$$q_{apogee}:=\arg \underset{q\in \mathcal{M}}{\max }\Vert q\Vert ,$$

$$\frac{dr}{dt}=\frac{f(q)}{r}.$$

$$q_0=\epsilon \vec{b}/\Vert b\Vert ,q_f=(1-\delta )q_{apogee}$$

$$\begin{array}[]{rcl}\dot{x}&=&b_{1}+\alpha_{1}x-uy\\ \dot{y}&=&b_{2}+\alpha_{2}y+ux\end{array}$$

$$\frac{d}{d\theta}g(\theta)\Big{|}_{\theta=\theta^{*}}=\frac{-b_{1}b_{2}(\alpha_{1}-\alpha_{2})\sqrt{(b_{1}/\alpha_{1})^{2}+(b_{1}/\alpha_{2})^{2}}}{\alpha_{2}b_{1}^{2}+\alpha_{1}b_{2}^{2}}.$$

$$\frac{d^2}{d{\theta }^2}g({\theta }^{\ast })=\frac{-\sqrt{b_1^2+b_2^2}}{|{\alpha }_1|}<0,$$

$$\frac{d^2}{d{\theta }^2}g({\theta }^{\ast })=\frac{-b_2^2|{\alpha }_1{\alpha }_2|(2{\alpha }_1-{\alpha }_2)}{{\alpha }_2^3|{\alpha }_1{b}_2|}$$

$$q^{*}=\frac{1}{u^{2}+\alpha_{1}\alpha_{2}}\left\{\begin{array}[]{c}-\alpha_{2}b_{1}-b_{2}u\\ -\alpha_{1}b_{2}+b_{1}u\end{array}\right\}.$$

$$\Vert b{\Vert }^2{u}^3+3b_1{b}_2({\alpha }_2-{\alpha }_1)u^2+(2{\alpha }_2^2{b}_1^2+2{\alpha }_1^2{b}_2^2-{\alpha }_1{\alpha }_2(b_2-b_1))u+{\alpha }_1{\alpha }_2{b}_1{b}_2({\alpha }_1-{\alpha }_2)=0.$$

$$\frac{-a+b-c+d}{a+b+c+d}\le 0$$

$$a+b+c+d=0.$$

$$a+b+c+d=0$$

$$\frac{-a+b-c+d}{a+b+c+d}\le 0$$

$$\frac{-a+b-c+d}{a+b+c+d}>0,\text{and}\min |u_{\pm }|<1$$

$$u_{\pm }=\frac{-2b\pm \sqrt{\Delta }}{6a}.$$

$$\min {\int }_0^{t_f}dt=\min {\int }_{x_0}^{x_f}\frac{dt}{dr}\frac{dr}{dx}dx,$$

$$I(q,q^{\prime })=\int L(q,q^{\prime })dt={\int }_{x_0}^{x_f}\frac{x+y{y}^{\prime }}{⟨q,b+Bq⟩}dx.$$

$$y(x)=y^0(x)+\sum _{i=1}^M{c}_i{y}^i(x),$$

$$y^i(x)=(x-x_0)(x-x_f{)}^i,y^0(x)=\frac{y_f-y_0}{x_f-x_0}(x-x_0)+y_0,i=1,\dots ,M$$

$$I_{c_i}={\int }_{x_0}^{x_f}\frac{\partial }{\partial c_i}Ldx=0.$$

$$\nu (c)={\left(\sum _{i=1}^M{\left(I_{c_i}\right)}^2\right)}^{1/2}.$$

Full text

The boundary value problem determined by the time-minimum singular optimal control problem is studied numerically. If controls are unbounded, simulations show that multiple local minimal solutions might exist. To find the unique globally minimal solution, we must repeat the algorithm for various initial conditions and find the best solution out of all of the candidates. If controls are bounded, optimal controls are given by bang-bang controls using the Pontryagin minimum principle. Using a switching map we construct optimal solutions consisting of singular arcs. If controls are bounded, the analysis of our model also implies classical analysis done previously for this problem.

Dedicated to Jürgen Scheurle

1. Introduction

Control of quantum conservative (Hamiltonian) systems has been widely studied in the last few decades from both theoretical and interdisciplinary points of view, see for example [7], [10], [14], [15], [16]. Recently, there has been a growing interest in control of open (dissipative, non-Hamiltonian) quantum systems because of their applications to physics, chemistry and quantum computing. For example, there has been interest in the control of the rotation of a molecule in the gaseous phase by using laser fields in dissipative media [19], where the dissipation is due to molecular collisions. Other applications include control of spin dynamics by magnetic fields in nuclear magnetic resonance [12], as well as the construction of quantum computers [20].

The aim of this work is to study time-minimal optimal control of two-level quantum systems in a dissipative environment, where we assume that the dissipation is Markovian (the dynamics depend only on the present state and not its history) and time-independent. In this case the evolution of the density matrix of the system can be described by a quantum dynamical semi-group and the Lindblad master equation [1], [8], [18]. The boundary value problem determined by the time-minimum singular optimal control problem is studied numerically with unbounded and bounded controls respectively. A preliminary version of this work, focusing mainly on the energy-minimization problem with unbounded control, appears in the conference paper [11]. We extend the results of [11] to the time-minimum optimal control problem including bounded controls.

We study time-minimal controls with $\lVert u\rVert_{\infty}\leq 1$. In our feedback classification of trajectories, firstly we reproduce the results in [3], [4], since their system fits in those studied in this work for particular choices of parameters. We also complete the study done in [3] and [4], where controls are constant and unbounded.

The state space for a closed finite-dimensional quantum system is an $n$-dimensional projective Hilbert space, $P(\mathcal{H})$. Typically one drops the requirement of projectiveness, and instead work with unit vectors. To preserve the length of the vectors, time evolution is unitary $\displaystyle{U(t_{1},t_{2})|\psi(t_{1})\rangle=|\psi(t_{2})\rangle}$, and the evolution is described by the Schrödinger equation

[TABLE]

where $H$ is the Hermitian Hamiltonian. Here the ket-bracket describes the vector associated with an observable state.

The density operator, $\rho$, describes a probabilistic ensemble of states. It is given by a positive semi-definite Hermitian operator $\rho$ with $\text{Tr}(\rho)=1$ and $\text{Tr}(\rho^{2})\leq 1$. The purity of a density operator describes how close $\rho$ is to a single state. It is usually defined as $\displaystyle{P_{2}(\rho)={\text{Tr}(\rho^{2})}\in[1/n,1]}$, where the unique operator that has a purity of $1/n$ is $\frac{1}{n}I_{n\times n}$, called the completely mixed state. The dynamics for the purity operator is described by the von Neumann equation

[TABLE]

Notice that this dynamical system preserves the purity of $\rho$ (because the system is iso-spectral). A consequence is that if the quantum system is controlled by its Hamiltonian, there is no control over its purity, since one cannot directly alter the eigenvalues.

Open quantum systems are quite different. For such systems, dissipation occurs when we allow the system to interact with the environment. The full picture is an integro-differential equation called the Nakajima-Zwanzig (NZ) equation. To make the dissipation purely of differential form, one usually makes two assumptions: the dissipation is Markovian (i.e. the dissipation only depends upon the current state, not past history) and the dissipation is time-invariant.

Under these assumptions, the dynamics of the density operator is given by the Lindblad master equation [18]

[TABLE]

where the $L_{j}$ are called the Lindblad operators, $N$ denotes the number of Lindblad operators, and $\{\cdot,\cdot\}$ is the anti-commutator: $\{A,B\}=AB+BA$. We assume the Lindblad operators are traceless.

The structure of the paper is as follow: Section $2$ introduces the Lindblad equation and the reformulation of the equations as a first-order dynamical system on the unit (Bloch) ball. The time-minimal control problem studied in this work is formulated after defining the apogee and escape chimney. Section $3$ is devoted to the study of the control dynamics when there is only one Lindblad term. The case of constant controls is also analyzed. Section $4$ studies the time-minimum optimal control problem. The resulting boundary value problem is studied numerically in two scenarios: unbounded controls and bounded controls.

2. Statement of the problem

2.1. The Lindblad equation

An open quantum system is described by a density operator $\rho$, which is a trace-one positive semi-definite Hermitian operator on an $n$-dimensional complex Hilbert space $\mathcal{H}$. If the dissipation is Markovian and time independent, the density operator obeys the Lindblad equation (see [8] and [18] for details)

[TABLE]

where $N$ denotes the number of Lindblad operators, $[\cdot,\cdot]$ denotes commutator of matrices, $H$ is the Hermitian Hamiltonian, $\dagger$ represents the Hermitian transpose and $\left\{L_{j}\right\}$ are the Lindblad operators. The purity of the system is defined as $P_{2}(\rho)={\text{Tr}(\rho^{2})}$. We assume that controls appear in the Hamiltonian $H$, and not in the Lindblad operator $\displaystyle{\mathcal{L}_{D}:=\sum_{j=1}^{N}L_{j}\rho L_{j}^{\dagger}-\frac{1}{2}\{L_{j}^{\dagger}L_{j},\rho\}}$.

When $n=2$, the density operator can be identified with a vector in the Bloch ball (the unitary sphere with its interior) [21], [22]. In this special case, we can change from the operator point-of-view to that of the Bloch ball by considering $\rho=\frac{1}{2}\left(I+q_{1}\sigma_{1}+q_{2}\sigma_{2}+q_{3}\sigma_{3}\right)$ where $q\in S^{2}$, $I$ is the $2\times 2$ identity matrix, and $\sigma_{j}$ are the Pauli matrices.

We can reformulate the Lindblad equation (3) as a first-order dynamical system on the unit ball. Using the derivation given in [21] (see appendix A therein), equation (3) is equivalent to

[TABLE]

where the bar represents the complex conjugate of matrices. The vectors $\vec{l}_{j}$, $\vec{u}$ are the traceless parts of $L_{j}$ and $H$ respectively with $H=\displaystyle{h_{0}I+\sum_{k=1}^{3}\!u_{k}\sigma_{k}}$ where $\sigma_{k}$ are the Pauli matrices. Notice that the matrix $A$ is positive semi-definite. We will define $B:=A-\text{tr}(A)$.

The controlled Hamiltonian dynamics cannot achieve purity one (see [25]) and in general cannot affect the purity of the state or transfer the states between unitary orbits. To control purity, one must use the dissipative dynamics to move between orbits as in [21], [22] and [23]. If $q\neq 0$, equation (4) can be written in terms of the radial component (i.e., $r=\lVert q\rVert$).

In [21] it is shown that the purity $P_{2}(\rho)$ is equal to $\displaystyle{\frac{(1+r^{2})}{2}}$, where $r=\lVert q\rVert$. Knowing that $r^{2}=\langle q,q\rangle$, we find that

[TABLE]

where $\hat{q}$ is the unit vector associated with $q$, $\hat{q}=q/||q||$. Therefore,

[TABLE]

So, we can control the purity via the orientation of the corresponding unit vector. Hereafter we refer to $\displaystyle{f(q)}$ as the purity derivative.

2.2. The apogee and the escape chimney

The purity derivative (5) is independent of the controls $u_{j}$. It is illuminating to examine the regions in the Bloch ball where the purity derivative is positive. To do this, we examine the zeros of $f$.

Define two sets, $\mathcal{U}=\{q|f(q)\geq 0\}$ and the ellipsoid $\mathcal{M}=\{q|f(q)=0\}$. To find $\mathcal{M}$, we define a new function $f_{q}(r):=f(rq)$ where $q\in S^{2}$. Finding the roots of $f_{q}$ will let us solve for $\mathcal{M}$ in spherical coordinates. Note that

[TABLE]

so, the nonzero root is

[TABLE]

Notice that $g$ is always defined, since $B$ is negative-definite. Also note that the maximum of $g$ must be bounded by $1$, so the Bloch ball is invariant under (4).

We define $q_{apogee}$ to be the apogee of the ellipsoid $\mathcal{M}$. i.e.

[TABLE]

which can be found by maximizing (7) on $S^{2}$. Therefore, the apogee of the ellipsoid $\mathcal{M}$ will be the state with the maximal achievable purity. We call the interior of the ellipsoid $\mathcal{U}$ the escape chimney.

The control problem consists of finding a trajectory of the state variables solving (5), starting at the completely mixed state (i.e., $r=0$) and ending at the apogee. It is important to note, however, that the dynamics (5) has a singularity at the origin since

[TABLE]

Also note that the apogee cannot be reached in finite time since it is not possible to reach equilibrium points in finite time. Since the purity derivative is independent of controls, the purity obeys an autonomous first-order dynamical system that cannot reach its fixed point.

By considering the boundary conditions representing initial and final states

[TABLE]

with $\varepsilon>0$ and $\delta>0$ sufficiently small, the optimal control problem can be stated as follows:

Problem: Let $J:S^{2}\times U\rightarrow\mathbb{R}$ be a cost functional dependent on the state as well as the controls. The time-minimum optimal control problem consists of finding a control, $u^{*}:[0,t_{f}]\rightarrow U$ satisfying the dynamics (5) such that $q(0)=q_{0}$, $q(t_{f})=q_{f}$, and $\displaystyle{u^{*}=\arg\min\!\int_{0}^{t_{f}}dt.}$

3. Dynamics for Lindblad equations

In the special case when $N=1$ in (3), $\vec{b}$ becomes an eigenvector of $B$. This fact lets us simultaneously diagonalize $B$ and rotate $\vec{b}$ into the first coordinate. By additionally taking $u_{1}=u_{2}=0$, the third component of $q$ in equation (4) becomes uncontrolled and exponentially decays to zero. Dropping this coordinate, our system collapses to a two-dimensional underactuated bi-linear control system:

[TABLE]

where $q=(x,y)$, $\vec{b}=(b_{1},b_{2})$, $u=u_{3}$, $\alpha_{1},\alpha_{2}<0$ are the coefficients of the matrix $A$, and the maximum of (6) is bounded by $1$. We assume without loss of generality that $\alpha_{2}<\alpha_{1}$.

Remark 3.1**.**

System (10) reduces in special cases to those studied in [5] and [24]. Our system (10) is of the form $\dot{q}=F(q)+uG(q)$, $F(q)=(b_{1}+\alpha_{1}x,b_{2}+\alpha_{2}y)^{T}$ and $G(q)=(-y,x)^{T}$ while in [5] and [24], after the change of coordinates $x=y$, $y=z$, the drift term and control vector field are $F(q)=(-\gamma x,\gamma_{-}-\gamma_{+}y)$ and $G(q)=(-y,x)$ with parameters $\Gamma,\gamma_{+},\gamma_{-}$. Therefore, by setting $b_{1}=0,\alpha_{1}=-\Gamma,b_{2}=\gamma_{-},\alpha_{2}=-\gamma_{+}$ in (10) we recover the situation in [5] and [24]. Note also that the constraints in [5] and [24] become $-2\alpha_{1}\geq-\alpha_{2}\geq|b_{2}|$, which always holds in our case as long as $\alpha_{2}<\alpha_{1}<0$. $\diamond$

Note that in the absence of controls, the (asymptotically stable) fixed point for the drift field is $q^{*}=(q_{1}^{*},q_{2}^{*}):=-B^{-1}b$. To determine when $q^{*}=q_{f}$ (when the apogee is at $q^{*}$), one can use polar coordinates and study the critical points of $g(\theta)$. Denoting $\theta^{*}:=\hbox{atan2}\left(\frac{q_{2}^{*}}{q_{1}^{*}}\right)$, and using MATLAB’s symbolic toolbox:

[TABLE]

Hence, $\theta^{*}$ is a critical point if $\alpha_{1}=\alpha_{2}$, or $b_{1},b_{2}=0$. Looking at the sign of the second-derivative for $\alpha_{1}=\alpha_{2}$:

[TABLE]

and therefore $g$ achieves its local maximum at $\theta^{*}$.

If $b_{1}=0$,

[TABLE]

and therefore $g$ achieves its local maximum at $\theta^{*}$ if $2\alpha_{1}-\alpha_{2}<0$. If $b_{2}=0$, we similarly obtain the same conclusion under the condition $2\alpha_{2}-\alpha_{1}<0$. Therefore we conclude that when controls are turned off, we can achieve maximum purity if (1) $\alpha_{1}=\alpha_{2}$, (2) $b_{1}=0,2\alpha_{1}-\alpha_{2}<0$, or (3) $b_{2}=0,2\alpha_{2}-\alpha_{1}<0$.

While using a constant control is not the optimal method, it offers the most simplistic choice of control for our problem. In the absence of controls, the trajectory hits the ellipsoid $\mathcal{M}$ either to the left or to the right of $q_{f}$. In the presence of controls the system introduces a swirl to the dynamics. We would like to create a constant swirl such that $q_{f}$ the end-state of the system.

Assume we have constant controls $u$. Fixed points for the dynamics determined by (10) are

[TABLE]

Letting $\theta^{*}=\text{atan2}(q^{*}_{2},q^{*}_{1})$, we want to find $u$ such that $g^{\prime}(\theta^{*})=0$. According to MATLAB’s symbolic toolbox, this is equivalent to finding a root for the following cubic polynomial (which is guaranteed to have a real root):

[TABLE]

Observe that $u=0$ is a solution to (11) exactly when $b_{1}=0$, $b_{2}=0$, or $\alpha_{1}=\alpha_{2}$, which is consistent with what we expect.

We are interested in studying time-minimum control problems with bounded controls, so a natural question to ask is which conditions are necessary for the root of (11) to be bounded by 1?

We rewrite (11) as $p(u)=au^{3}+bu^{2}+cu+d$ (note that $a>0$). We need to first determine how many real, distinct roots there are. Let $\Delta:=4b^{2}-12ac$ be the discriminant of $p^{\prime}(u)$.

1. (1)

$\Delta\leq 0$: Equation (11) can only have one real root. So necessary conditions for the root to be bounded by 1 is for $\text{sgn}(p(-1)/p(1))=-1$, or $p(1)=0$. In this case, we require:

[TABLE]

or

[TABLE] 2. (2)

$\Delta>0$: Equation (11) can have 1, 2, or 3 real roots. If the above condition is satisfied, then we have an odd number of roots in $[-1,1]$, but it will not detect an even number of roots. Assume there are two roots, $u_{1},u_{2}\in[-1,1]$. Then, by the mean value theorem, there must exist $4u^{*}\in(u_{1},u_{2})$ such that $p^{\prime}(u^{*})=0$ (which can be found by the quadratic formula).

So the conditions for a root in $[-1,1]$ are

[TABLE]

or

[TABLE]

or

[TABLE]

where

[TABLE]

4. Time-minimal optimal control problem

In this section we study the time-minimum optimal control problem. The resulting boundary value problem is studied numerically. We study two different scenarios: unbounded and bounded controls. If controls are unbounded, simulations show that multiple local minimal solutions might exist. If controls are bounded, optimal controls are given by bang-bang controls which are obtained by studying the associated adjoint system to the Pontryagin minimum principle.

4.1. Time-minimal optimal control problem with unbounded controls

We want to find (unbounded) controls that steer (10) with end points (9) in the minimal amount of time possible. That is, we want to find a minimal solution to the functional

[TABLE]

where $x(0)=x_{0}$ and $x(t_{f})=x_{f}$. The derivative $dr/dt$ is found in $(\ref{eq:redial})$, and $dr/dx=x+yy^{\prime}$. So we wish to minimize a functional with integrand

[TABLE]

This Lagrangian is not hyperregular, so the Euler-Lagrange equations will fail to yield meaningful results [2]. To get around this problem, we implement the Rayleigh-Ritz numerical algorithm [13]. We assume $y(x)$ is a sum of linearly independent functions

[TABLE]

where $y^{0}(x_{0})=y_{0}$, $y^{0}(x_{f})=y_{f}$, and $y^{i}(x_{0})=y^{i}(x_{f})=0$. Specifically, we will take the following functions as a basis of polynomials for our approximation:

[TABLE]

with $M$ an arbitrary integer. To minimize the functional (13), the following $M$ equations must hold:

[TABLE]

where $\displaystyle{I_{c_{i}}=\frac{\partial I}{\partial c_{i}}}.$ This can be done by symbolically computing $\partial L/\partial c_{i}$ and numerically integrating, using, for instance, a $4^{th}$ order Runge-Kutta method. In order to find the optimal values to the $c_{i}$’s, we construct the function $\nu:\mathbb{R}^{M}\to\mathbb{R}$ with:

[TABLE]

4.1.1. Numerical results:

We will work an example with parameter values $b_{1}=1$, $b_{2}=2$, $a_{1}=-3$ and $a_{2}=-4$. Additionally, we will take $\varepsilon=\delta=10^{-3}$.

Solving for the apogee (8) in polar coordinates, we get $q_{apogee}=[0.4079,0.4493]^{T}$. Figure $1$ shows a simulation of the trajectory for the $7^{th}$ order curve.

Since all the computations are done independent of time, we report a plot of $u$ versus $x$ in Figure $1$.

An expected downside of the simulations is that multiple local minimal solutions might exist. To find the unique global minimal solution, we repeat the algorithm for various initial conditions. We then find the best solution out of all of the candidates. For this case, we ran the algorithm 25 times and the initial conditions were uniformly randomly chosen in the $l^{\infty}$-ball of radius 2.

4.2. Time-minimal optimal control problem with bounded controls

The goal now is to perform time-minimal controls with $\lVert u\rVert_{\infty}\leq 1$. In our feedback classification of trajectories, we firstly reproduce the results in [3], [4], as a special case of our model by performing the time-minimal local syntheses [6]. We additionally present an example of our model that is not included in [3] or [4].

To find the optimal control, we use the Pontryagin minimum principle ([17] Section $5.4$). Consider the two-dimensional underactuated bi-linear control system

[TABLE]

From the Pontryagin minimum principle, the Hamiltonian is

[TABLE]

Here we denote $q=(x,y)$, $f(q)=b+Bq=(b_{1}+\alpha_{1}x,b_{2}+\alpha_{2}y)^{T}$ and $g(q)=(-y,x)^{T}$. Following the analysis performed in [17], the optimal controls are given by bang-bang controls as

[TABLE]

To represent trajectories with switching controls we follow the method of time-minimal local syntheses used in [3] and [6], which avoids using the adjoint system $\displaystyle{\dot{p}=-p\left(B+u\frac{\partial g}{\partial q}\right)}$, $u=\pm 1$, provided by Pontryagin principle. The main idea of the construction provided in [3] (Section $4.5$) is to construct the optimal control as the concatenation of arcs $\sigma_{+}$, $\sigma_{-}$ corresponding to $u=1$, $u=-1$, and singular arcs. The switching function is given by $\Phi(t)=\langle p(t),g(q(t))\rangle$ and the switching takes place when $\Phi(t)=0$. Note that $g$ rotates $\pi$ degrees when we follow an arc $\sigma_{+}$ or $\sigma_{-}$.

Assume $t_{0}=0$, $t_{1}=t>0$ are two consecutive switching times where the controls are $u=\varepsilon$, $\varepsilon=\pm 1$. Then by definition of $\Phi$, we must have

[TABLE]

As in [3] (Section $4.5.5$), let $v(\cdot)$ be the solution of the variational equation $\displaystyle{\dot{v}=\left(B+u\frac{\partial g}{\partial q}\right)v}$, $u=\pm 1$, such that $v(t)=g(q(t))$, where we integrate backwards from time $t_{1}$ to $t_{0}$. By construction, $\langle p(0),v(0)\rangle=0$. Thus $p(0)$ is orthogonal to $g(q(0))$ and $v(0)$, and therefore $v(0)$ and $g(q(0))$ are collinear. In other words, we want

[TABLE]

As it was shown in [6] (Chapter $2$), defining $\theta(t)$ to be the angle between $g(q(0))$ and $v(0)$, the switching occurs when $\theta(t)=0$ modulo $\pi$ due to reflections along the vertical axis and one has

[TABLE]

In the analytic case, this can be solved by a power series.

To explicitly compute (20), we lift our system to the semi-direct product Lie group $G=GL_{2}(\mathbb{R})\ltimes\mathbb{R}^{2}$. Its corresponding semi-direct product Lie algebra, $\mathfrak{g}=\mathfrak{gl}_{2}(\mathbb{R})\ltimes\mathbb{R}^{2}$, has a Lie bracket defined by

[TABLE]

for $A,B\in\mathfrak{gl}_{2}(\mathbb{R})$ and $a,b\in\mathbb{R}^{2}$.

We want to compute $\text{ad}(F+\varepsilon G)$ where $F=(B,b)$ and $G=(C,0)$ with $C(q)=g(q)$. To compute $\text{ad}(F+\varepsilon G)$, we take the canonical basis for $\mathfrak{gl}_{2}(\mathbb{R})\ltimes\mathbb{R}^{2}$ denoted by $\{e_{i}\}_{i=1}^{6}$ where

[TABLE]

$e_{5}=(1,0)^{T}$ and $e_{6}=(0,1)^{T}$.

Let $Z=F+\varepsilon G$. Computing the brackets, we get:

[TABLE]

Compiling (21), we get the adjoint representation

[TABLE]

We want to find $e^{-t\,\text{ad}(Z)}(-e_{2}+e_{3})$. Because this involes exponentiating a $6$ by $6$ matrix, we used MATLAB. To find the collinearity and the switching times, we want to solve for $\text{det}(G(q_{0}),e^{-t\text{ad}(Z)}(-e_{2}+e_{3})(q_{0}))=0$. Using numerical methods, we find $t$ such that

[TABLE]

and switch $\varepsilon$ to $-\varepsilon$ at that point.

In Figure 2 we show a simulation for a specific choice of parameters to mimic the results given in [3] and [4]. In the figure, the blue curves represent a control $u=1$ and the red curve if $u=-1$. The green curve and light blue curve are switching signs of the controls (as predicted by the above analysis). Finally, the yellow curve and purple curve corresponds to the case when one never switched signs and one continues with the original controls.

In the general case where $b\neq 0$ (see Figure 3), the switching times are no longer symmetric. The switching times for this example are

5. Final discussion

We studied time-minimum global optimal control problems for dissipative open quantum systems where the dynamics is described by the Lindblad equation and controls are both unbounded and bounded. We have transformed such a control system into a bi-linear singular control system in the Bloch ball and have come up with the construction of a numerical algorithm to design optimal paths to achieve a desired point given initial states close to the origin of the Bloch ball in both optimal control problems.

We would like to extend our results to higher-order dimensional systems, as well as study energy-minimum global optimal control problems for dissipative open quantum systems where the controls are bounded. It would also be interesting to determine the best basis of functions for the Rayleigh-Ritz methods as well as the best order of solutions to use.

Acknowledgments: This research was supported in part by grants from the National Science Foundation and Air Force Office of Scientific Research.