| Izvestiya: Mathematics |
|
|
|
|
|
This article is cited in 3 scientific papers (total in 3 papers) On optimization of coherent and incoherent controls for two-level quantum systems O. V. Morzhina, A. N. Pechenab a Steklov Mathematical Institute of Russian Academy of Sciences, Department of Mathematical Methods for Quantum Technologies, Moscow, Russia b National University of Science and Technology «MISIS», Moscow, Russia
Russian version:
Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 2023, Volume 87, Issue 5, DOI: https://doi.org/10.4213/im9372 
Bibliographic databases:
    Dedicated to the centenary of Academician V. S. Vladimirov § 1. IntroductionOptimal control of quantum systems (atoms, molecules, etc.) is an important research direction at the intersection of mathematics, physics, chemistry, and technology, which is crucial for development of various quantum technologies [1]–[12]. Quantum control aims to obtain desired properties of the controlled quantum system using external variable influence. For example, in the theory of quantum computing, the basic objects are qubits and quantum gates. Mathematical modeling of physical processes for generating quantum gates considers, for example, the Schrödinger equation with coherent control in the Hamiltonian and objective functional that characterizes how well a certain quantum gate is generated. In optimal quantum control, the Pontryagin maximum principle (PMP) [13] (about the PMP adaptation for various optimal control problems for quantum systems, see, for example, Chap. IV in [3] and the survey [14]), methods of the geometric control theory [15], iterative optimization methods, including gradient-based methods (for quantum control, see, for example, [16]), the Krotov method (see the survey [17] for references) are widely used. Control problems for two-level closed and open quantum systems were considered, for example in [18]–[24]. Various master, kinetic and transport equations are used for modelling dynamics of open classical and quantum systems. V. S. Vladimirov1 contributed to development of numerical schemes for solving kinetic and diffusion equations, including those for neutron transport [25] and [26]. An important problem is to study the dynamic landscape of the quantum control problem, that is geometric properties of the corresponding objective functional such as existence or absence of traps of various order, see [27]–[40]. On this topic, a rigorous proof of absence of traps was obtained for two-level coherently controlled closed quantum systems [29]–[31], which is the first and the only existing mathematically rigorous proof of trap-free behaviour for finite-level quantum systems, which is free of any assumptions. Absence of traps in the kinematic control landscape for open quantum systems was proved for the two-level case in [41], where for the first time parametrization of Kraus map open system evolution by points in the complex Stiefel manifolds was proposed and theory of optimization over complex Stiefel manifolds, including analytical computation of the gradient and Hessian of the quantum control objective, was developed for quantum control. The detailed theory for the multilevel case was developed in [42]. Control landscapes for open-loop and closed-loop control were analyzed in a unified framework [38]. A unified analysis of classical and quantum kinematic control landscapes was performed [39]. These results were used to explain unexpected, due to high dimensionality, efficiency of optimization in chemistry [8] and biology [10]. In [32], for the problem of generating the phase shift quantum gate, analytical study of the Hessian for the objective functional was performed at zero control, which is a stationary point of the objective functional. Also numerical exploration of the control landscape was performed using GRAPE (gradient ascent pulse engineering, see [43]) and stochastic methods of global optimization. The conditions of optimality of zero control and ability to exactly generate the phase shift quantum gate, including in minimal time, were found. To continue the analysis performed in [32], in § 3 of the present work, we numerically show using GRAPE, which is the most suitable to reveal the structure of a quantum control landscape, and also differential evolution [44] and dual annealing [45], [46] (using the stochastic methods implementations available in SciPy) that zero coherent control which satisfies the PMP for arbitrary magnitude constraints is not optimal; it thereby gives an example of subtle point for studying quantum control problems. In [47], [48], incoherent control method for quantum systems was proposed using controlled drive in the dissipator of the Gorini–Kossakovsky–Sudarshan–Lindblad (GKSL) master equation. Such approach exploits incoherent action of the environment surrounding the system to efficiently manipulate the system dynamics via diffusion and decoherence. Incoherent control in this approach can be supplied by standard coherent control. In [47], two general classes of master equations were studied for quantum control, which correspond to weak coupling limit (WCL) and collisional based low density limit (LDL) [49]–[53] in the theory of open quantum systems. Higher order corrections to the WCL [54] can be included when necessary. Diffusion in collision-free models of quantum particles is also studied, for example, in a quantum Poincaré model realizing behaviuor of ideal gas of noninteracting quantum Bolztman particles [55]. The papers [47], [48] present the general theory for $N$-level quantum systems and examples for the cases with $N=2, 4$. The work [48] contains fundamental results about approximate controllability of $N$-level quantum systems, whose dynamics is determined by coherent and incoherent controls. In [56], an analytical description of reachable sets for two-level systems driven by coherent and incoherent controls was obtained using geometric control theory and surprisingly, unreachable states in the Bloch ball were discovered. For the case of coherent and incoherent control of a two-level open quantum system, various quantum control problems were studied (generation of a target density matrix, including in minimal time, maximizing Hilbert–Schmidt scalar product between final and target density matrices, maximizing the Uhlmann–Jozsa fidelity, estimating reachable and controllability sets, etc.) for various classes of controls and using various optimality conditions and optimization methods [56]–[60]. Quantum coherent and incoherent control of other systems has been studied, such as two-qubit systems (for example, [61]) and quantum harmonic oscillator [62]. A model of incoherent control for the optical signals using the physics of a quantum heat engine was proposed [63]. A particular result of [48] is a two-stage method for approximate generation of a given target density matrix which was proposed and developed for generic $N$-level quantum systems and provided an explicit analytical form of incoherent control on the first stage. At the first stage, coherent control is zero and constant in time incoherent control is applied, which is analytically computed using eigenvalues of the target density matrix; this incoherent control drives the system state to a density matrix which approximately coincides with the intermediate target density matrix $\widetilde{\rho}_{\mathrm{target}}$ whose eigenvalues are equal to eigenvalues of the target density matrix $\rho_{\mathrm{target}}$. At the second stage, incoherent control is zero and optimized coherent control is applied to obtain the final density matrix $\rho(T)$ which with some sufficient accuracy coincides with $\rho_{\mathrm{target}}$. This method is capable of approximately creating arbitrary density matrices of generic $N$-level quantum systems, but it does not consider time optimality, which determines the importance of the problem of minimizing time necessary for the first stage which is typically most time consuming. In this regard, in § 4 of the present paper, we propose and compare with the original case the following modification of the first stage of the two-stage method for two-level systems: instead of using constant incoherent control, optimization in the class of piecewise constant incoherent controls is performed. For this modification, we give exact analytical formulas for the quantum system state at the end of the first stage, the objective functions, and their gradients depending on the parameters defining piecewise constant controls. Then we use these formulas for implementation of the two-step gradient projection method (GPM-2) – a projection version [64] of the heavy ball method [65], [66]. § 2. Formulation of the optimal control problemsThe paper [32] considers the problem of generating the phase shift quantum gate for a two-level quantum system whose dynamics is determined by the Schrödinger equation for evolution operator $U(t)$ and with objective functional $J_W$,
$$
\begin{equation}
\dfrac{dU(t)}{dt} = -i (H_0+v(t) V) U(t), \qquad U(0)=\mathbb{I}_2, \quad t \in [0, T],
\end{equation}
\tag{2.1}
$$
$$
\begin{equation}
J_W(v) = \frac{1}{4} |{\operatorname{Tr}(W^{\unicode{8224}} U(T))}|^2 \to\max, \qquad W=e^{i \varphi_W \sigma_z},
\end{equation}
\tag{2.2}
$$
where $U(t)$ and $W$ are $2 \times 2$ unitary matrices; $H_0$, $V$ are $2 \times 2$ Hermitian matrices; $v$ is coherent control being, in general, a measurable function; $i$ is the imaginary unit. We consider the system of units where the Planck constant is $\hbar=1$. The goal is to find a unitary operator (matrix) $U(T)$ for some $T$ such that $U(T)$ coincides with $W$ or is as close as possible to $W$. The parameter $\varphi_W$, which appears in $W$, and the final time $T$ allow one to parameterize the optimal control problem. The work [32] considers various domains on the coordinate plane of $\varphi_W$ and $T$. In the present paper, we consider the set $D := \{(\varphi_W, T)\colon \varphi_W \in (0,\pi/2), \, T \in (0,\pi/2] \}$. The articles [47], [48] contain the general theory of controlling open $N$-level quantum systems driven by coherent and incoherent controls, and examples for $N=2, 4$. In [57], the case $N=2$ is considered with piecewise continuous coherent $v$ and incoherent $n$ controls. The master equation for the density matrix $\rho(t)$, which is $2 \times 2$ Hermitian positive semidefinite matrix with unit trace ($\rho(t)= \rho^{\unicode{8224}}(t) \geqslant 0$, $\operatorname{Tr}\rho(t)=1$), has the form
$$
\begin{equation}
\frac{d \rho(t)}{dt}=-i[ H_0+v(t) V, \rho(t)]+\mathcal{L}_{n(t)}(\rho(t)),\qquad \rho(t) \in \mathbb{C}^{2 \times 2}, \qquad \rho(0)=\rho_0.
\end{equation}
\tag{2.3}
$$
The initial density matrix $\rho_0$ is fixed (in contrast, when studying controllability sets [60], $\rho_0$ was not fixed). The dissipation superoperator $\mathcal{L}_{n(t)}$ acts on the density matrix as
$$
\begin{equation}
\begin{aligned} \, \mathcal{L}_{n(t)} (\rho(t)) &= \gamma(n(t)+1) \biggl( \sigma^- \rho(t) \sigma^+- \frac{1}{2}\{\sigma^+ \sigma^-, \rho(t)\} \biggr) \nonumber \\ &\qquad+\gamma n(t) \biggl( \sigma^+ \rho(t) \sigma^--\frac{1}{2} \{ \sigma^- \sigma^+, \rho(t)\} \biggr), \qquad \gamma>0, \end{aligned}
\end{equation}
\tag{2.4}
$$
and describes the controlled interactions between the system and its environment (reservoir), where control $n$ is scalar. The notations $[A,B]= AB-BA$ and $\{A,B\}=AB+ BA$ mean commutator and anticommutator of operators $A,B$; the matrices
$$
\begin{equation*}
\sigma^+=\begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix},\qquad \sigma^-=\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}
\end{equation*}
\notag
$$
determine the transitions between two energy levels, control $n$ is used to regulate this interaction. The following constraint is obligatory by the physical meaning of incoherent control: In the closed (2.1) and open (2.3) systems, the following forms of free Hamiltonian $H_0$ and interaction Hamiltonian $V$ are considered:
$$
\begin{equation*}
\sigma_x=\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix},\qquad \sigma_y=\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix},\qquad \sigma_z=\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}.
\end{equation*}
\notag
$$
We set $\alpha_x=\cos\theta$ and $\alpha_y=\sin\theta$. As noted in Lemma 5 of [30], one can take $\theta=0$ in $V$, taking into account the invariance of $J_W$ with respect to $\theta$. Then $\alpha_x=1$, $\alpha_y=0$, and $V=\sigma_x$. Thus, in the closed and open systems under consideration, the free Hamiltonians differ from each other additively by physically unrelevant $\mathbb{I}_2$ and by the scale factor (which is equivalent to time rescaling), while the interaction Hamiltonians, taking into account $\theta=0$, are almost the same, they differ by the factor $\mu$ that is equivalent to rescaling the magnitudes of coherent control. Hence, up to rescaling time and magnitudes of coherent control, they are the basically same. If we set $\gamma= 0$, then system (2.3) becomes a closed system whose dynamics is described by the von Neumann equation (as in [3], Chap. IV). Moreover, taking into account the relation $\rho(t)=U(t) \rho_0 U^{\unicode{8224}}(t)$, system (2.3) for $\gamma=0$ can be described by (2.1) with the corresponding Hamiltonian. The present paper considers the following classes of controls. In § 3, for the closed system (2.1), we consider the class of piecewise continuous controls when discussing the stationary point of the objective functional $J_W$ and the PMP, and then we consider the following class of piecewise constant controls for reduction to a finite-dimensional optimization problem
$$
\begin{equation}
v(t)=\sum_{k=1}^N a_k \chi_{[t_k, t_{k+1})}(t), \qquad v(T)=v(T-), \quad \Delta t= \frac{T}{N},
\end{equation}
\tag{2.6}
$$
where $a_k \in \mathbb{R}$ is parameter representing control $v$ during a partial interval $[t_k, t_{k+1})$ (“shelf”); $\chi_{[t_k, t_{k+1})}$ is the characteristic function of this interval; moreover, as an option, the constraint $|a_k| \leqslant \nu$ with some $\nu>0$ can be used. In § 4, for the open system (2.3), when we formulate the modification of the first stage of the two-stage method, the coherent control is $v=0$ and incoherent control is piecewise constant
$$
\begin{equation}
n(t)=\sum_{k=1}^N a_k \chi_{[t_k,t_{k+1})}(t), \qquad a_k \geqslant 0, \quad \Delta t= \frac{\widehat{t}}{N},
\end{equation}
\tag{2.7}
$$
where $\widehat{t}>0$ is the end time of the first stage. Here, constraint (2.5) is taken into account, and there may be the additional constraint $a_k \leqslant n_{\max}$, where $n_{\max}>0$ is some given value. In § 4, to implement the second stage with zero incoherent control, we consider (by analogy with [48]) the coherent control of the form
$$
\begin{equation}
v(t)=A \cos(\Omega t), \qquad t \in [\widehat{t}, T],
\end{equation}
\tag{2.8}
$$
where $A$ is amplitude and $\Omega$ is frequency. Moreover, in addition to the steering problem, we consider the condition $v(\widehat{t}\,)=v(T)=0$ for coherent control $v$, that is, smooth switching on and off, correspondingly, at the moments $\widehat{t}$ and $T$. In this regard, we consider the following class of controls:
$$
\begin{equation}
v(t)=A \sin\frac{\pi d (t-\widehat{t}\,)}{T-\widehat{t}}, \qquad d \in \{ 1, \dots, \overline{d}\}, \quad T>\widehat{t}, \quad t \in [\widehat{t}, T], \quad \overline{d} \in \mathbb{N}.
\end{equation}
\tag{2.9}
$$
In § 4, we present a version of the two-stage method. At the first stage, we consider two variants of the problem of minimizing the squared Hilbert–Schmidt distance between $\rho(\widehat{t}\,)$ and $\widetilde{\rho}_{\mathrm{target}}$:
§ 3. Nonoptimality of zero control satisfying the first order necessary optimality conditions for the quantum gate generation3.1. Realification of $U(t)$, satisfying the first order necessary optimality conditions at the control $v=0$Various equivalent representations can be used here. Using the representation
$$
\begin{equation*}
U(t)=\begin{pmatrix} x_1(t)+i x_2(t) & x_3(t)+i x_4(t) \\ -x_3(t)+i x_4(t) & x_1(t)-i x_2(t) \end{pmatrix},\qquad x_j(t) \in \mathbb{R},\quad j=1,\dots,4
\end{equation*}
\notag
$$
(by analogy with [67]), system (2.1) with $H_0=\sigma_z$ and $V= \sigma_x$ is rewritten as the bilinear system
$$
\begin{equation}
\dot x(t)=(A+B v(t))x(t), \qquad x(0)=(1,0,0,0)^\top.
\end{equation}
\tag{3.1}
$$
Here, $A=\left(\begin{smallmatrix} A_2 & 0_2 \\ 0_2 & A_2 \end{smallmatrix}\right)$, $B=\left(\begin{smallmatrix} 0_2 & B_2 \\ -B_2 & 0_2 \end{smallmatrix}\right)$, where $A_2=\left(\begin{smallmatrix} 0 & 1\\ -1 & 0 \end{smallmatrix}\right)$, $B_2=\left(\begin{smallmatrix} 0 & 1 \\ 1 & 0 \end{smallmatrix}\right)$; $0_2$ is the $2 \times 2$ zero matrix; “$\top$” means transpose. The objective functional $J_W$ (2.2) to be maximized is rewritten as
$$
\begin{equation}
J_W(v)=\langle x(T), L x(T) \rangle=\bigl( x_1(T) \cos\varphi_W+x_2(T) \sin\varphi_W \bigr)^2,
\end{equation}
\tag{3.2}
$$
where the final state $x(T)$ of system (3.1) is computed for a certain control $v$; symmetric matrix $L=\left(\begin{smallmatrix} L_2 & 0_2 \\ 0_2 & 0_2 \end{smallmatrix}\right)$, where $L_2=\left(\begin{smallmatrix} \cos^2\varphi_W & \cos\varphi_W \sin\varphi_W \\ \cos\varphi_W \sin\varphi_W & \sin^2 \varphi_W \end{smallmatrix}\right)$. For this reformulation of the problem of maximizing the objective functional $J_W$ in terms of (3.2) and (3.1), we consider the first-order necessary optimality conditions known in the theory of optimal control [13], [68], [69]. The Pontryagin function is $H(p,x,v)=\langle p, (A+Bv)x\rangle$, where the variables $x,p \in \mathbb{R}^4$ and $v \in \mathbb{R}$, and the conjugate system
$$
\begin{equation}
\dot{p}(t)=-\bigl(A^\top+B^\top v(t)\bigr) p(t), \qquad p(T)=2Lx(T)
\end{equation}
\tag{3.3}
$$
for some admissible process $(x, v)$. For any pair $(\varphi_W, T)$, the control $v= v^0= 0$ satisfies the stationarity condition, which is a first-order necessary optimality condition for piecewise continuous controls without constraints on their magnitudes:
$$
\begin{equation}
\frac{\partial H(p,x,v)}{\partial v}\bigg|_{p=p^0(t), \, x=x^0(t)}=\langle p^0(t), B x^0(t)\rangle \equiv 0, \qquad t \in [0,T],
\end{equation}
\tag{3.4}
$$
where $x^0$ is the solution of system (3.1) for $v=v^0$, and $p^0$ is the solution of system (3.3) for the process $(x, v)=(x^0,v^0)$. If the pointwise constraint $v(t) \in Q := [-\nu, \nu]$ $\forall t \in [0, T]$ with some arbitrarily given $\nu>0$ is imposed, then we note that the pointwise maximum condition for the Pontryagin function in the PMP holds for any $\nu>0$ for the control $v^0=0$ in the class of piecewise continuous coherent controls (to the left of the “$\equiv$” sign, the letter “$v$” means variable $v \in \mathbb{R}$):
$$
\begin{equation}
\max_{v \in [-\nu, \nu]} H\bigl(p^0(t), x^0(t), v\bigr) \equiv H\bigl(p^0(t), x^0(t), v^0(t)\bigr), \qquad t \in [0,T].
\end{equation}
\tag{3.5}
$$
3.2. Results of the finite-dimensional optimizationIn connection with the conditions (3.4) and (3.5), we numerically study whether the control $v^0=0$ is a global solution to the problem of maximizing the objective functional $J_W$ for a given $T$ in the class of piecewise constant controls of the form (2.6) for various pairs $(\varphi_W^j, T^i) \in D$, where the indices $j,i$ specify points in the grid $\widehat{D} := \{\varphi_W^j=(\pi/20) j, \, j= 1,\dots,9, \ T_i=(\pi/20) i, \, i=1,\dots,10\} \subset D$ consisting of $90$ nodes (see Figure 1(a)); all the nodes on the line $T=\pi/2-\varphi_W$ are marked with another marker. For piecewise constant control of the form (2.6), we denote $\mathbf{a}=(a_1,\dots,a_N)\in\mathbb R^N$.  Figure 1.Visualization of the numerical optimization results obtained for $(\varphi_W, T)\in\mathcal{D}$: (a) grid consisting of $90$ nodes; (b) (results by the GRAPE-type method) interpolation surface obtained based on the data from Table 1 showing the differences between the numerically maximized $J_W$ and $J_W(v=0)$ for each node of the grid $\widehat{D}$; (c) (results by stochastic methods) the table-defined function $\{\widehat{J}_W^{\,l}\}_{l=1}^{90}$ shown by markers, representing the computed maximal values of the objective functional $J_W$, for each $l$th node of the grid on $\mathcal{D}$, and the interpolation surface corresponding to this table-defined function; (d) (results by stochastic methods) table-defined function $\{\Delta_l\}_{l=1}^{90}= \{\widehat{J}_W^{\,l}-J^l_W(v^0) \}_{l=1}^ {90}$, where $J^l_W(v^0)$ is the value of the objective functional for the control $v^0=0$, and the interpolation surface corresponding to this table-defined function Table 1.For the quantum gate generation problem: computed using analytical expression for gradient differences between numerically maximized $J_W$ and $J_W(v=0)$, that is, $J_W^{\max}-J_W(v=0)$, for the grid $\widehat{D}$
For each node $(\varphi_W^j,T^i) \in \widehat{D}$, the number $N$ of uniform partitions in $t$ is determined when defining piecewise constant controls of the form (2.6) so that $[0, T^i]$ is divided into $N$ equal parts. As in [32], we use the following two approaches: 1) GRAPE-type gradient search algorithm and exploiting exact analytically computed expression for gradient of the objective for piecewise constant controls which does not need solving differential evolution equation; 2) dual annealing and differential evolution methods together. For implementing GRAPE, we take $N=4+i$ and apply built in MATLAB function fminunc for unconstrainded optimization via quasi-Newton method using exact formula for gradient of the objective ($-J_W[\mathbf{a}]$). In numerical simulations, for each node, we run algorithm starting with $10$ random initial controls, where every control is produced by choosing $a_k$ randomly and uniformly distributed in the interval $[-A;A]$ ($A=1$ is chosen), and then select the best result among $10$. The explicit expression for gradient of the objective used for optimization is [32]
$$
\begin{equation*}
\nabla J_W[\mathbf{a}]=\frac{1}{2}\operatorname{Im} [\operatorname{Tr} (Y^\unicode{8224})\operatorname{Tr}(YV_{k})].
\end{equation*}
\notag
$$
Here,
$$
\begin{equation*}
\begin{gathered} \, Y= W^\unicode{8224} U_T^\mathbf{a},\qquad U^\mathbf{a}_T=U_N\cdots U_k \cdots U_1, \qquad V_{k}=U^\unicode{8224}_k\cdots U^\unicode{8224}_1 V U_1\cdots U_k, \\ U_k=\cos\alpha_k-i\delta t(\sigma_z+a_k\sigma_x)\frac{\sin\alpha_k}{\alpha_k},\qquad \alpha_k=\delta t\sqrt{1+a_k^2},\quad \delta t=\frac{T}{N}. \end{gathered}
\end{equation*}
\notag
$$
For each run, maximum number of iterations allowed is $10^{6}$, maximum number of function evaluations allowed is $10^{6}$, and the algorithm stops if $\|\nabla J_W\|_\infty\equiv \max_i|(\nabla J_W)_i|< 10^{-8}$. The results of optimization with the same stopping criteria are shown in Table 1 and Figure 1(b). Gradient-based search algorithms are among the most suitable methods for the analysis of quantum control landscape, especially in the case they experience some difficulties that may indicate some fine details of the underlying control landscape. The second approach uses two runs of the differential evolution method and two runs of the dual annealing method to minimize the objective function $f(\mathbf{a}) :=-J_W(v)$, where $\mathbf{a}=(a_k)_{k=1}^N$, $|a_k| \leqslant 50$, with comparison of their results; system (3.1) is solved numerically here. The minimization problem is considered here instead of the maximization problem, because the SciPy implementations of these methods are designed for minimization. In Figure 1(c), via the markers, the table-defined function $\{\widehat{J}_W^{\,l}\}_{l=1}^{90}$, which represents the computed optimal values of the objective functional $J_W(v)$ to be maximized for each $l$th node grid, and also the interpolation surface corresponding to this table-defined function are shown. Figure 1(d) shows the table-defined function $\{\Delta_l\}_{l=1}^{90}= \{\widehat{J}_W^{\,l}-J ^l_W(v^0) \}_{l=1}^{90}$, where $J^l_W(v^0)$ is the value of the objective functional on the control $v^0=0$, and also the interpolation surface corresponding to this table-defined function. As is known (see (6) in [32]), $J_W(v^0)=\cos^2(\varphi_W+T)$. According to the computations, we have $\min_{1 \leqslant l \leqslant 90} \{ \widehat{J}_W^{\,l} \}_{l=1}^{90} \approx 0.036$, $\max_{1 \leqslant l \leqslant 90} \{ \widehat{J}_W^{\,l} \}_{l=1}^{90} \approx 1$, $(1/90) \sum_{l=1}^{90} \widehat{J}_W^{\,l} \approx 0.843$ (see Figure 1(c)); we also have $\min_{1 \leqslant l \leqslant 90} \{ \Delta_l \}_{l=1}^{90} \approx 0.024$, $\max_{1 \leqslant l \leqslant 90} \{ \Delta_l \}_{l=1}^{90} \approx 1$, $(1/90) \sum_{l=1}^{90} \Delta_l \approx 0.565$ (see Figure 1(d)). Thus, it is shown that after reduction to the class of piecewise constant controls (2.6), the finite-dimensional optimization (GRAPE and stochastic methods) gives such piecewise constant controls for different pairs $(\varphi_W^j, T^i) \in D$ that the values of the objective functional $J_W(v)$ are substantially larger than for the control $v=v^0=0$. Since the class of piecewise constant controls of the form (2.6) is nested in the class of piecewise continuous controls for the same $T$, then we have the situation, when such controls are found in the class of piecewise continuous controls that are significantly better in terms of the values of $J_W(v)$ than for the control $v=v^0= 0$ satisfying the stationarity condition. In other words, the control $v=v^0=0$ is not globally optimal here. In addition, numerical results show that, for some part of the nodes of the grid introduced above, we have $J_W$ equal to $1$, that is, there is full generation of the quantum gate here. Note that in [30] and [32], Theorem 3, for a general class of measurable controls, it was shown that the control $v=v^0=0$ is a saddle point of the objective functional $J_W(v)$ for $(\varphi_W, T) \in D \setminus \{(\varphi_W, T)\colon \varphi_W+T \neq \pi/2\}$ (the Hessian of the functional $J_W$ has negative and positive eigenvalues). For a comparison, see [70], where, via an example of direct and inverse quantum Fourier transforms gates for qutrit with unitary dynamics (that is, for a 3-level closed quantum system), a relationship between the phase factor of a quantum gate and the arrangement of energy levels of the effective Hamiltonian was shown. § 4. Modification of the first stage of the two-stage method for approximate steering4.1. Reduction to real states ($x(t)$), and objective functions in terms of intermediate final state $x(\widehat{t}\,)$Following the two-level system example given in [48], the paper [57] uses the Bloch parameterization
$$
\begin{equation*}
\rho=\frac{1}{2} \biggl( \mathbb{I}_2+\sum_{j \in \{x, y, z\}} x_j \sigma_j \biggr) = \frac{1}{2} \begin{pmatrix} 1+x_3 & x_1-i x_2 \\ x_1+i x_2 & 1-x_3 \end{pmatrix}.
\end{equation*}
\notag
$$
Here, $x=(x_1, x_2, x_3) \in \mathcal{B} := \{ x \in \mathbb{R}^3 \colon \| x \| \leqslant 1\}$ is a Bloch vector representing some point in the Bloch ball; $\sigma_x$, $\sigma_y$, and $\sigma_z$ are the Pauli matrices; $x_j=\operatorname{Tr}(\rho \sigma_j)$, $j \in \{x,y,z\}$. The initial point $x_0$ has the coordinates $x_{0,j}=\operatorname{Tr}(\rho_0 \sigma_j)$, $j \in \{x,y,z\}$. Due to the bijection between density matrices and Bloch vectors, the Bloch parametrization allows one to carry out the study in terms of the equivalent dynamical control system whose state $x(t)$ at time $t$ is the Bloch vector. Various equivalent parametrizations of quantum systems can be used [71]. The following system is equivalent to (2.3):
$$
\begin{equation}
\frac{dx(t)}{dt} =\bigl(A+B^v v(t)+B^n n(t) \bigr)x(t)+d,\qquad x(0)=x_0 \in \mathcal{B};
\end{equation}
\tag{4.1}
$$
here
$$
\begin{equation*}
\begin{gathered} \, A=\begin{pmatrix} -\dfrac{\gamma}{2} & \omega & 0 \\ -\omega & -\dfrac{\gamma}{2} & 0 \\ 0 & 0 & -\gamma \end{pmatrix},\qquad B^v=\begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -2 \mu \\ 0 & 2\mu & 0 \end{pmatrix}, \\ B^n=\begin{pmatrix} -\gamma & 0 & 0 \\ 0 & -\gamma & 0 \\ 0 & 0 & -2\gamma \end{pmatrix},\qquad d=\begin{pmatrix} 0 \\ 0 \\ \gamma \end{pmatrix}. \end{gathered}
\end{equation*}
\notag
$$
In terms of the Bloch parametrization and taking into account the parameterization (2.7) of incoherent control, the following problems of finite-dimensional constrained optimization are considered in accordance with (2.10) and (2.11):
Here, $\widetilde{\rho}_{{\mathrm{target}},j}= \operatorname{Tr}(\widetilde{\rho}_{\mathrm{target}} \sigma_j)$, $j \in \{x, y, z\}$. 4.2. Exact analytical formulas for the system intermediate final state, objective functions and their gradientsBelow we show that the evolution equation (4.1) can be solved exactly for any initial state from the Bloch ball, any piecewise constant control $n$ under zero control $v$. This exact solution is used to represent the state $x(\widehat{t}\,)$ as some function of the parameters $\{a_k\}$, which determine the control $n$, and the moment $\widehat{t}$, if it is not fixed. In addition, the exact analytical formulas are obtained that explicitly show how the objective functions considered in the problems (4.2) and (4.3) depend on their arguments. Theorem 1. For a non-fixed or fixed final time $\widehat{t}$ and a vector $\mathbf{a}$ that define an incoherent control of the form (2.7) at the first stage of the method, the corresponding final state of system (4.1) is $x(\widehat{t}, \mathbf{a}) := \bigl(x_1(\widehat{t}, \mathbf{a}), x_2(\widehat{t}, \mathbf{a}), x_3(\widehat{t}, \mathbf{a}) \bigr)$, where The theorem will be proved in the next subsection. Due to this theorem we do not need to numerically solve system (4.1). Corollary 1. If $\mathbf{a}=(p, \dots, p)$ with some $p \geqslant 0$, then $x_i$ in Theorem 1 are as follows: Formula (4.10) is obtained from (4.7) using the factorization formula
$$
\begin{equation*}
1-y^N=(1-y) \sum_{s=1}^N y^{N-s}=(1-y) \biggl(1+\sum_{s=1}^{N-1} y^{N-s}\biggr).
\end{equation*}
\notag
$$
Moreover, § 5 of [61] considers the case when the control $v=0$ and incoherent control $n$ is constant over the entire time interval. Taking $t=\widehat{t}$ in (5.1) in [60], we also get (4.8)–(4.10). Theorem 2. For the functions $x_j(\widehat{t}, \mathbf{a})$, $j=1,2,3$, which are defined in (4.4), (4.5), and (4.7), their first order partial derivatives are as follows: 4.2.1. Proofs of Theorems 1, 2 Proof of Theorem 1. Consider equation (4.1) sequentially on partial intervals of length $\Delta t= \widehat{t}/N$ with initial conditions $x(t_1=0)=x_0$, $x(t_k)=x(t_k-0)$, $k=2, \dots, N$. For each $k$th initial condition, the corresponding Cauchy problem for the system of ordinary differential equations has the exact solution
$$
\begin{equation*}
\begin{aligned} \, x_1(t) &= \exp\biggl(-\frac{\gamma}{2}(1+2a_k)(t-t_k) \biggr) \bigl(x_1(t_k)\cos(\omega(t-t_k))+x_2(t_k)\sin(\omega(t-t_k)) \bigr), \\ x_2(t) &= \exp\biggl(-\frac{\gamma}{2}(1+2a_k)(t-t_k) \biggr) \bigl(x_2(t_k)\cos(\omega(t-t_k))-x_1(t_k)\sin(\omega(t-t_k)) \bigr), \\ x_3(t) &= \exp\bigl(-\gamma (1+2a_k)(t-t_k) \bigr) x_3(t_k) + \frac{1-\exp(-\gamma(1+2a_k)(t-t_k))}{1+2a_k}. \end{aligned}
\end{equation*}
\notag
$$
Thus, for $k=1, \dots, N$, we have
$$
\begin{equation}
x_1(t_{k+1}) = E^{\mathrm{I}}(\widehat{t}, a_k) \bigl(x_1(t_k)\cos(\omega \Delta t)+ x_2(t_k)\sin(\omega \Delta t) \bigr),
\end{equation}
\tag{4.23}
$$
$$
\begin{equation}
x_2(t_{k+1}) = E^{\mathrm{I}}(\widehat{t}, a_k) \bigl(x_2(t_k)\cos(\omega \Delta t)- x_1(t_k)\sin(\omega \Delta t) \bigr),
\end{equation}
\tag{4.24}
$$
$$
\begin{equation}
x_3(t_{k+1}) = E^{\mathrm{II}}(\widehat{t}, a_k) x_3(t_k)+\frac{1- E^{\mathrm{II}}(\widehat{t}, a_k)}{1+2a_k},
\end{equation}
\tag{4.25}
$$
where
$$
\begin{equation*}
E^{\mathrm{I}}(\widehat{t}, a_k) := \exp\biggl(- \frac{\gamma \widehat{t} (1+2a_k)}{2 N} \biggr), \qquad E^{\mathrm{II}}(\widehat{t}, a_k) := \exp\biggl(-\frac{\gamma \widehat{t} (1+2a_k)}{N} \biggr).
\end{equation*}
\notag
$$
Equations (4.23)–(4.25) are linear difference equations [72]. With the initial conditions $x_j(t_1=0)=x_{x,j}$, $j=1,2,3$, consider the Cauchy problem for these difference equations. Equations (4.23), (4.24) do not depend on (4.25). If $k=1$, then (4.23), (4.24) give
$$
\begin{equation*}
\begin{aligned} \, x_1(t_2) &= E^{\mathrm{I}}(\widehat{t}, a_1) \bigl(x_1(t_1=0)\cos(\omega \Delta t)+ x_2(t_1=0)\sin(\omega \Delta t) \bigr), \\ x_2(t_2) &= E^{\mathrm{I}}(\widehat{t}, a_1) \bigl(x_2(t_1=0)\cos(\omega \Delta t)- x_1(t_1=0)\sin(\omega \Delta t) \bigr). \end{aligned}
\end{equation*}
\notag
$$
Further, using (4.23) and (4.24) sequentially for $k=2,3$ and substituting the right-hand sides of the equations instead of $x_j(t_k)$, $j=1,2$, we obtain the following non-recurrent (that is, only with $x_j(0)$, $j=1,2$) formulas for $x_1(t_{k+1})$, $x_2(t_{k+1})$, $k=2, 3$:
$$
\begin{equation*}
\begin{aligned} \, x_1(t_3) &= E^{\mathrm{I}}(\widehat{t}, a_2) \bigl(x_1(t_2)\cos(\omega \Delta t)+ x_2(t_2)\sin(\omega \Delta t) \bigr) \\ &= E^{\mathrm{I}}(\widehat{t}, a_1) E^{\mathrm{I}}(\widehat{t}, a_2) \bigl(x_1(0)\cos(2\omega \Delta t)+x_2(0)\sin(2 \omega \Delta t) \bigr), \\ x_2(t_3) &= E^{\mathrm{I}}(\widehat{t}, a_2) \bigl(x_2(t_2)\cos(\omega \Delta t)- x_1(t_2)\sin(\omega \Delta t) \bigr) \\ &= E^{\mathrm{I}}(\widehat{t}, a_1) E^{\mathrm{I}}(\widehat{t}, a_2) \bigl(x_2(0)\cos(2\omega \Delta t)-x_1(0)\sin(2 \omega \Delta t) \bigr), \\ x_1(t_4) &= E^{\mathrm{I}}(\widehat{t}, a_3) \bigl(x_1(t_3)\cos(\omega \Delta t)+ x_2(t_3)\sin(\omega \Delta t) \bigr) \\ &= \prod_{s=1}^3 E^{\mathrm{I}}(\widehat{t}, a_s) \bigl(x_1(0)\cos(3\omega \Delta t)+ x_2(0)\sin(3 \omega \Delta t) \bigr), \\ x_2(t_4) &= E^{\mathrm{I}}(\widehat{t}, a_3) \bigl(x_2(t_3)\cos(\omega \Delta t)- x_1(t_3)\sin(\omega \Delta t) \bigr) \\ &= \prod_{s=1}^3 E^{\mathrm{I}}(\widehat{t}, a_s) \bigl(x_2(0)\cos(3\omega \Delta t)- x_1(0)\sin(3 \omega \Delta t) \bigr). \end{aligned}
\end{equation*}
\notag
$$
Based on these results, we write
$$
\begin{equation}
x_1(t_{k+1}) = \prod_{s=1}^k E^{\mathrm{I}}(\widehat{t}, a_s) \bigl(x_1(0) \cos(k \omega \Delta t)+x_2(0) \sin(k \omega \Delta t) \bigr),
\end{equation}
\tag{4.26}
$$
$$
\begin{equation}
x_2(t_{k+1}) = \prod_{s=1}^k E^{\mathrm{I}}(\widehat{t}, a_s) \bigl(x_2(0) \cos(k \omega \Delta t)-x_1(0) \sin(k \omega \Delta t) \bigr),
\end{equation}
\tag{4.27}
$$
where $k=1, \dots, N$. Using (4.26), (4.27) for $k=N$ with the condition $N \omega \Delta t=\omega \widehat{t}$ and rewriting
$$
\begin{equation*}
\prod_{s=1}^N E^{\mathrm{I}}(\widehat{t}, a_s)=\exp\biggl( -\gamma \widehat{t}\, \biggl(\frac{1}{2}+\frac{1}{N} \sum_{s=1}^N a_s \biggr) \biggr),
\end{equation*}
\notag
$$
we obtain (4.4)–(4.6).
Equation (4.25) does not depend on (4.23), (4.24). Using sequentially (4.25) for $k=1,\dots, 4$, we obtain
$$
\begin{equation*}
\begin{aligned} \, x_3(t_2) &= E^{\mathrm{II}}(\widehat{t}, a_1) x_3(t_1=0)+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_1)}{1+2a_1}, \\ x_3(t_3) &= E^{\mathrm{II}}(\widehat{t}, a_2) x_3(t_2)+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_2)}{1+2a_2} \\ &= E^{\mathrm{II}}(\widehat{t}, a_1) E^{\mathrm{II}}(\widehat{t}, a_2) x_3(0)+\frac{E^{\mathrm{II}}(\widehat{t}, a_2)(1-E^{\mathrm{II}}(\widehat{t}, a_1) )}{1+2a_1} + \frac{1-E^{\mathrm{II}}(\widehat{t}, a_2)}{1+2a_2}, \\ x_3(t_4) &= E^{\mathrm{II}}(\widehat{t}, a_3) x_3(t_3)+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_3)}{1+2a_3} \\ &= \prod_{s=1}^3 E^{\mathrm{II}}(\widehat{t}, a_s) x_3(0)+\frac{E^{\mathrm{II}}(\widehat{t}, a_2) E^{\mathrm{II}}(\widehat{t}, a_3) (1-E^{\mathrm{II}}(\widehat{t}, a_1))}{1+2a_1} \\ &\qquad+ \frac{E^{\mathrm{II}}(\widehat{t}, a_3)(1-E^{\mathrm{II}}(\widehat{t}, a_2))}{1+ 2a_2}+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_3)}{1+2a_3}, \\ x_3(t_5) &= E^{\mathrm{II}}(\widehat{t}, a_4) x_3(t_4)+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_4)}{1+2a_4} \\ &= \prod_{s=1}^4 E^{\mathrm{II}}(\widehat{t}, a_s) x_3(0)+\frac{\prod_{m=2}^4 E^{\mathrm{II}}(\widehat{t}, a_m)(1-E^{\mathrm{II}}(\widehat{t}, a_1))}{1+2a_1} \\ &\qquad +\frac{\prod_{m=3}^4 E^{\mathrm{II}}(\widehat{t}, a_m)(1-E^{\mathrm{II}}(\widehat{t}, a_2))}{1+2a_2} \\ &\qquad + \frac{E^{\mathrm{II}}(\widehat{t}, a_4)(1-E^{\mathrm{II}}(\widehat{t}, a_3))}{1+2a_3}+ \frac{1-E^{\mathrm{II}}(\widehat{t}, a_4)}{1+2a_4}. \end{aligned}
\end{equation*}
\notag
$$
Based on these relations, we write
$$
\begin{equation}
\begin{aligned} \, x_3(t_{k+1}) &= x_3(0) \prod_{s=1}^k E^{\mathrm{II}}(\widehat{t}, a_s) \nonumber \\ &\qquad +\sum_{s=1}^{k-1} \frac{\prod_{m=s+1}^k E^{\mathrm{II}}(\widehat{t}, a_m)(1-E^{\mathrm{II}}(\widehat{t}, a_s))}{1+2 a_s}+\frac{1-E^{\mathrm{II}}(\widehat{t}, a_k)}{1+2 a_k}, \end{aligned}
\end{equation}
\tag{4.28}
$$
where $k=1, \dots, N$. Using (4.28) for $k=N$ and transforming the products $\prod_{s=1}^N E^{\mathrm{II}}(\widehat{t}, a_s)$ and $\prod_{m=s+1}^N E^{\mathrm{II}}(\widehat{t}, a_m)$, we obtain (4.7).
The proof is complete. Proof of Theorem 2. The derivatives (4.11)–(4.14) and (4.19)–(4.22) are easily obtained by the rules of differentiation. To obtain the derivative (4.15)–(4.18), we take into account the difference between the indices in the second term in (4.7). The proof is complete. 4.2.2. Gradient projection method using the exact forms for the objective functions and their gradientsConsider the set $Q$ as $Q=Q_{\infty} := [0, \infty)$ or $Q=Q_{n_{\max}} := [0, n_{\max}]$ with a given $n_{\max}>0$. In finite-dimensional unconstrained optimization, the heavy ball method [65], [66] is known. Also its projection version [64], GPM-2, is known in finite-dimensional constrained optimization. GPM-2 (see the formula (3.3) in [64]) is adapted in the present work and has the following form as applied to the problem of minimizing the objective function $g_1(\mathbf{a})$ taking into account (4.20):
$$
\begin{equation}
\qquad\mathbf{a}^{(0)} \in Q \quad \text{is a given admissible initial vector},
\end{equation}
\tag{4.29}
$$
$$
\begin{equation}
a_k^{(1)} = \overline{a}_k^{\,(0)}(\beta^{(0)}), \qquad k =1,\dots,N,
\end{equation}
\tag{4.30}
$$
$$
\begin{equation}
\overline{a}_k^{\,(0)}(\beta) := \operatorname{Pr}_Q\biggl(a_k^{(0)} - \beta \frac{\partial g_1(\mathbf{a})}{\partial a_k}\bigg|_{\mathbf{a}=\mathbf{a}^{(0)}} \biggr),
\end{equation}
\tag{4.31}
$$
$$
\begin{equation}
a_k^{(m+1)} = \overline{a}_k^{\,(m)}(\beta^{(m)};\lambda), \qquad k =1,\dots,N, \quad m=1,\dots,M,
\end{equation}
\tag{4.32}
$$
$$
\begin{equation}
\overline{a}_k^{\,(m)}(\beta;\lambda) := \operatorname{Pr}_Q \biggl(a_k^{(m)} - \beta \frac{\partial g_1(\mathbf{a})}{\partial a_k}\bigg|_{\mathbf{a}=\mathbf{a}^{(m)}} + \lambda \bigl(a_k^{(m)}-a_k^{(m-1)} \bigr) \biggr),
\end{equation}
\tag{4.33}
$$
$$
\begin{equation}
\operatorname{Pr}_Q(z) = \begin{cases} 0 &\text{if }z<0, \\ z &\text{if }Q=Q_{\infty}\text{ and }z \in Q_{\infty}, \\ z &\text{if } Q=Q_{n_{\max}}\text{ and }z \in Q_{\max}, \\ n_{\max} & \text{if }Q=Q_{\max} \text{ and } z>n_{\max}, \end{cases}
\end{equation}
\tag{4.34}
$$
where the parameters $\beta^{(m)}>0$ ($m=0,\dots,M$) can, as a variant, be sought as fixed for all iterations and providing convergence of the method; also it is required to fix some appropriate value of the parameter $\lambda \in (0, 1)$ for the all iterations. As a condition for the end of computations, consider the inequality $g_1(\mathbf{a})<\varepsilon$ with some sufficiently small $\varepsilon>0$ along with some “ceiling” for numbers of such iterations. The first iteration uses the one-step GPM [66]. By analogy with (4.29)–(4.34), we adapt GPM-2 for minimizing the objective function $g_{\Phi}(\widehat{t}, \mathbf{a}; P')$ taking into account (4.21), (4.22):
$$
\begin{equation*}
\begin{aligned} \, (\widehat{t}, \mathbf{a}^{(0)}) &\in [\widehat{t}_{\min}, \widehat{t}_{\max}] \times Q \quad \text{is a given admissible initial pair}, \\ a_k^{(1)} &= \overline{a}_k^{\,(0)}(\beta^{(0)}), \qquad k =1,\dots, N, \\ \overline{a}_k^{\,(1)}(\beta) &= \operatorname{Pr}_Q\biggl(a_k^{(0)} - \beta \frac{\partial g_{\Phi}(\widehat{t}, \mathbf{a}; P')}{\partial a_k}\bigg|_{(\widehat{t}, \mathbf{a})=(\widehat{t}^{\,(0)}, \mathbf{a}^{(0)})} \biggr), \\ \widehat{t}^{\,(1)} &=\overline{\widehat{t}}^{\,(0)}(\beta^{(0)}), \\ \overline{\widehat{t}}^{\,(0)}(\beta) &= \operatorname{Pr}_{[\widehat{t}_{\min}, \widehat{t}_{\max}]}\biggl(\widehat{t}^{\,(0)} - \beta \frac{\partial g_{\Phi}(\widehat{t}, \mathbf{a}; P')}{\partial \widehat{t}}\bigg|_{(\widehat{t}, \mathbf{a})=(\widehat{t}^{\,(0)}, \mathbf{a}^{(0)})} \biggr), \\ a_k^{(m+1)} &= \overline{a}_k^{\,(m)}(\beta^{(m)},\lambda), \qquad k =1,\dots,N, \quad m=1,\dots,M, \\ \overline{a}_k^{\,(m)}(\beta,\lambda) &= \operatorname{Pr}_Q \biggl(a_k^{(m)} - \beta \frac{\partial g_{\Phi}(\widehat{t}, \mathbf{a}; P')}{\partial a_k}\bigg|_{(\widehat{t}, \mathbf{a})=(\widehat{t}^{\,(0)}, \mathbf{a}^{(0)})} + \lambda \bigl(a_k^{(m)}-a_k^{(m-1)} \bigr) \biggr), \\ \widehat{t}^{\,(m+1)} &= \overline{\widehat{t}}^{\,(m)}(\beta^{(m)},\lambda), \qquad m= 1,\dots, M, \\ \overline{\widehat{t}}^{\,(m)}(\beta,\lambda) &= \operatorname{Pr}_{[\widehat{t}_{\min}, \widehat{t}_{\max}]} \biggl(\widehat{t}^{\,(m)} {-}\, \beta \frac{\partial g_{\Phi}(\widehat{t}, \mathbf{a}; P')}{\partial \widehat{t}}\bigg|_{(\widehat{t}, \mathbf{a})=(\widehat{t}^{\,(0)}, \mathbf{a}^{(0)})} \!{+}\, \lambda \bigl(\widehat{t}^{\,(m)}-\widehat{t}^{\,(m-1)} \bigr) \biggr), \end{aligned}
\end{equation*}
\notag
$$
where $\widehat{t}_{\min}$ and $\widehat{t}_{\max}$ are fixed such that $0<\widehat{t}_{\min}<\widehat{t}_{\max}$; the parameters $\beta^{(m)}>0$ ($m=0,\dots,M$) should be adjusted. Further, in addition, one can restrict variations of piecewise constant controls via the constraints $0 \leqslant (a_{k+1}-a_k)^2 \leqslant \delta_a$, $\delta_a>0$, $k=1, \dots, N-1$. To do this, consider an objective function that includes the function $g_1(\mathbf{a})$ and the additive penalty function $R^{\alpha}(\mathbf{a})$ by the method of external penalty functions (see Chap. 9 in [66]):
$$
\begin{equation*}
g_1^{\alpha}(\mathbf{a})=g_1(\mathbf{a})+R^{\alpha}(\mathbf{a}) \to \inf_a,\qquad R^{\alpha}(\mathbf{a})=\alpha \sum_{k=1}^N \bigl(\max\{ (a_{k+1}-a_k)^2-\delta_a, 0 \} \bigr)^2,
\end{equation*}
\notag
$$
where the cutoff functions are squared for continuity of the partial derivatives, and the parameter $\alpha>0$. The gradient of the penalty function is formed from the partial derivatives
$$
\begin{equation*}
\begin{aligned} \, \frac{\partial R^{\alpha}(\mathbf{a})}{\partial a_1} &=4 \max\{ (a_2-a_1)^2-\delta_a, 0 \}(a_1-a_2), \\ \frac{\partial R^{\alpha}(\mathbf{a})}{\partial a_N} &=4 \max\{ (a_N-a_{N-1})^2-\delta_a, 0 \}(a_N-a_{N-1}), \\ \frac{\partial R^{\alpha}(\mathbf{a})}{\partial a_{k}} &=4 \max\{ (a_k-a_{k-1})^2-\delta_a, 0 \}(a_k- a_{k-1}) \\ &\qquad+4 \max\{ (a_{k+1}-a_k)^2-\delta_a, 0 \}(a_k-a_{k+1}),\qquad 1<k<N. \end{aligned}
\end{equation*}
\notag
$$
For a comparison, in [57] for system (2.3), the time-minimal control problem was considered with the terminal constraint $\rho(T)=\rho_{\mathrm{target}}$. For this problem, the approach was considered with a series of optimal control problems without terminal constraint, with minimizing the square of the Hilbert–Schmidt distance, that is, $\|\rho(T_i)-\rho_{\mathrm{target}} \|^2=\operatorname{Tr}(\rho(T_i)- \rho_{\mathrm{target}})^2$, on descending sequences of fixed final times $\{T_i\}$ under piecewise continuous coherent and incoherent controls, which in the computer implementation of the algorithm were considered as piecewise linear functions. For each such auxiliary optimization problem, the two-parameter one-step GPM with the gradient of the functional was applied (for more general class of optimal control problems with free final state, the theory of optimal control gives more general formula for the corresponding gradient [68], [69], which was specified in [57]). 4.3. Comparing the results of the various first stagesThe values $\omega=1$, $\gamma=0.002$, $\mu=0.01$, and $n_{\max}=100$ are taken. Example 1. Consider the initial state $x_0=(1, 0, 0)$ corresponding to the initial density matrix $\rho_0=\frac{1}{2} \left(\begin{smallmatrix} 1 & 1\\ 1 & 1 \end{smallmatrix}\right)$, and take the target state $x_{\mathrm{target}}=(0, 0, -1/2)$ corresponding to the target density matrix $\rho_{\mathrm{target}}=\left(\begin{smallmatrix} 1/4 & 0\\ 0 & 3/4 \end{smallmatrix}\right)$. The eigenvalues sorted in descending order are $p_1=3/4$ and $p_2=1/4$. As in the original two-stage method, consider the intermediate target density matrix $\widetilde{\rho}_{\mathrm{target}} = \frac{1}{4} |0\rangle \langle 0|+\frac{3}{4} |1\rangle \langle 1| = \left(\begin{smallmatrix} 3/4 & 0\\ 0 & 1/4 \end{smallmatrix}\right)$, the corresponding intermediate target state $\widetilde{x}_{\mathrm{target}}=(0, 0, 1/2)$, and the constant incoherent control $\overline{n}(t) \equiv p_2/(p_1-p_2)=1/2$. The corresponding solution $\overline{x}$ of (4.1) is found from (4.8)–(4.10): $\overline{x}_1(t)=e^{-\gamma t} \cos(\omega t)$, $\overline{x}_2(t)=-e^{-\gamma t} \sin(\omega t)$, $\overline{x}_3(t)=(1/2) (1-e^{-2 \gamma t})$, $t \in [0, \widehat{t}]$. We take into account that the GKSL master equation in the WCL is an approximate for describing relevant physical phenomena. In this regard, when we consider the equation for obtaining $\widehat{t}$, we restrict the consideration to $\varepsilon_{\mathrm{1stage}} \in \{10^{-2}, 10^{-3} \}$. For numerical solving of this equation, the function NSolve in Wolfram Mathematica is used. Thus, for $\varepsilon_{\mathrm{1stage}}=10^{-2}$ and $\varepsilon_{\mathrm{1stage}}=10^{-3}$, we obtain, correspondingly, $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-2}) \approx 2303$ and $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-3}) \approx 3454$. Because when we consider piecewise constant incoherent controls, we expand the class of constant controls, we expect that $\widehat{t}$, which is sufficient for satisfying the condition
$$
\begin{equation}
J_1(n)=g_1(\mathbf{a})= \| x(\widehat{t}\,)-\widetilde{x}_{\mathrm{target}} \|_2^2 \leqslant \varepsilon_{\mathrm{1stage}}^2,
\end{equation}
\tag{4.36}
$$
can be decreased; here $x(\widehat{t}\,)$ is obtained by solving the system with an admissible piecewise constant incoherent control and zero coherent control; $\varepsilon_{\mathrm{1stage}} \in \{10^{-2}, 10^{-3} \}$. For minimizing the objective function $g_1(\mathbf{a})$ for a given $\widehat{t}$, which is less than the values of $\widehat{t}$ obtained above via the unmodified method, we use GPM-2 in the form (4.29)–(4.34). The method was implemented in Python 3 language. We do not pretend to minimize $\widehat{t}$ among such values of $\widehat{t}$ that are sufficient for providing $J_1(n)=g_1(\mathbf{a}) \leqslant \varepsilon_{\mathrm{1stage}}^2=10^{-4}$ or $J_1(n)=g_1(\mathbf{a}) \leqslant \varepsilon_{\mathrm{1stage}}^2=10^{-6}$ with GPM-2. We fix $\widehat{t}=450$ that is significantly less than $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-2}) \approx 2303$ and $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-3}) \approx 3454$. For considering piecewise constant controls, we divide $[0, 450]$ into $N=225$ parts. For visualization of the resulting trajectory, we compute the functions $x_j$, $j=1,2,3$, in a larger number of time instants. In this example, we consider the initial vector $\mathbf{a}^{(0)}=(0, 0, \dots, 0)$ that gives $g_1(\mathbf{a}^{(0)}) \approx 0.42$.  Figure 2.For Example 1. Realization of the modified first stage of the two-stage method: steering from $x_0=(1, 0, 0)$ to the point close to the intermediate target state $\widetilde{x}_{\mathrm{target}}=(0, 0, 0.5)$ for $\widehat{t}=450$ using GPM-2. The figures show the numerically optimized process $(x,n)$: (a) incoherent control $n$ for $t \in [0, 450]$; (b) functions $x_j$ of $t \in [0, 450]$, $j=1,2,3$; (c) trajectory $x$ in the Bloch ball. Here, the distance $\| x(\widehat{t}\,)-\widetilde{x}_{\mathrm{target}} \|_2 \approx 10^{-3}$ For $\widehat{t}=450$ and $\varepsilon_{\mathrm{1stage}}=10^{-3}$, GPM-2, which uses $\beta^{(m)}=10$ $\forall\,m \geqslant 0$ and $\lambda=0.999$, provides the condition $J_1(n)=g_1(\mathbf{a}) \leqslant \varepsilon_{\mathrm{1stage}}^2=10^{-6}$ after $264$ iterations by reaching $g_1(\mathbf{a}) \approx 9.9 \cdot 10^{-7} \approx 10^{-6}$. Thus, the distance $\| x(\widehat{t}\,)-\widetilde{x}_{\mathrm{target}} \|_2 \approx 10^{-3}$. Figure 2 visualizes the results. For comparing, we use GPM-1 (that is, always $\lambda=0$) and after $1000$ iterations we see that $g_1 \approx 0.0269$ which is significantly worse than the result of GPM-2 which uses $\lambda > 0$, while the other algorithmic parameters are the same. Further, for $\widehat{t}=450$ and $\varepsilon_{\mathrm{1stage}}=10^{-2}$, we consider the iterations, which were carried out above with GPM-2 for $\varepsilon_{\mathrm{1stage}}=10^{-3}$, and find that the condition $J_1(n)=g_1(\mathbf{a}) \leqslant \varepsilon_{\mathrm{1stage}}^2=10^{-4}$ is satisfied after $82$ iterations by reaching $g_1(\mathbf{a}) \approx 9.7 \cdot 10^{-5} \approx 10^{-4}$. Moreover, consider $\widehat{t}=400$, $N=200$ for $\varepsilon_{\mathrm{1stage}}=10^{-2}$. Then GPM-2 with $\beta^{(m)}=10$ and $\lambda=0.999$ reaches $g_1(\mathbf{a}) \approx 10^{-4}$ after $120$ iterations. Comparing, on one hand, the results obtained via the modified first stage using GPM-2 and, from the other hand, the results of the unmodified first stage, we see that the modified first stage can reach the same accuracy $\varepsilon_{\mathrm{1stage}}$ significantly faster than the unmodified first stage. If $\varepsilon_{\mathrm{1stage}}=10^{-3}$ and GPM-2 is used when $\widehat{t}= 450$, then the duration of the modified first stage is $\widehat{t}(\varepsilon_{\mathrm{1stage}} = 10^{-3}) / 450 \approx 3454 / 450 \approx 7.7$ times less than the duration of the unmodified first stage. If $\varepsilon_{\mathrm{1stage}}=10^{-2}$ and GPM-2 is used when $\widehat{t}=450$, then $\widehat{t}(\varepsilon_{\mathrm{1stage}}= 10^{-2}) / 450 \approx 2303 / 450 \approx 5.1$ times. If $\varepsilon_{\mathrm{1stage}}=10^{-2}$ and GPM-2 is used when $\widehat{t}=400$, then we find $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-2}) / 450 \approx 2303 / 400 \approx 5.8$ times. However, this acceleration is achieved at the cost of considering the class of piecewise constant controls and the GPM-2 iterations. Example 2 (related to the example from [48], § V). Consider the initial state $x_0=(0, 0, -1)$ (the Bloch ball’s south pole) and the target state $x_{\mathrm{target}}=(0, 0, 1/2)$ corresponding to the target density matrix $\rho_{\mathrm{target}}=\left(\begin{smallmatrix} 3/4 & 0 \\ 0 & 1/4 \end{smallmatrix}\right)$. The intermediate target density matrix is $\widetilde{\rho}_{\mathrm{target}} = \frac{3}{4} |0\rangle \langle 0|+\frac{1}{4} |1\rangle \langle 1| = \left(\begin{smallmatrix} 1/4 & 0\\ 0 & 3/4 \end{smallmatrix}\right)$ that corresponds to the state $\widetilde{x}_{\mathrm{target}}=(0, 0, -1/2)$. In the unmodified first stage of the method, we use the same constant control $\overline{n}=1/2$ that is used in Example 1. Because $x_0$ is different from that in Example 1, here we find the corresponding solution $\overline{x}$ of system (4.1) from (4.8)–(4.10): $\overline{x}_1(t)=0$, $\overline{x}_2(t)=0$, $\overline{x}_3(t)=(1/2)(1-3 e^{-2 \gamma t})$, $t \in [0, \widehat{t}]$. For this solution, consider equation (4.35): For $\varepsilon_{\mathrm{1stage}}=10^{-2}$ and $\varepsilon_{\mathrm{1stage}}=10^{-3}$, we obtain, correspondingly, $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-2}) \approx 98.9$ and $\widehat{t}(\varepsilon_{\mathrm{1stage}}=10^{-3}) \approx 101.1$ in the unmodified first stage.  Figure 3.For Example 2. Realization of the modified first stage: (a) functions $x_j$ of $t$; $j=1,2,3$; (b) trajectory $x$ in the Bloch ball In the modified first stage, consider $\widehat{t}=10$, $N=1$ (that is, without subintervals, see Corollary 1). We use GPM-2 with $\mathbf{a}^{(0)}=a^{(0)}_1=0$, $\beta^{(m)}=1$, and $\lambda= 0.999$. The condition (4.36) with $\varepsilon_{\mathrm{1stage}}=10^{-2}$ is satisfied via GPM-2 after $35$ iterations. The obtained control is $n\,{\equiv}\, 16.205$. This gives $J_1\,{\approx}\, 4 \cdot 10^{-6}\,{<}\, \varepsilon_{\mathrm{1stage}}^2= 10^{-4}$. Considering only $N=1$ significantly simplifies the computations, but at the same time this significantly decreases the first stage’s duration in comparison with $\widehat{t}(\varepsilon_{\mathrm{1stage}}= 10^{-2})$ of the unmodified first stage. Figure 3 shows the results. Example 3. Consider the initial state $x_0=(1, 0, 0)$ as in Example 1 and different target state $x_{\mathrm{target}}=(1/2, 0, 0)$ that corresponds to the target density matrix $\rho_{\mathrm{target}}=\left(\begin{smallmatrix} 1/2 & 1/4 \\ 1/4 & 1/2 \end{smallmatrix}\right)$ whose eigenvalues are $p_1=3/4$ and $p_2=1/4$. The intermediate target density matrix $\rho_{\mathrm{target}}=\left(\begin{smallmatrix} 3/4 & 0 \\ 0 & 1/4 \end{smallmatrix}\right)$, which is considered in Example 1, has the same eigenvalues. In the current example, this matrix $\rho_{\mathrm{target}}$ and the results, which were obtained in Example 1 for the (un)modified first stage, are suitable. 4.4. Realization of the second stage of the two-stage methodExample 4 (in continuation of Example 1). The second stage of the method is implemented in the range $[\widehat{t}, T]$. At the moment $t=\widehat{t}$, we consider as the initial state $x(\widehat{t}\,)$ the intermediate final state $x(\widehat{t}=450)$ computed at the modified first stage using the class of controls (2.7), GPM-2, and $\varepsilon_{\mathrm{1stage}}=10^{-3}$. This state differs only very little from the intermediate target state $\widetilde{x}_{\mathrm{target}}=(0, 0, 0.5)$. At the second stage, with zero incoherent control ($n=0$), we look for such final moment $T$ and coherent control $v$ that allow one to obtain the final state $x(T)$ close to the given target state $x_{\mathrm{target}}=(0, 0, -0.5)$ with some suitable accuracy.  Figure 4.For Example 4. In continuation of Example 1, realization of the second stage of the two-stage method: steering from the point, which is close to the intermediate target state $\widetilde{x}_{\mathrm{target}}\,{=}\,(0, 0, 0.5)$, to some point, which is close to the target state $x_{\mathrm{target}}=(0, 0, -0.5)$, during the time $T-\widehat{t}=4.99$ under zero incoherent control and coherent control $v=-61.8 \sin(2\pi (t- 450)/4.99)$ of the type (2.9). We show the resulting (a) control $v$; (b) functions $x_j$ of $t \in [450, 454.99]$; $j=1,2,3$; (c) trajectory $x$ in the Bloch ball The implementation of the second stage is performed in analogy with the example from [48]. Namely, we consider coherent control $v$ of the type (2.8), at that we use the constraint $|A| \leqslant \nu=100$ for the amplitude, and we set the frequency $\Omega=1$. The incoherent control is zero. The objective functional is $J_{2, \mathrm{trig.}}(v)=\| x(T)-x_{\mathrm{target}} \|_2=\bigl( \sum_{j=1}^3 (x_j(T)- x_{{\mathrm{target}},j})^2 \bigr)^{1/2} $. The final time $T$ is defined as described below. On the range $[-\nu, \nu]=[-100, 100]$, we introduce a uniform grid $\{A_j\}$ with step $\Delta A=0.05$. For each node $A_j$: 1) we numerically solve the dynamical system with the controls $v=A_j \cos t$, $n=0$; 2) the corresponding vector function $x$ is stored on a uniform grid with nodes $\{t_k\}$ in the range $[\widehat{t}, \overline{T}]$ with step $\Delta t=0.01$, where $\overline{T}$ is some sufficiently large number. After that, we look for such $\{t_k\}$ and $\{A_j\}$ that the distance $\| x(t_k)-x_{\mathrm{target}} \|_2 \leqslant \varepsilon_{\mathrm{2stage}}$ with some small $0< \varepsilon_{\mathrm{2stage}} \ll 1$, if, of course, such time nodes exist. We take $\varepsilon_{\rm 2stage}=10^{-2}$. From the set of such nodes $\{t_k\}$, we select the node with the smallest value. We set $\overline{T}= 40$, that is, the time grid is introduced in the range $[450, 490]$. By the described way, we obtain the amplitude $A=A_{j^{\ast}}=-67.6$ and the final time $T=t_{k^{\ast}}=455.38$. The first stage duration is $\widehat{t}=450$, the second stage duration is $455.38-450= 5.38$. Here, $\|x(T)-x_{\mathrm{target}}\|_2 \approx 4.5 \cdot 10^{-3}<0.5\varepsilon_{\mathrm{2stage}}$. Note that $v(\widehat{t}\,) \approx 49$ and $v(T) \approx 67$ here. Further, consider the condition $v(\widehat{t}\,)=v(T)=0$ and the class of controls (2.9). With $\varepsilon_{\mathrm{2stage}}=10^{-2}$, we obtain $d=2$, $T=454.99$, and $A=-61.8$ that give the distance $\|x(T)-x_{\mathrm{target}}\|_2 \approx 5 \cdot 10^{-3} \approx 0.5\varepsilon_{\mathrm{2stage}}$. Figure 4 illustrates the resulting control process for the second stage. Another version for realization of the second stage can be the following. With a given $T$, consider the objective functional
$$
\begin{equation*}
J_2^{\alpha}(v)=\| x(T)-x_{\mathrm{target}} \|^2_2+\alpha \int_{\widehat{t}}^T \frac{v^2(t)}{S(t)}\, dt
\end{equation*}
\notag
$$
to be minimized, where $v$ is piecewise continuous coherent control; $x$ is the solution of (4.1) with an admissible $v$ and $n=0$; the integral term with $\alpha>0$, the function
$$
\begin{equation*}
S(t)=\exp\biggl(-b \biggl(\frac{t- \widehat{t}}{T-\widehat{t}}-\frac{1}{2} \biggr)^2 \biggr)
\end{equation*}
\notag
$$
with $b>0$ is for approximate providing the condition $v(\widehat{t}\,)=v(T)=0$ (this term is taken by analogy with such the possibility in the Krotov method’s technique [17]). With respect to numerical optimization of $v$ via a gradient method, consider the gradient of $J_2^{\alpha}(v)$ using the general form for such gradient known in the theory of optimal control [68], [69]. The increment formula for $J_2^{\alpha}(v)$ on a given $v^{(m)}$ ($m \geqslant 0$) and arbitrary $v$ is
$$
\begin{equation*}
\begin{aligned} \, J_2^{\alpha}(v)-J_2^{\alpha}(v^{(m)}) &=\int_{\widehat{t}}^T \bigl\langle \operatorname{grad} J_2^{\alpha}(v^{(m)})(t),\, v(t)-v^{(m)}(t) \bigr\rangle\, dt +r, \\ \operatorname{grad} J_2^{\alpha}(v^{(m)})(t) &=-2 \mu \bigl(p_3^{(m)}(t) x_2^{(m)}(t)-p_2^{(m)}(t) x_3^{(m)}(t) \bigr)+2\alpha \frac{v^{(m)}(t)}{S(t)}. \end{aligned}
\end{equation*}
\notag
$$
Here, the function $x^{(m)}$ is the solution of (4.1) with the initial state $x(\widehat{t}\,)$ and controls $v = v^{(m)}$ and $n=0$, and $p^{(m)}$ is the solution of the conjugate system
$$
\begin{equation*}
\frac{dp^{(m)}(t)}{dt}=-\bigl(A^\top+(B^v)^\top v^{(m)}(t) \bigr) p^{(m)}(t),\qquad p^{(m)}(T)=-2(x^{(m)}(T)-x_{\mathrm{target}}).
\end{equation*}
\notag
$$
§ 5. ConclusionsThe present paper considers control problems for some closed and open two-level quantum systems. The dynamics of the closed system is determined by the Schrödinger equation with coherent control, while the dynamics of the open system is determined by the GKSL master equation whose Hamiltonian depends on coherent control and superoperator of dissipation depends on incoherent control. For the closed system, the problem of generating the phase shift quantum gate for some set of phases and final times is studied. It is shown numerically (via the GRAPE-type approach and the approach using the dual annealing and differential evolution methods) that zero coherent control, which is a stationary point of the objective functional, is not optimal; that gives an example of subtle point for practice in solving quantum control problems. The search based on the exact analytical expression for the gradient is a suitable choice for revealing fine details of the quantum control landscape. For the open system, a modified version of the two-stage method which was proposed and developed for approximate generation of target density matrices for generic $N$-level quantum systems in [48], in this work is studied for two-level systems with the goal of decreasing time necessary for the first (incoherent) stage: at the first stage, instead of constant incoherent control, which is explicitly specified using eigenvalues of the target density matrix, numerical optimization of piecewise constant incoherent control is carried out. Exact analytical formulas are obtained for the quantum system state at the end of the first stage, the objective functions and their gradients depending on the parameters of piecewise constant controls. These formulas are then used in the two-step GPM. For some values of the system parameters, the initial and target density matrices, and various precision for the distance, in the numerical experiments it is shown that the durations of the modified first stage can be significantly less (up to several times) than that of the unmodified first stage, but at the cost of optimizing in the class of piecewise constant controls and using time-dependent incoherent controls. If the simplicity and efficiency of obtaining the constant incoherent control is more important than the decrease in the duration of the first stage, then the unmodified method should be used. The choice of the unmodified method can also be caused by the fact that the constant incoherent control is simpler than some piecewise constant control. Otherwise if there is no experimental difficulty in working with piecewise constant controls and with applying an appropriate optimization method (for example, GPM-2) and if it is necessary to decrease the duration of the first stage, then one can use the modified version.
Citation in format AMSBIB
\Bibitem{MorPec23} |
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Contact us: |
Terms of Use
|
Registration to the website |
Logotypes |
|