††thanks: Electronic address<EMAIL_ADDRESS>Electronic
address<EMAIL_ADDRESS>
# Einstein-Podolsky-Rosen steering based on semi-supervised machine learning

Lifeng Zhang, Zhihua Chen  School of Science, Jimei University, Xiamen
361021,China  Shao-Ming Fei  School of Mathematical Sciences, Capital Normal
University, Beijing 100048, China,
Max Planck Institute for Mathematics in the Sciences, 04103 Leipzig, Germany

###### Abstract

Einstein-Podolsky-Rosen(EPR)steering is a kind of powerful nonlocal quantum
resource in quantum information processing such as quantum cryptography and
quantum communication. Many criteria have been proposed in the past few years
to detect the steerability both analytically and numerically. Supervised
machine learning such as support vector machines and neural networks have also
been trained to detect the EPR steerability. To implement supervised machine
learning, one needs a lot of labeled quantum states by using the semidefinite
programming, which is very time consuming. We present a semi-supervised
support vector machine method which only uses a small portion of labeled
quantum states in detecting quantum steering. We show that our approach can
significantly improve the accuracies by detailed examples.

##  I Introduction

Einstein-Podolsky-Rosen (EPR) steering is an intermediate quantum correlation
between quantum entanglement and quantum nonlocality Wiseman , wherein one
party can steer the state of another distant party by making measurements on
one party of a shared bipartite state. EPR-steering was first observed by
Schrödinger in the famous EPR paradox in 1935 Schrodinger ; Einstein .
However, the rigorous definition of EPR-steering was not proposed till 2007.
EPR-steering is proven to be asymmetric in general, that is, one party can
steer the another but not vice versa Bowles ; A-V-N ; Bowles2 . It has been
shown that steering has many applications in quantum information processing
such as one-sided device independent quantum key distribution Branciard ,
randomness certification Passaro ; Coyle , subchannel discrimination Piani ;
Sun , secret sharing Xiang , quantum teleportation He , coupled qubits and
magnetoreception Ku , no-cloning of quantum steering Chiu  and quantum
networks Chen .

Many methods have been proposed to detect EPR-steering, such as linear and
nonlinear steering inequalities ineq1 ; ineq2 ; ineq3 ; ineq4 ; ineq5 ; ineq6
, uncertainty relations unc1 ; unc2 ; unc3 , moment matrix approach mom  and
all-versus-nothing method A-V-N . Many steering measures were also proposed,
which can be estimated based on semidefinite programming SDP . Nevertheless,
it is generally very time consuming to detect steerability of a state
numerically, as a huge amount of measurement directions has to be run over.
Efficient criteria of steerability are still far from being satisfied.

Machine learning is a useful tool for classification problems, which has
already successful applications in quantum information processing tasks such
as entanglement classification, nonlocality discriminant, Hamiltonian learning
and phase transition identification svm1 ; svm2 ; svm3 ; svm4 ; svm5 ; svm6 .
Machine learning methods such as support vector machine (SVM), artificial
neural networks and decision trees have been also applied to the EPR-steering
detection and quantification Changliang ; Zhang . However, these methods are
supervised machine learning ones, in which a lot of labeled results are
needed. While quantum state labeling is also a very time consuming task.

Semi-supervised machine learning was proposed to combine labeled and unlabeled
data to improve the learning behavior Zhu X . Semi-supervised support vector
machine combines support vector machine and semi-supervised learning, which
can make use of the scarce labeled data and sufficient unlabeled data to
improve the classification performance Ding . The safe semi-supervised SVM
(S4VM) was proposed to exploit multiple candidate low-density separators and
to select the representative separators, which is never significantly inferior
to the inductive SVM S4VM . Semi-supervised SVMs have been successfully
applied to text and image classification, face recognition and many other
areas ssvm1 ; ssvm2 ; ssvm3 ; ssvm4 .

In this paper, we apply semi-supervised SVM-S4VM to deal with quantum steering
problems. We need only a small portion of labeled states generated randomly by
the semidefinite programming, together with a large portion of unlabeled
quantum states used in S4VM. Compared with the inductive SVM, the accuracies
can be significantly improved. Moreover, less time is needed to label the
quantum states, while high accuracy is still attained.

##  II Supervised Machine Learning of EPR-Steering

A two-qubit quantum state,
$\rho_{AB}=\frac{1}{4}(I_{4}+\sum\limits_{i=1}^{3}r_{i}\sigma_{i}\otimes
I_{2}+\sum\limits_{i=1}^{3}s_{i}I_{2}\otimes\sigma_{i}+\sum\limits_{i,j=1}^{3}t_{ij}\sigma_{i}\otimes\sigma_{j})$,
where $I_{d}$ is the $d\times d$ identity matrix and $\sigma_{i}$ ($i=1,2,3$)
are the standard Pauli matrices, is said to admit a local hidden state (LHS)
model if the probability $P(a,b|\textbf{A},\textbf{B})$ that Alice and Bob
performs the measurements A and B with measurement outcomes $a$ and $b$,
respectively, satisfy

$\displaystyle P(a,b|\textbf{A},\textbf{B})=$
$\displaystyle\rm{tr}[(M_{\textbf{A}}^{a}\otimes
M_{\textbf{B}}^{b}).\rho_{AB}]$ $\displaystyle=$ $\displaystyle\sum
p(\lambda)p(a|\textbf{A},\lambda)p_{Q}(b|\textbf{B},\lambda)$ (1)

where $P_{Q}(b|\textbf{B},\lambda)=\rm{tr}[\rho_{\lambda}M_{\textbf{B}}^{b}]$,
$\rho_{\lambda}$ are qubit states specified by the parameter $\lambda$. State
$\rho_{AB}$ is said to be not steerable if it satisfies relation (II).
Otherwise we say that $\rho_{AB}$ is steerable from Alice to Bob.
Equivalently, one may say that Alice can not steer Bob if there exist
probability $p_{\lambda}$ and a set of quantum states $\rho_{\lambda}$ such
that

$\rho^{a}_{\textbf{A}}=\rm{tr}_{A}((M^{a}_{\textbf{A}}\otimes
I_{2}).\rho_{AB})=\sum p_{\lambda}p(a|\textbf{A},\lambda)\rho_{\lambda},$ (2)

where $p(a|\textbf{A},\lambda)$ is a probability distribution given by $a$ and
$\lambda$.

For any $\rho_{AB}$ if Alice performs $m$ measurements
$\textbf{A}=\\{0,1,2,\cdots,m-1\\}$ with $q$ outcomes
$a\in\\{1,2,\cdots,q\\}$, an assemblage
$\\{\rho_{\textbf{A}}^{a}\\}_{\textbf{A},a}$ of $m$ ensembles is obtained.
Denote $D(a|\textbf{A},\lambda)=\delta_{a,\lambda(\textbf{A})}$, that is,
$D(a|\textbf{A},\lambda)=1$ if $\lambda(\textbf{A})=a$ and
$D(a|\textbf{A},\lambda)=0$ if $\lambda(\textbf{A})\neq a$. The following SDP
algorithm can be used to detect if Alice can steer Bob for any given
$\rho^{a}_{\textbf{A}}$ and $D(a|\textbf{A},\lambda)$ SDP ,

$\displaystyle\min\limits_{F_{a|\textbf{A}}}\rm{tr}\sum\limits_{a,\textbf{A}}F_{a|\textbf{A}}\rho^{a}_{\textbf{A}}$
(3) $\displaystyle\rm{such}\
\rm{that}.\sum\limits_{a,\textbf{A}}F_{a|\textbf{A}}D(a|\textbf{A},\lambda)\geq
0,$
$\displaystyle\rm{tr}(\sum\limits_{a,\textbf{A},\lambda}F_{a|\textbf{A}}D(a|\textbf{A},\lambda))=1.$

If the objective value of (3) is negative for some measurements $\textbf{A},$
$\rho_{AB}$ is steerable from Alice to Bob. Otherwise, there exists an LHS
model.

A state $\rho_{AB}$ is labeled $-1$ if the objective value is negative after
using the SDP program 100 times with different values of measurements.
Otherwise, the state is labeled $+1$. In Changliang  a SVM is used for the
classifier. For $m=2,3,\cdots,8$, $5000$ samples with label $+1$ and $5000$
samples with $-1$ were obtained. The last $1000$ positive samples and $1000$
negative samples were reserved for tests, the remaining $4000$ positive and
$4000$ negative samples were kept as training set to learn a classifier.

SVM is a supervised learning model used for classification and regression
analysis. Given $(\textbf{x}_{i},y_{i})$, $i=1,2,\cdots,n$, where $n$ is the
number of samples, $\textbf{x}_{i}$ $(i=1,2,\cdots,n)$ are the feature vectors
of the samples and $y_{i}$ $(i=1,2,\cdots,n)$ are the labels of the samples,
$y_{i}=1$ or $y_{i}=-1$. For linearly separable problems, SVM aims to find a
hyperplane with the largest distance to the nearest feature vectors of each
class. The hyperplanes $H$ are defined as $\omega\textbf{x}_{i}+b\geq 1$ when
$y_{i}=1$, and $\omega\textbf{x}_{i}+b\leq-1$ when $y_{i}=-1$. Two marginal
hyperplanes are defined by $H_{1}:$ $\omega\textbf{x}_{i}+b=1$ and $H_{2}:$
$\omega\textbf{x}_{i}+b=-1.$ The vectors on $H_{1}$ and $H_{2}$ are the
support vectors. The distance $\frac{2}{\|\omega\|}$ between $H_{1}$ and
$H_{2}$ needs to be maximized. The final separating hyperlane $H$ is selected
to be the middle one between the two marginal hyperplanes.

Hence, the problems to be solved are given as follows,

$\displaystyle\min\limits_{\omega}\frac{1}{2}\omega^{T}\omega$ (4)
$\displaystyle\rm{such}\ \rm{that}$ $\displaystyle
y_{i}(\omega\textbf{x}_{i}+b)\geq 1.$

When the errors $\\{{\xi_{i}}\\}_{i=1}^{n}$ in classification are allowed, the
soft marginal hyperplane problem is attained,

$\displaystyle\min\limits_{\omega,b,\xi}\frac{1}{2}\omega^{T}\omega+C\sum\limits_{i=1}^{n}\xi_{i}$
(5) $\displaystyle\rm{such}\ \rm{that}$ $\displaystyle
y_{i}(\omega\textbf{x}_{i}+b)\geq 1-\xi_{i},$ $\displaystyle\xi_{i}\geq 0,$

where $\xi_{i}$ are slack variables, $\xi_{i}=0$ if there is no error for
$\textbf{x}_{i}$, $C$ is a tradeoff parameter between the error and the
margin. For non-linearly separable problems, by using the kernel trick the
problem needing to be solved becomes

$\displaystyle\min\limits_{\omega,b,\xi}\frac{1}{2}\omega^{T}\omega+C\sum\limits_{i=1}^{n}\xi_{i}$
(6) $\displaystyle\rm{such}\ \rm{that}$ $\displaystyle
y_{i}(\omega\phi(\textbf{x}_{i})+b)\geq 1-\xi_{i},$ $\displaystyle\xi_{i}\geq
0,$

where $\phi(\textbf{x})$ is the kernel function svmp1 ; svmp2 .

Concerning the steering detection problem, the training set
$\\{\textbf{x}_{i},y_{i}\\}$ $(i=1,2,\cdots,8000)$ and the radial basis
function
$K(\phi(\textbf{x}_{i}),\phi(\textbf{x}_{j}))=e^{-\gamma\|\textbf{x}_{i}-\textbf{x}_{j}\|^{2}}$
are used with the parameters $C$ and $\gamma$ determined by a grid search
approach. The classifier is attained by solving (6) Changliang .

##  III Semi-supervised Machine Learning of EPR-Steering

A well-trained SVM needs a lot of labeled samples. For steering detection it
takes much time to generate the labeled samples and to obtain the data set
with the increase of measurements Changliang . Semi-supervised SVM can be used
in classifying problems with a small portion of labeled samples and a large
number of unlabeled samples.

The S3VM was first proposed to simultaneously learn the optimal hyperplane and
the labels for unlabeled instances s3vm ; Zhu . Given a set of $l$ labeled
data $\\{\textbf{x}_{i},y_{i}\\}_{i=1}^{l}$ and a set of $u$ unlabeled data
$\\{\hat{x}_{j}\\}_{j=l+1}^{l+u}$, $y_{i}\in\\{1,-1\\}$, one needs to solve
the following optimization problem:

$\displaystyle\min\limits_{\omega,b,\hat{y}\in{\emph{B}}}(\frac{1}{2}\|\omega\|^{2}+C_{1}\sum\limits_{i=1}^{l}\xi_{i}+C_{2}\sum\limits_{j=l+1}^{l+u}\hat{\xi}_{j})$
(7) $\displaystyle\rm{such}\ \rm{that}\quad
y_{i}(\omega^{\prime}\phi(\textbf{x}_{i})+b)\geq 1-\xi_{i},~{}~{}\xi_{i}\geq
0$ $\displaystyle\hat{y}_{j}(\omega^{\prime}\phi(\hat{\textbf{x}}_{j})+b)\geq
1-\hat{\xi}_{j},~{}~{}\hat{\xi}_{j}\geq 0$ $\displaystyle\forall
i=1,\cdots,l,~{}\forall j=l+1,\cdots,l+u,$

where $\emph{B}=\\{\hat{y}\in\\{\pm
1\\}^{u}|-\beta\leq\frac{\sum\limits_{j=l+1}^{l+u}\hat{y}_{j}}{u}-\frac{\sum\limits_{i=1}^{l}y_{i}}{l}\leq\beta\\}$
is the balanced constraint to avoid the trivial solution.

Unlike S3VM, the S4VM was proposed to construct diverse large-margin
separators since multiple large-margin low density separators may coincide
with the limited labeled data, and then the label assignment for unlabeled
instances is optimized such that the worst-case performance improvement over
inductive SVM is maximized S4VM . Recalling the S4VM in detail, one first has
the following minimization problem:

$\displaystyle\min\limits_{\\{\omega_{t},b_{t},\hat{y}_{t}\in{\emph{B}}\\}_{t=1}^{T}}\sum\limits_{t=1}^{T}(\frac{1}{2}\|\omega_{t}\|^{2}+C_{1}\sum\limits_{i=1}^{l}\xi_{i}+C_{2}\sum\limits_{j=l+1}^{l+u}\hat{\xi}_{j})$
(8) $\displaystyle+G\sum\limits_{1\leq t\neq\tilde{t}\leq
T}\delta(\frac{\hat{y}^{\prime}_{t}\hat{y}_{\tilde{t}}}{u}\geq 1-\varsigma),$
$\displaystyle\rm{such}\ \rm{that}\
y_{i}(\omega_{t}^{\prime}\phi(\textbf{x}_{i})+b_{t})\geq
1-\xi_{i},~{}~{}\xi_{i}\geq 0$
$\displaystyle\hat{y}_{t,j}(\omega_{t}^{\prime}\phi(\hat{\textbf{x}}_{j})+b_{t})\geq
1-\hat{\xi}_{j},~{}~{}\hat{\xi}_{j}\geq 0$ $\displaystyle\forall
i=1,\cdots,l,~{}\forall j=l+1,\cdots,l+u,~{}\forall t=1,\cdots,T.$

where $\emph{B}=\\{\hat{y}_{t}\in\\{\pm
1\\}^{u}|-\beta\leq\frac{\sum\limits_{j=l+1}^{l+u}\hat{y}_{t,j}}{u}-\frac{\sum\limits_{i=1}^{l}y_{i}}{l}\leq\beta\\},$
$\hat{y}_{t,j}$ is the $jth$ entry of $\hat{y}_{t},$ $\delta$ is the indicator
function which represents a quantity of penalty about the diversity of
separators and $\varsigma\in[0,1]$ is a constant, $G$ is a large constant
enforcing large diversity. $C_{1}$ and $C_{2}$ are regularization parameters
trading off the complexity and the empirical error on label and unlabeled
data. $T$ is the number of the separators.

Second, from the global simulated annealing search with a deterministic local
search scheme method or the sampling strategy used in solving minimization
problems Eq.(8), we have the following problem,

$\displaystyle\bar{y}=arg\max\limits_{y\in\\{\pm
1\\}^{u}}\min\limits_{\hat{y}\in M_{0}}J(y,\hat{y},y^{\rm{SVM}}),$ (9)

where $y^{\rm{SVM}}$ is the predictive labels of inductive SVM on unlabeled
instances, $M_{0}=\\{\hat{y}_{t}\\}_{t=1}^{T}$ is obtained from (8),

$\displaystyle J(y,\hat{y},y^{\rm{SVM}})=$
$\displaystyle\rm{gain}(y,\hat{y},y^{\rm{SVM}})-\lambda\rm{loss}(y,\hat{y},y^{\rm{SVM}})$
$\displaystyle=$ $\displaystyle c_{t}^{\prime}y+d_{t}$ (10)

with

$\displaystyle\rm{gain}(y,\hat{y},y^{\rm{SVM}})=\sum\limits_{j=l+1}^{l+u}\frac{1+y_{j}\hat{y}_{j}}{2}\frac{1-y^{\rm{SVM}}_{j}\hat{y}_{j}}{2},$
(11)
$\displaystyle\rm{loss}(y,\hat{y},y^{\rm{SVM}})=\sum\limits_{j=l+1}^{l+u}\frac{1-y_{j}\hat{y}_{j}}{2}\frac{1+y^{\rm{SVM}}_{j}\hat{y}_{j}}{2},$
$\displaystyle
c_{t}=\frac{1}{4}[(1+\lambda)\hat{y}_{t}+(\lambda-1)y^{\rm{SVM}}],$
$\displaystyle
d_{t}=\frac{1}{4}[-(1+\lambda)\hat{y}^{\prime}_{t}y^{\rm{SVM}}+(1-\lambda)u].$

Then the minimization problem in (9) can be reformulated as the following
maximization problem,

$\displaystyle\max\limits_{y,\tau}\tau\hskip 128.0374pt$ (12)
$\displaystyle\rm{such}\ \rm{that}$ $\displaystyle\tau\leq
c_{t}^{\prime}y+d_{t},~{}\forall t=1,2,\cdots,T;~{}y\in\\{\pm 1\\}^{u}.$

Convex linear programming problems can be solved by relaxing $y\in\\{\pm
1\\}^{u}$ to $y\in[-1,1]^{u},$ and then project it back to integer solutions
with minimum distance. If the function value of the resulting integer solution
is smaller than that of $y^{\rm{SVM}}$, $y^{\rm{SVM}}$ is output as the final
solution instead.

For any two-qubit quantum state $\rho_{AB},$ set
$\rho_{0}=({I}_{2}\otimes\rho_{B}^{\frac{1}{2}})\rho_{AB}({I}_{2}\otimes\rho_{B}^{\frac{1}{2}}),$
where $\rho_{B}=\rm{tr}_{A}\rho_{AB}$. The coefficients
$\tau_{kl}=\rm{tr}(\rho_{0}(\sigma_{k}\otimes\sigma_{l}))$ can be reformulate
as a nine-dimensional feature vector
$\textbf{x}=[\tau_{11},\tau_{12},\cdots,\tau_{33}]^{T}$, where $T$ denotes
transpose. Generating randomly quantum states by SDP SDP ; Changliang , we get
a set of $l$ feature vectors of labeled quantum states,
$\\{\textbf{x}_{i}\\}$, $i=1,2,\cdots,l$, with $y_{i}=-1$ for steerable states
and $y_{i}=1$ for unsteerable states. Denote
$\\{\hat{\textbf{x}}_{j}\\}_{j=l+1}^{l+u}$ the set of feature vectors of
unlabeled quantum states whose labels can be determined by the above S4VM S4VM
method. Radial basis function kernels are used here and
$K(\phi(\textbf{x}_{i}),\phi(\textbf{x}_{j}))=e^{-\gamma\|\textbf{x}_{i}-\textbf{x}_{j}\|^{2}}$.
We have three hyperparameters $C_{1}$, $C_{2}$, and $\gamma$ which are also
determined by the grid research method.

##  IV Numerical Results

We generate randomly quantum states and get the class of steerable states from
such random states by SDP. For balance, we choose half positive quantum states
and half negative random states.

We implement SVMs by using $l$ labeled quantum states, the radial basis
function kernel, ten-folded cross validation, and the grid search approach.
$u$ unlabeled quantum states are considered as the test sets and the test
errors are obtained. For S4VM, sample strategy is implemented here, three
hyperparameters $C_{1}$, $C_{2}$, and $\gamma$ are determined by the grid
search approach. We set $\beta=0.1,\lambda=3,$ the number of clusters $T$ is
10 and the number of samples is 100 in the sample strategy for the examples.
$u$ unlabeled quantum states are divided into $M$ different sets. Taking $M=2$
as an example, we have the following steps:

$(1)$ The first set of unlabeled quantum states can be predicted by using $l$
labeled quantum states. Then, we implement five-fold cross validation for
these $l+u/M$ quantum states by using SVMs. The best cross-validation accuracy
and the best hyperparameters are then obtained. The labels of the first set of
unlabeled quantum states from the best hyperparameters can be regarded as the
real labels and the classification accuracies of these unlabeled quantum
states are obtained.

$(2)$ After labeling the first set of unlabeled quantum states, $l+u/M$
labeled quantum states and the second set of unlabeled quantum states can be
utilized to implement S4VM the second time. Then the second set of unlabeled
quantum states can also be labeled. Implementing five-fold cross validation
for these $l+2*u/M,$ quantum states by using SVMs, the best cross-validation
accuracy and the best hyperparameters are obtained. The labels of the second
unlabeled quantum states from the best hyperparameters are considered as the
real labels and the classification accuracies of the second set of the
unlabeled quantum states are obtained.

$(3)$ The average accuracy of the two sets of unlabeled quantum states is
considered as the classification accuracy of the unlabeled set.

We compute the errors of $4000$ unlabeled quantum states when $30$ labeled
quantum states are used to implement S4VM for $m=2$ and $M=1,2,4,8$. Except
for the third labeled set, the errors are smaller when $M=1$ and $M=2$
compared to the errors when $M=4$ and $M=8$ for the most cases but the errors
when $M=2$ are smaller than that when $M=1$ for $m=8$, see Figs. 1 and 2. In
addition, the computational speed for $M=2$ is three times faster than that
for $M=1.$

Figure 1: Classification errors of the unlabeled quantum states by using ten
different sets of $l$ labeled quantum states for $m=2,$ $l=30,$ and $u=4000$.
The errors of the unlabeled states by S4VM for $M=1,2,4,8$ are represented by
solid line with star, thick dashed line with x, dashed line with square and
solid line with $+$, respectively.  Figure 2: Classification errors of the
unlabeled quantum states by using ten different sets of 30 labeled quantum
states for $m=8$ and $u=4000$. The errors of the unlabeled states by S4VM for
$M=1$ and $2$ are represented by solid line with star and thick dashed line
with x, respectively.

Based on the above errors and computational speed, in the following numerical
experiments we set $M=2$ and use $10$, $30$, or $50$ labeled states and $4000$
unlabeled quantum states to implement S4VM for $m=2,4,6,8$. When Alice
performs $m$ $(m=2,4,6,8)$ measurements, S4VM is implemented ten times by
using ten different sets of $l$ labeled quantum states. The errors for these
$u$ unlabeled quantum states are shown in the Figs. 3-5.

Figure 3: Classification errors of SVM and S4VM by using ten different sets of
$l$ labeled quantum states for $m=2,4,6,$ and $8,$ $l=10,$ $u=4000$ and $M=2$.
The accuracies of the unlabeled states by SVM and S4VM are represented by the
lines with squares and $\star$, respectively. Figure 4: Classification errors
of SVM and S4VM by using ten different sets of $l$ labeled quantum states for
$m=2,4,6,$ and $8,$ $l=30,$ $u=4000$ and $M=2$. The accuracies of the
unlabeled states by SVM and S4VM are represented by the lines with squares and
$\star$, respectively,  Figure 5: Classification errors of SVM and S4VM by
using ten different sets of $l$ labeled quantum states for $m=2,4,6,$ and $8,$
$l=50,$ $u=4000$ and $M=2$. The accuracies of the unlabeled states by SVM and
S4VM are represented by lines with squares and $\star$, respectively.

When two measurements are performed and ten labeled quantum states are used,
the errors by SVM are larger than $0.2$ in many cases, which can be reduced to
about 0.05 by S4VM in most cases. When $l=30$, the errors by S4VM are less
than 0.06, while the errors by SVM are larger than 0.09 in most cases. When
$l=50,$ the errors by S4VM are less than 0.05, while the errors by SVM are
larger than 0.07 in most cases. In Table 1, we list the maximum differences
between errors from SVM and S4VM by using 10, 30, and 50 labeled quantum
states and 4000 unlabeled quantum states for $m=2,\,4,\,6$, and 8. It can be
seen that the accuracies are significantly improved.

Table 1: The maximum differences between the errors from SVM and S4VM by using ten different $l$ labeled quantum states and 4000 unlabeled quantum states for $l=10,30,$ and $50$ and $m=2,4,6,$ and $8$. $\Delta E_{\rm{max,l}}$ is the maximum difference, the errors from SVM minus the errors from S4VM, by using $l$ labeled quantum states. $m$ | $\Delta E_{\rm{max,10}}$ | $\Delta E_{\rm{max,30}}$ | $\Delta E_{\rm{max,50}}$
---|---|---|---
2 | 0.216 | 0.129 | 0.055
4 | 0.22 | 0.088 | 0.159
6 | 0.36 | 0.12 | 0.13
8 | 0.104 | 0.028 | 0.01

To investigate the performance on positive and negative quantum states by
S4VM, we show the errors of positive class and negative class for $l=10$, $30$
and $50$ in FIGs. 6, 7 and 8.

Figure 6: Classification errors of $2000$ positive quantum states and $2000$
negative quantum states by SVM and S4VM with ten different labeled sets for
$l=10,$ $m=2,4,6,$ and $8$. The errors for positive and negative quantum
states by S4VM are represented by lines with squares and $\star$,
respectively, while the errors for positive and negative quantum states by SVM
are represented by lines with diamond and $+$, respectively.  Figure 7:
Classification errors of $2000$ positive quantum states and $2000$ negative
quantum states by SVM and S4VM with ten different labeled sets for $l=30,$
$m=2,4,6,$ and $8$. The errors for positive and negative quantum states are
represented by lines with squares and $\star$, respectively, while the errors
for positive and negative quantum states by SVM are represented by lines with
diamond and $+$, respectively.  Figure 8: Classification errors of $2000$
positive quantum states and $2000$ negative quantum states by SVM and S4VM
with ten different labeled sets for $l=50,$ $m=2,4,6,$ and $8$. The errors for
positive and negative quantum states are represented by lines with squares and
$\star$, respectively, while the errors for positive and negative quantum
states by SVM are represented by lines with diamond and $+$, respectively.

The errors for positive and negative classes are also reduced in most cases,
i.e. the errors for either the positive or negative class by SVM are also
larger than the corresponding ones by S4VM for $l=10$, $30$, and $50$,
respectively. In general the errors for the negative class are smaller than
those for the positive class, maybe due to the errors of SDP. When $m=2$ and
$m=8$ the errors by S4VM for the negative class are approximately zero in most
cases.

In Table 2 (Table 3), we also list the maximum differences between errors of
positive (negative) states from S4VM and SVM by using 10, 30, or 50 labeled
quantum states and 4000 unlabeled quantum states for $m=2,4,6,$ and 8.

Table 2: The maximum differences between the errors of positive states from SVM and S4VM by using ten different $l$ labeled quantum states and 4000 unlabeled quantum states for $l=10,30,$ and $50$ and $m=2,4,6,$ and $8$. $\Delta E^{+}_{\rm{max,l}}$ is the maximum difference, errors of positive states from SVM minus the errors from S4VM, by using $l$ labeled quantum states. $m$ | $\Delta E^{+}_{\rm{max,10}}$ | $\Delta E^{+}_{\rm{max,30}}$ | $\Delta E^{+}_{\rm{max,50}}$
---|---|---|---
2 | 0.337 | 0.14 | 0.11
4 | 0.24 | 0.13 | 0.11
6 | 0.43 | 0.19 | 0.11
8 | 0.04 | 0.007 | 0.004
Table 3: The maximum differences between the errors of negative states from SVM and S4VM by using ten different $l$ labeled quantum states and 4000 unlabeled quantum states for $l=10,30,$ and $50$ and $m=2,4,6,$ and $8$. $\Delta E^{-}_{\rm{max,l}}$ is the maximum difference, the errors of negative states from SVM minus the errors from S4VM, by using $l$ labeled quantum states. $m$ | $\Delta E^{-}_{\rm{max,10}}$ | $\Delta E^{-}_{\rm{max,30}}$ | $\Delta E^{-}_{\rm{max,50}}$
---|---|---|---
2 | 0.34 | 0.10 | 0.03
4 | 0.47 | 0.12 | 0.07
6 | 0.15 | 0.10 | 0.07
8 | 0.029 | 0.032 | 0.031

Compared with the results from SVMs, semi-supervised SVMs can improve the
accuracies in most cases. The differences between average errors from S4VM and
SVM by using ten different sets of labeled quantum states are shown in Fig. 9.
All the average errors can be improved. The smaller the $l$ (the number of
labeled samples), the better the improved accuracies.

Figure 9: The differences between the average errors from S4VM and SVM by
using ten different $l$ labeled quantum states and $u$ unlabeled quantum
states for $m=2,4,6,$ and $8$. The differences between the average errors for
$l=50,30,$ and $10$ are represented by the lines with $\diamond$, $\star$ and
circles, respectively.

To show the validity of cross validation, the relationship between the cross
rate and the accuracies are shown in FIG. 10 by taking 10 labeled quantum
states. The vertical coordinates represent the accuracies of the unlabeled
states, while the horizontal coordinates represent the the accuracies of cross
validation.

Figure 10: The relationship between the cross validation accuracies and the
accuracies of unlabeled quantum states by using ten labeled quantum states and
$u$ unlabeled quantum states for $m=2,4,6,$ and $8$.

To investigate the validity of the S4VM, we also study the classification of
the generalized Werner states,

$\rho_{w}=p|\psi\rangle\langle\psi|+(1-p)\rho_{A}\otimes\frac{\rm{I}_{2}}{2},$
(13)

where $|\psi\rangle=\cos\xi|00\rangle+\sin\xi|11\rangle,$
$\rho_{A}=\rm{tr}_{B}(|\psi\rangle\langle\psi|).$ The state is unsteerable
from Alice to Bob if

$\cos^{2}2\xi\geq\frac{2p-1}{(2-p)p^{3}}.$ (14)

Let the unlabeled set be of $2500$ generalized Werner states from $p=0$ to
$p=1.$ By using 30 random labeled quantum states, the classification errors of
the unlabeled generalized Werner states for $\psi=\frac{\pi}{8}$ by SVM and
S4VM are shown in Fig. 11. Different from the above examples, we take $M=1$
here to show the performance of S4VM. It can be seen that the accuracies of
the unlabeled quantum states can be improved to 0.95 or higher. The accuracies
by S4VM are 0.15 higher than those by SVM in the best case when $m=4$ and
$m=8$. The average accuracies by S4VM are 0.969 and 0.979 for $m=4$ and $m=8$,
respectively, which are again 0.076 and 0.036 higher than those by SVM.

(a) m=8

(b) m=4

Figure 11: Classification errors of SVM and S4VM with ten different sets of
$l$ labeled quantum states for $m=4,$ and $m=8,$ $l=30$, $N_{2}=2500,$ and
$M=1$. The errors of SVM for $l=30$ are represented by line with $\star$, the
errors of S4VM for $l=30$ are represented by line with square.

We only need a small portion of the labeled data to implement the semi-
supervised learning algorithm, which can indeed save a lot of time. However,
the algorithm S4VM may still cost much time since the global simulated
annealing search, the sample strategy, and the grid search approach in S4VM
cost time. Besides, the accuracies by S4VM cannot be improved significantly
when the number of the labeled samples becomes larger, while the accuracies by
SVM are greater than 0.95.

##  V Conclusion

Semi-supervised SVMs can be used in the situation that the labeled samples are
scarce or very difficult to obtain. We have implemented a special semi-
supervised SVM-S4VM to detect the EPR steering with a small portion of labeled
quantum states and a large portion of unlabeled quantum states. Compared with
inductive SVM, the detection errors can be significantly decreased in most
cases. Our approach makes a useful step to solve the quantum correlation
detection problems based on the semi-supervised machine learning method, as it
costs time to get enough labeled quantum states. This approach may be applied
similarly to detect other quantum correlations such as quantum entanglement
and bell non-locality. Our results may also highlight the applications of
other more efficient semi-supervised machine-learning methods such as
artificial neural networks to quantum correlation detections.

ACKNOWLEDGEMENTS This work is supported by the National Natural Science
Foundation of China (NSFC) under Grant Nos. 11571313, 12071179, 12075159 and
12171044; Beijing Natural Science Foundation (Grant No. Z190005); Academy for
Multidisciplinary Studies, Capital Normal University; the Academician
Innovation Platform of Hainan Province; and Shenzhen Institute for Quantum
Science and Engineering, Southern University of Science and Technology (No.
SIQSE202001).

## References

  * (1) H. M. Wiseman, S.J. Jones , A. C. Doherty, Steering, entanglement, nonlocality, and the Einstein-Podolsky-Rosen paradox. _Phys. Rev. Lett._ , 98:140402(2007).
  * (2) E. Schrodinger, Discussion of probability relations between separated systems. _Math. Proc. Camb. Philos. Soc_. 31:555(1935).
  * (3) A. Einstein, B. Podolsky, N. Rosen, Can quantum-mechanical description of physical reality be considered complete? _Phys. Rev._ 47:777(1935).
  * (4) J. Bowles, T. Vertesi, M. T. Quintino, N. Brunner: One-way Einstein-Podolsky-Rosen steering. _Phys. Rev. Lett._ 112, 200402 (2014).
  * (5) J. L. Chen, X. J. Ye, C. F. Wu, H. Y. Su, A. Cabello, L. C. Kwek, C. H. Oh, All-versus-nothing proof of Einstein-Podolsky-Rosen steering, _Sci. Rep._ 3,2143(2013).
  * (6) J. Bowles, F. Hirsch, M. T. Quintino, N. Brunner: Sufficient criterion for guaranteeing that a two-qubit state is unsteerable. _Phys. Rev. A_ 93, 022121 (2016).
  * (7) C. Branciard, E. G. Cavalcanti, S. P. Walborn, V. Scarani, and H. M. Wiseman, One-sided device-independent quantum key distribution: Security, feasibility, and the connection with steering, _Phys. Rev. A_ 85, 010301(R)(2012).
  * (8) E. Passaro, D. Cavalcanti, P. Skrzypczyk, A. Ac$\acute{i}$n, Optimal randomness certification in the quantum steering and prepare-and-measure scenarios, _New J. Phys._ 17,113010(2015).
  * (9) B. Coyle, M. J. Hoban, E. Kashefi, One-sided device-independent certification of unbounded random numbers, _EPTCS_ 273, 14-26(2018)
  * (10) M. Piani, J. Watrous, Necessary and sufficient quantum information characterization of Einstein-Podolsky-Rosen steering, _Phys. Rev. Lett._ 114, 060404(2015).
  * (11) K. Sun, X.-J. Ye, Y. Xiao, X.-Y. Xu, Y.-C. Wu, J.-S. Xu, J.-L. Chen, C.-F. Li, and G.-C. Guo, Demonstration of Einstein-Podolsky-Rosen steering with enhanced subchannel discrimination, _npj Quantum Inf._ 4, 12(2018).
  * (12) Y. Xiang, I. Kogias, G. Adesso, and Q. He, Multipartite gaussian steering: monogamy constraints and quantum cryptography applications,_Phys. Rev. A_ 95, 010101(R)(2017).
  * (13) Q. He, L. Rosales-Z$\acute{a}$rate, G. Adesso, M. D. Reid, Secure continuous variable teleportation and Einstein-Podolsky-Rosen steering,_Phys. Rev. Lett._ 115, 180502(2015).
  * (14) H. Y. Ku, S. L. Chen, H.B Chen, N. Lambert, Y. N. Chen, F. Nori, Temporal steering in four dimensions with applications to coupled qubits and magnetoreception. _Phys. Rev. A_ 94, 062126(2016).
  * (15) C. Y. Chiu, N. Lambert, T.L. Liao, F. Nori, C. M. Li, No-cloning of quantum steering. _NPJ Quantum Inf._ , 2, 16020(2016).
  * (16) S. L. Chen, N. Lambert, C. M. Li, G. Y. Chen, Y. N. Chen, A. Miranowicz, F. Nori, Spatio-temporal steering for testing nonclassical correlations in quantum networks. _Sci. Rep._ , 7, 3728(2017).
  * (17) E. G. Cavalcanti, S. J. Jones, H. M. Wiseman, M. D. Reid, Experimental criteria for steering and the Einstein-Podolsky-Rosen paradox. _Phys. Rev. A_ 2009, 80, 032112.
  * (18) M. F. Pusey, Negativity and steering: A stronger Peres conjecture. _Phys. Rev. A_ 2013, 88, 032313.
  * (19) C .L. Ren, H. Y. Su, H. F. Shi, J. L. Chen, Maximally steerable mixed state based on the linear steering inequality and the Clauser-Horne-Shimony-Holt-like steering inequality. _Phys. Rev. A_ 2018, 97, 032119.
  * (20) Y. N. Chen, C. M. Li, N. Lambert, S. L. Chen, Y. Ota, G. Y. Chen, F. Nori, Temporal steering inequality. _Phys. Rev. A_ 2014, 89, 032112.
  * (21) M. Zukowski, A. Dutta, Z. Yin, Geometric Bell-like inequalities for steering. _Phys. Rev. A_ 2015, 91, 032107.
  * (22) H. Zhu, M. Hayashi, L. Chen, Universal steering inequalities. _Phys. Rev. Lett._ 2016, 116, 070403.
  * (23) S. P. Walborn, A. Salles, R. M. Gomes, F. Toscano, P. H. Souto Ribeiro, Revealing hidden Einstein-Podolsky-Rosen nonlocality. _Phys. Rev. Lett._ 2011, 106, 130402.
  * (24) J. Schneeloch, C. J. Broadbent, S. P. Walborn, E. G. Cavalcanti, J.C. Howell, Einstein-Podolsky-Rosen steering inequalities from entropic uncertainty relations. _Phys. Rev. A_ 2013, 87, 062103.
  * (25) T. Pramanik, M. Kaplan, A. S. Majumdar, Fine-grained Einstein-Podolsky-Rosen steering inequalities. _Phys. Rev. A_ 2014, 90, 050305.
  * (26) I. Kogias, P. Skrzypczyk, D. Cavalcanti, A. Acin, and G. Adesso, Hierarchy of Steering Criteria Based on Moments for All Bipartite Quantum Systems 2015, _Phys. Rev. Lett._ 115, 210401.
  * (27) D. Cavalcanti, P. Skrzypczyk, Quantum steering: a review with focus on semidefinite programming, _Rep. Prog. Phys._ 80 024001 (2017).
  * (28) Y. C. Ma, M. H. Yung, Transforming bell’s inequalities into state classifiers with machine learning, _npj, Quantum. Info.,_ 4,34(2018).
  * (29) S. R. Liu, S. L. Huang, K. R. Li, et al, A separability-entanglement classifier via machine learning, _Phys. Rev. A_ 98, 012315 (2018).
  * (30) D. L. Deng, X. Li and S. D. Sarma, Quantum entanglement in neural network states, _Phys. Rev. X_ 7, 021021 (2017).
  * (31) M. Yang, C. L. Ren, Y. C. Ma, Y. Xiao, X. J. Ye, L. L. Song, J. S. Xu, M. H. Yung, C. F. Li, G. C. Guo, Experimental Simultaneous Learning of Multiple Non-Classical Correlations,_Phys. Rev. Lett._ 123, 190401 (2019).
  * (32) N. Wiebe, C. Granade, C. Ferrie, and D. G. Cory, Hamiltonian learning and certification using quantum resources, _Phys. Rev. Lett._ 112, 190501 (2014).
  * (33) S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, A structural approach to relaxation in glassy liquids, _Nat. Phys._ 12, 469 (2016).
  * (34) C. L. Ren,C. B. Chen, Steerability detection of arbitrary 2-qubit state via machine learning, _Phys. Rev. A_ 100, 022314 (2019).
  * (35) Y.Q.Zhang, L.J. Yang, Q.L. He, L. Chen, Machine learning on quantifying quantum steerability, _Quan. inf. Process._ , 19(8):263(2020)
  * (36) X. Zhu, AB. Goldberg, Introduction to semi-supervised learning. _Synth. Lect. Artif. Intell. Mach. Learn._ , 3(1):1-130(2009).
  * (37) S.F. Ding, Z.B. Zhu, X. K. Zhang, An overview on semi-supervised support vector machine, _Neural. Comput. Applic,_ 28(5):969-978(2017).
  * (38) Y.F. Li, Z.H. Zhou, Towards making unlabeled data never hurt, _IEEE Transactions on Pattern Analysis and Machine Intelligence_ ,37(1):175-188(2015).
  * (39) T. Joachims, Transductive inference for text classification using support vector machines. _In Proc. 16th Int. Conf. On Mach. Learn., Bled, Slovenia,_ 200-209(1999).
  * (40) L.M. Yang, Q.Y. Su, B.Y. Yang, D. Tong, and X. Xiao, A new semi-supervised support vector machine classifier based on wavelet transform and its application in the iris image recognition.,_Int J Appl Math Stat_ 52(5):86-93(2014).
  * (41) K. Lu, X. He, J. Zhao, Semi-supervised support vector learning for face recognition. _Advances in neural networks-ISNN_ , Springer, Berlin, 104-109(2006).
  * (42) J. C. Ang, H. Haron, H. N. A. Hamed , Semi-supervised SVM-based feature selection for cancer classification using microarray gene expression data. _In M. Ali., Y. Kwon, CH. Lee , J. Kim, K. Kim (eds) Current Approaches in Applied Artificial Intelligence IEA/AIE 2015, Seoul, South Korea,Lecture Notes in Artificial Intelligence-ISNN_ , Springer, Cham., 9101,468-477(2015).
  * (43) C. Cortes, V. Vapnik, Support-vector networks, _Mach. Learn._ , 20(3), 273-297(1995).
  * (44) S. Abe, Support Vector Machines for Pattern Classification. Berlin, Germany: Springer-Verlag(2005).
  * (45) O. Chapelle, A. Zien, Semi-supervised learning by low density separation. In Proc. of the 10th Int. Wksh. on Arti. Intell. and Stat., Savannah,U.S.A., 57-64(2005).
  * (46) X. Zhu, Semi-Supervised Learning Literature Survey, _technical report_ , Univ. of Wisconsin-Madison, 2007


# Time-optimal control of two-level quantum systems by piecewise constant
pulses

E. Dionis111Laboratoire Interdisciplinaire Carnot de Bourgogne (ICB), UMR 6303
CNRS-Université Bourgogne, 9 Av. A. Savary, BP 47 870, F-21078 Dijon Cedex,
France, D. Sugny222Laboratoire Interdisciplinaire Carnot de Bourgogne (ICB),
UMR 6303 CNRS-Université Bourgogne, 9 Av. A. Savary, BP 47 870, F-21078 Dijon
Cedex, France<EMAIL_ADDRESS>
###### Abstract

We apply an extension of the Pontryagin Maximum Principle to derive time-
optimal controls of two-level quantum systems by means of piecewise constant
pulses. Global optimal solutions are obtained for state-to-state transfer in
the cases with one and two controls. Exact quantum speed limits are
established as a function of the sampling period. We observe numerically an
exponential convergence towards the minimum time in the continuous limit when
this period goes to zero. We show that this convergence is only polynomial for
a linearized quantum system. We discuss the experimental impact of this
result.

##  1 Introduction

Designing fast and accurate control sequences while taking into account
specific experimental constraints is a central task in quantum control [1, 2,
3, 4, 5, 6]. In this framework, optimal control theory (OCT) [8, 7] is a
powerful and versatile tool in order to propose efficient answers to the
different issues that can be encountered in experimental setups [1, 3]. For
complex quantum systems, optimal controls are computed by using iterative
optimization algorithms [9, 11, 10, 12] such as the well-known GRAPE procedure
[9].

A key limitation of pulse-shaping devices is their ability to generate time-
continuous controls. They are therefore usually replaced by piecewise constant
pulses in experiments with a minimum sampling period over which the controls
remain constant. This hybrid situation is characterized by a digital control
which evolves in a discrete way, while the state of the system varies
continuously. This constraint was one of the motivations for the development
of GRAPE, a gradient-based algorithm directly optimizing piecewise constant
functions that can be implemented experimentally. Even if the efficiency of
numerical algorithms has been demonstrated in a variety of domains [1, 3],
their application is not completely satisfactory because only local optimal
protocols are achieved, therefore with no certitude about the global
optimality of the control process. This problem can be solved from the
Pontryagin Maximum Principle (PMP), which has the decisive property of
transforming the initial infinite-dimensional control problem into a
generalized Hamiltonian system subject to a maximization condition and
boundary constraints. The optimal solution corresponds then to a Hamiltonian
trajectory both reaching the target state and minimizing a cost functional
such as the control time. The low and finite dimension of this new dynamic
space simplifies greatly the search for globally optimal controls [7]. This
method has been applied recently to a series of quantum control problems [14,
13, 18, 15, 17, 16, 19, 20]. However, as with the original formulation of
Pontryagin [8], such studies consider that the control can be modified at any
time, leading to continuous functions or possibly having discontinuities in
the case of bang-bang protocols [7]. As this assumption is too strong for
standard wave-form generators, the PMP only gives mathematical optima which
cannot necessarily be physically achieved in practice. An interesting example
is given by time-optimal control processes in which the goal is to perform a
given task as quickly as possible. In this framework, quantum speed limits
provide lower bounds on the minimum time required for a given control protocol
[21, 22, 23]. Limiting the shaping to piecewise constant controls necessarily
increases the minimum time. A key point for the experimental performance of
the process is the quantitative impact of the control digitalization on this
time.

In this study, we aim to take a step towards solving this control problem by
applying a recent mathematical extension of the PMP to piecewise constant
pulses [25, 26, 24]. In this new and rigorous formulation (called _discrete_
below to differentiate it from the standard _continuous_ one), the overall
structure of the PMP is not modified. The main difference between the
continuous and the discrete versions relies on the maximization condition in
which the instantaneous condition is replaced by an integral one. This
integral condition can be interpreted as an average over the sampling period
of the standard maximization condition. We apply this general procedure to the
control of two-level quantum systems. We consider two benchmark control
protocols with respectively one and two control parameters and for which the
time-optimal solutions are bang-bang and continuous. We find in the two cases
the optimal dynamics for any sampling period and we show that the
corresponding minimum time converges exponentially toward its continuous limit
when this period goes to zero. We show that this convergence is only
polynomial for a specific linearized quantum system. Finally, we describe the
link between this new version of the PMP and GRAPE and we propose a new
formulation of the algorithm based on the PMP with an exact gradient. A
comparison with different versions of GRAPE is performed for a specific
control problem. We consider for GRAPE both the split-operator [9] and the
auxiliary matrix [28] approaches. In the first formulation, the gradient of
the fidelity with respect to the control is approximated by replacing the
elementary evolution operator over a time step into a product of two or more
factors, while in the second the gradient is exactly computed by extending the
size of the Hilbert space and computing the matrix exponential in this new
space. We show the superiority of the PMP-based formulation for the example
under study.

The paper is organized as follows. Section 2 introduces the PMP for piecewise
constant pulses. The time-optimal solutions for controlling two-level quantum
systems with two and one control parameters are respectively presented in Sec.
3 and 4. The optimal solution for the control of a linear system is described
in Sec. 5. Section 6 focuses on the link between this new formulation of the
PMP and gradient-based optimal algorithms. Conclusion and prospective views
are given in Sec. 7. Technical computations are reported in Appendix A.

##  2 Methodology

We briefly recall in this section the main steps of the application of the PMP
for bilinear differential systems. We refer the interested reader to a recent
paper about the PMP in quantum control for details [7]. We describe the
differences in the PMP between the continuous and the discrete cases. We
present a numerical method called the shooting algorithm to solve the
optimization equations for low-dimensional dynamical systems.

We consider the control of a system described by the state $\displaystyle
X\in\mathbb{R}^{n}$, $\displaystyle n\geq 1$, whose dynamics are described by
the following differential system

$\displaystyle\dot{X}=H_{0}X+u_{1}H_{1}X+u_{2}H_{2}X$

where $\displaystyle H_{0}$, $\displaystyle H_{1}$ and $\displaystyle H_{2}$
are $\displaystyle n\times n$ skew-symmetric real matrices, and $\displaystyle
u_{1}(t)$, $\displaystyle u_{2}(t)$ the two real control parameters. Note that
the approach detailed below can be straightforwardly extended to any number of
control parameters. The goal of the control process is to bring in minimum
time $\displaystyle t_{f}$ the state of the system from $\displaystyle X_{0}$
to the target state $\displaystyle X_{f}$ with the constraint $\displaystyle
u_{1}(t)^{2}+u_{2}(t)^{2}\leq 1$ at any time $\displaystyle t\in[0,t_{f}]$.

We first state the PMP in the standard continuous situation where
$\displaystyle u_{1}$ and $\displaystyle u_{2}$ are assumed to be sufficiently
smooth functions [7]. We only consider the regular case in which
$\displaystyle u_{1}$ and $\displaystyle u_{2}$ are not simultaneously equal
to 0, except in isolated points. The Pontryagin Hamiltonian $\displaystyle
H_{P}$ can be expressed as

$\displaystyle
H_{P}(X,P,u_{1},u_{2})=P^{\intercal}[H_{0}X+u_{1}H_{1}X+u_{2}H_{2}X]$

with the adjoint state $\displaystyle P\in\mathbb{R}^{n}$. Since the final
time $\displaystyle t_{f}$ is free, $\displaystyle H_{P}$ is a constant fixed
to 1 (see [7] for technical details). The adjoint state $\displaystyle P$
fulfills the same equation as $\displaystyle X$

$\displaystyle\dot{P}=H_{0}P+u_{1}H_{1}P+u_{2}H_{2}P$

The candidate controls $\displaystyle u_{1}^{*}$ and $\displaystyle u_{2}^{*}$
for optimality satisfy the following maximization condition at any time
$\displaystyle t$

$\displaystyle
H(X,P,u_{1}^{*}(X,P),u_{2}^{*}(X,P))=\max_{u_{1}^{2}+u_{2}^{2}\leq
1}H_{P}(X,P,u_{1},u_{2})$

A straightforward computation shows that the solutions can be written as

$\displaystyle\displaystyle
u_{1}^{*}=\frac{P^{\intercal}H_{1}X}{\sqrt{(P^{\intercal}H_{1}X)^{2}+(P^{\intercal}H_{2}X)^{2}}},$
$\displaystyle\displaystyle
u_{2}^{*}=\frac{P^{\intercal}H_{2}X}{\sqrt{(P^{\intercal}H_{1}X)^{2}+(P^{\intercal}H_{2}X)^{2}}}.$

Hamiltonian trajectories of $\displaystyle H$ can then be computed from the
expressions of $\displaystyle u_{1}^{*}$ and $\displaystyle u_{2}^{*}$. The
shooting approach consists in finding the initial adjoint state $\displaystyle
P(0)$ such that the corresponding trajectory reaches exactly the target state
$\displaystyle X(t_{f})=X_{f}$ in minimum time. The optimal controls are then
deduced from the time evolution of $\displaystyle X$ and $\displaystyle P$. If
the optimization equations cannot be integrated analytically, a Newton
algorithm is generally used to estimate the initial adjoint state with a very
high numerical precision.

We study now the case where $\displaystyle u_{1}$ and $\displaystyle u_{2}$
are piecewise constant pulses with the same constraint $\displaystyle
u_{1}(t)^{2}+u_{2}(t)^{2}\leq 1$. We denote by $\displaystyle T$ the sampling
period. The controls are characterized by $\displaystyle N$ values
$\displaystyle u_{1}^{(k)}$, $\displaystyle u_{2}^{(k)}$ in the time intervals
$\displaystyle[kT,(k+1)T]$, $\displaystyle k\in\\{0,1,2,\cdots,(N-1)\\}$, with
$\displaystyle(N-1)T+\delta T=t_{f}$ and $\displaystyle 0<\delta T\leq T$.
Note that the last time interval of length $\displaystyle\delta T$ can be
lower or equal to $\displaystyle T$. This point is discussed in the different
examples. We consider here for simplicity a regular time discretization but
the same approach can be used if the sampling period is not constant.

The application of the PMP is very similar to the continuous version [25, 26].
The Pontryagin Hamiltonian is defined the same way (except that $\displaystyle
H_{P}$ is no more constant in time but still equal to +1 at the final time
$\displaystyle t_{f}$) and the state $\displaystyle X$ and adjoint state
$\displaystyle P$ satisfy the same differential equations. The only difference
corresponds to the formulation of the maximization condition. Introducing the
integrals

$\begin{split}&H_{1}^{(k)}=\int_{kT}^{(k+1)T}P(\tau)^{\intercal}H_{1}X(\tau)d\tau\\\
&H_{2}^{(k)}=\int_{kT}^{(k+1)T}P(\tau)^{\intercal}H_{2}X(\tau)d\tau\end{split}$
(1)

the maximization condition can be expressed for all $\displaystyle k$ as

$\displaystyle
H_{1}^{(k)}(v_{1}-u_{1}^{(k)})+H_{2}^{(k)}(v_{2}-u_{2}^{(k)})\leq 0$

for any couple $\displaystyle(v_{1},v_{2})$ such that $\displaystyle
v_{1}^{2}+v_{2}^{2}\leq 1$. We then deduce that the extremal controls are
given by

$\begin{split}&u_{1}^{(k)*}=\frac{H_{1}^{(k)}}{\sqrt{(H_{1}^{(k)})^{2}+(H_{2}^{(k)})^{2}}}\\\
&u_{2}^{(k)*}=\frac{H_{2}^{(k)}}{\sqrt{(H_{1}^{(k)})^{2}+(H_{2}^{(k)})^{2}}}.\end{split}$
(2)

If the last time interval is not of length $\displaystyle T$ then the
integrals (1) are taken from $\displaystyle(N-1)T$ to $\displaystyle t_{f}$.
From this modified maximization condition, the same shooting approach can be
used to find the initial value of $\displaystyle P$ bringing the system to the
target state. Note that $\displaystyle\delta T$ is also computed by the
shooting procedure when this parameter is not fixed. Finally, we point out
that in the limit $\displaystyle T\to 0$, the piecewise constant controls
converge towards their continuous limit given by the standard PMP [25, 26].

##  3 The case of two controls

###  3.1 The general formalism

As a first application, we consider the time optimal control of a two-level
quantum system by means of two resonant controls [7]. In a given rotating
frame, the dynamics of the system is given by the Schrödinger equation written
in units such that $\displaystyle\hbar=1$

$\displaystyle
i|\dot{\psi}\rangle=H|\psi\rangle=\frac{1}{2}(\omega_{x}\sigma_{x}+\omega_{y}\sigma_{y})|\psi\rangle$

where $\displaystyle\omega_{x}$ and $\displaystyle\omega_{y}$ are the two real
control amplitudes along the $\displaystyle x$ and $\displaystyle y$\-
directions with the constraint $\displaystyle\omega_{x}^{2}+\omega_{y}^{2}\leq
1$, and $\displaystyle\sigma_{x,y}$ the two Pauli matrices. We limit our study
to state-to-state transfer. In the continuous case, it has been shown that the
optimal protocol consists in saturating the bound at any time, that is
$\displaystyle\omega_{x}^{2}+\omega_{y}^{2}=1$. This saturation can be
qualitatively explained by the fact that a fastest control protocol requires
generally maximum control intensity. The same assumption is made in this
section for the discrete formulation of the PMP. We introduce to this aim the
Bloch vector $\displaystyle X=(x,y,z)$ which satisfies

$\displaystyle\dot{X}=(\omega_{x}M_{x}+\omega_{y}M_{y})X$

with the matrices $\displaystyle M_{x}$, $\displaystyle M_{y}$ and
$\displaystyle M_{z}$ defined as

$\displaystyle M_{x}=\begin{pmatrix}0&0&0\\\ 0&0&-1\\\
0&1&0\end{pmatrix},~{}M_{y}=\begin{pmatrix}0&0&1\\\ 0&0&0\\\
-1&0&0\end{pmatrix},~{}M_{z}=\begin{pmatrix}0&-1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}$

The Pontryagin Hamiltonian $\displaystyle H_{P}$ can be expressed as

$\displaystyle H_{P}=P^{\intercal}(\omega_{x}M_{x}+\omega_{y}M_{y})X$

where the adjoint state $\displaystyle P=(p_{x},p_{y},p_{z})$ fulfills the
differential equation

$\displaystyle\dot{P}=(\omega_{x}M_{x}+\omega_{y}M_{y})P.$

We consider piecewise constant pulses described by a phase
$\displaystyle\varphi$ such that $\displaystyle\omega_{x}=\cos\varphi$ and
$\displaystyle\omega_{y}=\sin\varphi$. We denote $\displaystyle\varphi_{k}$
the amplitude of the control parameter in the interval
$\displaystyle[kT,(k+1)T]$, $\displaystyle 0\leq k\leq(N-1)$. We also
introduce $\displaystyle X_{k}$ and $\displaystyle P_{k}$ the values of the
state and adjoint state at time $\displaystyle t=kT$, and the quantities
$\displaystyle I_{x,y,z}^{(k)}=P_{k}M_{x,y,z}X_{k}$. Using Eq. (2), the
maximization condition of $\displaystyle H_{P}$ leads to the condition

$\tan(\varphi_{k})=\frac{\int_{kT}^{(k+1)T}P_{k}^{\intercal}(\tau)M_{y}X_{k}(\tau)d\tau}{\int_{kT}^{(k+1)T}P_{k}^{\intercal}(\tau)M_{x}X_{k}(\tau)d\tau}$
(3)

with $\displaystyle X_{k}(\tau)=R(\tau)X_{k}$ and $\displaystyle
P_{k}(\tau)=R(\tau)P_{k}$, $\displaystyle\tau\in[kT,(k+1)T]$, and
$\displaystyle
R(\tau)=\exp[(\cos(\varphi_{k})M_{x}+\sin(\varphi_{k})M_{y})\tau]$. The
rotation matrix $\displaystyle R(\tau)$ can be explicitly computed as

$\displaystyle
R(\tau)=\begin{pmatrix}\sin^{2}(\varphi_{k})\cos(\tau)+\cos^{2}(\varphi_{k})&(1-\cos(\tau))\sin(\varphi_{k})\cos(\varphi_{k})&\sin(\tau)\sin(\varphi_{k})\\\
(1-\cos(\tau))\sin(\varphi_{k})\cos(\varphi_{k})&\sin^{2}(\varphi_{k})+\cos(\tau)\cos^{2}(\varphi_{k})&-\sin(\tau)\cos(\varphi_{k})\\\
-\sin(\tau)\sin(\varphi_{k})&\sin(\tau)\cos(\varphi_{k})&\cos(\tau)\end{pmatrix}$

Setting $\displaystyle h_{k}=\tan(\varphi_{k}/2)$ and computing the integrals
in Eq. (3), we arrive at a polynomial of order 2 in $\displaystyle h_{k}$

$\displaystyle\displaystyle[I_{y}^{(k)}\sin(T)+(1-\cos(T)I_{z}^{(k)})]h_{k}^{2}+2I_{x}^{(k)}\sin(T)h_{k}$
$\displaystyle\displaystyle+[(1-\cos(T)I_{z})-I_{y}\sin(T)]=0.$

The two real solutions (if they exist) are candidates to be the optimal value
of $\displaystyle\varphi_{k}$. As in the continuous case, they lead to two
symmetric trajectories with respect to the equator. We select in Sec. 3.2 the
positive root which corresponds to $\displaystyle p_{z}(0)<0$ and to a dynamic
in the southern hemisphere. Finally, we point out that in the continuous limit
$\displaystyle T\to 0$, we recover the extremal solution
$\displaystyle\tan\varphi_{k}=\frac{I_{y}}{I_{x}}$.

###  3.2 Numerical results

We apply the shooting method to a state-to-state transfer on the Bloch sphere.
The initial and target states are respectively of coordinates
$\displaystyle(1,0,0)$ and $\displaystyle(0,1,0)$. As mentioned in Sec. 2, we
consider two cases for which $\displaystyle t_{f}$ is either equal to
$\displaystyle NT$ or to $\displaystyle(N-1)T+\delta T$. In the first
situation, for a fixed time step value, $\displaystyle N$, we only need to
find the initial adjoint state $\displaystyle P(0)$ and the sampling period
$\displaystyle T$ such that the trajectory reaches the target state in a
minimum time $\displaystyle t_{f}$. In the numerical simulations, we therefore
optimize four parameters, $\displaystyle p_{x}(0)$, $\displaystyle p_{y}(0)$,
$\displaystyle p_{z}(0)$ and $\displaystyle T$ such that $\displaystyle
x(t_{f}=NT)=0$, $\displaystyle y(t_{f})=1$, $\displaystyle z(t_{f})=0$ and
$\displaystyle H_{P}(t_{f})=1$. In the second case, $\displaystyle T$ is fixed
and we find the coordinates of $\displaystyle P(0)$, $\displaystyle N$ and
$\displaystyle t_{f}$ such that the four final conditions are met. The
parameter $\displaystyle\delta T$ is then given by $\displaystyle\delta
T=t_{f}-(N-1)T$. We point out that we use exactly $\displaystyle N$ time steps
in the first approach, while the last step may be different from the others in
the second case.

The continuous limit can be integrated exactly and leads to the minimum time
$\displaystyle t_{f}^{(c)}=\frac{\pi\sqrt{3}}{2}\simeq 2.7207$ (see Ref. [7]
for the derivation). The corresponding initial adjoint state $\displaystyle
P(0)$ has the coordinates $\displaystyle(p_{x}(0),\frac{1}{\sqrt{3}},\pm 1)$
where the first coordinate is not fixed. We stress that the same control
process is obtained for any value of $\displaystyle p_{x}(0)$. Note that a
similar behavior is observed in the discrete case. The two possible values of
$\displaystyle p_{z}(0)$ correspond to two symmetric trajectories with respect
to the equator on the Bloch sphere [7]. Details about the results of the
shooting method and a comparison with the continuous limit are given in
Appendix A.

Figure 1 compares the results obtained by the two different approaches. We
observe a rapid convergence of the minimum time towards the continuous limit
since the ratio between the two times is respectively of the order of
$\displaystyle 10^{-3}$ and $\displaystyle 10^{-5}$ when $\displaystyle N=10$
and 100. The results of Fig. 1 clearly show that this convergence is
exponential. As could be expected, releasing the constraint $\displaystyle
t_{f}=NT$ allows to decrease the minimum time. In this second case, the curve
$\displaystyle t_{f}$ as a function of $\displaystyle T$ has a quasi-quadratic
shape when $\displaystyle T$ is not too small and horizontal tangents when
$\displaystyle\delta T=T$.

Figure 1: Evolution of the minimum time $\displaystyle t_{f}$ as a function of
the number of time steps $\displaystyle N$ (a) (with $\displaystyle\delta
T=T$) and of the sampling period $\displaystyle T$ (b) when
$\displaystyle\delta T$ is free. In panel (a), the minimum value of
$\displaystyle N$ is 2. The points are computed for integer values of
$\displaystyle N$. The solid line is just to guide the reading. Quantities
plotted are dimensionless.

An example of optimal solutions is given in Fig. 2 for $\displaystyle N=3$ and
$\displaystyle\delta T\neq T$. We observe that the discrete and the continuous
trajectories are very close to each other on the Bloch sphere while the two
controls are quite different even if they have a similar evolution.

Figure 2: Panel (a) displays the optimal trajectories for $\displaystyle N=3$
and $\displaystyle\delta T\neq T$ (solid lines) and in the continuous limit
(dashed lines) for the coordinates $\displaystyle x$ (black), $\displaystyle
y$ (red or light gray) and $\displaystyle z$ (blue or solid gray). The black
vertical line corresponds to the minimum continuous time $\displaystyle
t_{f}^{(c)}$. The corresponding controls are represented in panel (b) with the
same color code. The dots indicate the position of the switching times
$\displaystyle t_{k}$. Quantities plotted are dimensionless.

As a quantitative example, we briefly analyze the control of a spin 1/2
particle in Nuclear Magnetic Resonance [39]. We consider the experimental
parameters given in [48], i.e. a time digitization of 0.5 $\displaystyle\mu$s
and a maximum amplitude of $\displaystyle\nu=100$ kHz. The experimental time
$\displaystyle t_{\textrm{exp}}$ is then given by $\displaystyle t=2\pi\nu
t_{\textrm{exp}}$ where $\displaystyle t$ is the time in normalized unit used
in the theoretical part of Sec. 3. For $\displaystyle N=9$ and
$\displaystyle\delta T\neq T$, the minimum time is equal to 4.34
$\displaystyle\mu$s, while this time is 4.33 $\displaystyle\mu$s in the
continuous limit. Due to the exponential convergence of the optimization
procedure, we observe that the difference is very small between the two times
even for relatively low values of $\displaystyle N$.

##  4 The case of one control

###  4.1 Description of the optimal control problem

We consider in this section a similar control problem for a two-level quantum
system, but with only one control $\displaystyle\omega(t)$ such that
$\displaystyle|\omega(t)|\leq 1$ [14, 7]. The dynamics of the Bloch vector are
governed by the following differential equation

$\dot{X}=(\Delta M_{z}+\omega M_{x})X$ (4)

where the offset $\displaystyle\Delta$ is a constant parameter and we
introduce $\displaystyle\Omega=\sqrt{1+\Delta^{2}}$. The goal of the control
is to steer the system from the north pole to the south one in minimum time.
This corresponds to a complete population transfer from the ground state to
the excited state in the two-level quantum system. The adjoint state
$\displaystyle P$ fulfills the same equation as $\displaystyle X$ and the
Pontryagin Hamiltonian $\displaystyle H_{P}$ is given by

$\displaystyle H_{P}=\Delta(p_{y}x-p_{x}y)+\omega(p_{z}y-p_{y}z)$

The continuous limit can be solved exactly with the PMP [7, 14]. In the case
$\displaystyle|\Delta|\leq 1$, the optimal solution is a bang-bang control
with a switching time equal to $\displaystyle
t_{1}=\frac{1}{\Omega}(\pi-\arccos(\Delta^{2}))$ or $\displaystyle
t_{2}=\frac{1}{\Omega}(\pi+\arccos(\Delta^{2}))$. These two symmetric
solutions lead to the same total time $\displaystyle t_{f}=2\pi/\Omega$. The
discrete version of the PMP is applied in Appendix B. The numerical results
are described in Sec. 4.2.

###  4.2 Numerical results

The two shooting methods described in Sec. 2 can be applied to this control
problem. We only present the results corresponding to $\displaystyle\delta
T=T$. In this case, we observe numerically that the solution is very sensitive
to the initial conditions of the adjoint vector. The initial guess was adapted
to each value of $\displaystyle N$. Despite these numerical difficulties, we
show in Fig. 3 that the minimum time $\displaystyle t_{f}$ converges very
quickly towards its continuous limit as a function of $\displaystyle N$, the
difference being of the order of $\displaystyle 10^{-4}$ for $\displaystyle
N=20$. The cloud of points that can be seen in Fig. 3b is a signature of the
numerical instabilities of the shooting algorithm, even if a clear trend can
be observed with an approximate exponential convergence. An example of optimal
trajectories and controls is plotted in Fig. 4. The bang-bang optimal solution
in the continuous limit is replaced by another piecewise constant control but
with fixed switching times. While the former switches between -1 and 1
(singular control is not optimal in this case), the latter can take values in
the interval $\displaystyle(-1,1)$ as allowed by the discrete version of the
PMP. This degree of freedom compensates for the fact that switching cannot
take place at any time $\displaystyle t\in[0,t_{f}]$. For most values of
$\displaystyle N$, we observe that the algorithm finds the optimal control
displayed in Fig. 4 with three constant parts of respective amplitudes -1,
$\displaystyle\omega_{0}$ and +1. The control $\displaystyle\omega=\omega_{0}$
is applied during only one time step around the switching time in the
continuous limit. This parameter $\displaystyle\omega_{0}\in[-1,1]$ depends in
a rather complex way on $\displaystyle N$.

Figure 3: Evolution of the minimum time $\displaystyle t_{f}$ as a function of
the number of steps $\displaystyle N$ (crosses). The solid line is just to
guide the eye. The panel (b) corresponds to the same data as panel (a) but in
a logarithmic scale. Quantities plotted are dimensionless.

Figure 4: Panel (a): Optimal trajectories for $\displaystyle N=4$ and
$\displaystyle\delta T=T$ (solid lines) and in the continuous limit (dashed
lines) for the coordinates $\displaystyle x$ (black), $\displaystyle y$ (red
or light gray) and $\displaystyle z$ (blue or solid gray). The black vertical
line corresponds to the minimum continuous time $\displaystyle t_{f}^{(c)}$.
Panel (b): Optimal control in the discrete case. The dots indicate the
position of the times $\displaystyle t_{k}=kT$. The vertical red line
represents the position of the switching time for the bang-bang solution in
the continuous limit. Quantities plotted are dimensionless.

###  4.3 Optimal control of a Landau-Zener-type Hamiltonian

We analyze in this section a closely related control of a two-level quantum
system that has been widely investigated both theoretically [15, 43, 44] and
experimentally [40, 41, 42] in recent years. The two levels are coupled
through a Landau-Zener Hamiltonian leading to the dynamics of Eq. (4) where
$\displaystyle\omega$ is a constant parameter (set to 0.5 in the numerical
computation) describing the coupling between the two levels and
$\displaystyle\Delta(t)$ the control parameter. In the case of a Bose Einstein
Condensate (BEC) in an optical lattice [40], this latter corresponds to the
quasimomentum. The initial and final states of the control are respectively
the ground adiabatic levels for $\displaystyle\Delta=+1$ and
$\displaystyle\Delta=-1$. The two states are close to the north and south
poles of the Bloch sphere. The time optimal solution in the continuous limit
has been established in [15]. With the constraint
$\displaystyle|\Delta(t)|\leq\Delta_{\textrm{max}}$ (with
$\displaystyle\Delta_{\textrm{max}}=2$ below), the control law is the
concatenation of two bang pulses of maximum intensity at the start and end of
the sequence with a zero singular control in the middle. It can be shown that
this protocol reaches the QSL when the bang controls are replaced by Dirac
pulses. Using the equations introduced in Sec. 4.1, where this time the role
of control is played by $\displaystyle\Delta$, the same control problem can be
solved with the constraint of piecewise constant pulses. This study is
illustrated in Fig. 5 for the case $\displaystyle N=5$ and
$\displaystyle\delta T=T$. We obtain a ratio between the two control times
equal to 1.0554. We surprisingly observe that the discrete solution does not
use maximum intensity pulses for which $\displaystyle\Delta(t)=2$, but only
less strong controls. As expected, the structure of the two optimal sequences
are very close to each other.

Figure 5: Plot of the optimal controls for the Landau-Zener Hamiltonian. The
black and blue lines depict respectively the continuous and the discrete
solutions for $\displaystyle N=5$. Quantities plotted are dimensionless.

##  5 Time-optimal control in a linear system

We have observed numerically in Sec. 3 and 4 that the minimum control time by
means of piecewise constant pulses converges exponentially towards the time in
the continuous limit when the number of time steps $\displaystyle N$ tends to
infinity. The aim of this section is to investigate the same issue but for a
linearized quantum system. We show on a specific example that the
linearization process is responsible for the change in convergence which is no
more exponential, but polynomial with respect to the number of time steps.

A linearization of a quantum system can be achieved around an equilibrium
point. This corresponds to the north or south pole of the Bloch sphere for a
two-level quantum system. The control of quantum dynamics is then reduced to
that of a spring which provides information on the initial problem [46, 45,
47]. This approach is well-known in Nuclear Magnetic Resonance and usually
called the _small flip angle approximation_ [39].

As a non-trivial example, we consider a linearization around the north pole of
the quantum system investigated in Sec. 3, but with non resonant controls.
This slight modification is due to the fact that the optimal control is
constant in the resonant case. In a given rotating frame, the dynamics of the
two-level quantum system can be expressed in Bloch coordinates
$\displaystyle(x,y,z)$ as

$\begin{split}&\dot{x}=-\omega y+\sin(\varphi)z\\\ &\dot{y}=\omega
x-\cos(\varphi)z\\\ &\dot{z}=-\sin(\varphi)x+\cos(\varphi)y\end{split}$

where $\displaystyle\omega$ is the offset term and $\displaystyle\varphi(t)$
the control parameter. Starting from the point $\displaystyle(0,0,1)$, the
goal of the control is to reach a state in a neighborhood of the north pole.
Since $\displaystyle z\simeq 1$ during the control process, we make this
approximation in the dynamical system by replacing $\displaystyle z$ by 1 and
by considering only the first two equations. The linearized system on
$\displaystyle\mathbb{R}^{2}$ are then governed by the following differential
equations

$\begin{split}&\dot{x}=-\omega y+\sin\varphi\\\ &\dot{y}=\omega
x-\cos\varphi\end{split}$

In this case, the objective of the control is to steer in minimum time the
system from the origin of coordinates $\displaystyle(0,0)$ to the point of
coordinates $\displaystyle(1,0)$. We point out that the conclusion of this
section would be the same with other initial and final states, as long as the
linear approximation is valid. We first solve this problem in the continuous
limit. The Pontryagin Hamiltonian $\displaystyle H_{P}$ is given by

$\displaystyle H_{P}=\omega(p_{y}x-p_{x}y)+p_{x}\sin\varphi-p_{y}\cos\varphi$

where the adjoint state $\displaystyle P=(p_{x},p_{y})$ is solution of
$\displaystyle\dot{p}_{x}=-\omega p_{y}$ and $\displaystyle\dot{p}_{y}=\omega
p_{x}$. We get

$\begin{split}p_{x}(t)=p_{x}(0)\cos(\omega t)-p_{y}(0)\sin(\omega t)\\\
p_{y}(t)=p_{x}(0)\sin(\omega t)+p_{y}(0)\cos(\omega t)\end{split}$

Note that $\displaystyle P^{2}=p_{x}^{2}+p_{y}^{2}$ is a constant of the
motion. We introduce the angle $\displaystyle\vartheta$ such that
$\displaystyle p_{x}=p\cos\vartheta$ and $\displaystyle p_{y}=p\sin\vartheta$,
with $\displaystyle p$ the norm of $\displaystyle P$. The optimal control is
given by $\displaystyle\frac{\partial H_{P}}{\partial\varphi}$ and fulfills
$\displaystyle\tan\varphi=-\frac{p_{x}(t)}{p_{y}(t)}$. We arrive at
$\displaystyle\tan\varphi=-\tan(\frac{\pi}{2}-\vartheta-\omega t)$ i.e.

$\displaystyle\varphi(t)=\vartheta+\omega t-\frac{\pi}{2}$

The optimal dynamics for the state of the system are then given by

$\displaystyle\dot{Z}=i\omega Z-e^{i\vartheta}e^{i\omega t}$

with $\displaystyle Z=x+iy$. The solution can be integrated exactly

$\displaystyle Z(t)=e^{i\omega t}[Z(0)-e^{i\vartheta}t]$

We stress that the optimal trajectories depend only on one free parameter
$\displaystyle\vartheta$ which is fixed by the initial adjoint state. In the
case $\displaystyle Z(0)=0$ and $\displaystyle Z_{f}=1$ under study, we deduce
that

$\displaystyle e^{-i\omega t_{f}}=-t_{f}e^{i\vartheta}$

which leads to the minimum time $\displaystyle t_{f}^{(c)}=1$, and
$\displaystyle\vartheta=\pi-\omega t_{f}$.

We consider now a piecewise constant control wit $\displaystyle t_{f}=NT$. The
maximization condition is given by

$\displaystyle\tan\varphi_{k}=-\frac{\int_{kT}^{(k+1)T}p_{x}}{\int_{kT}^{(k+1)T}p_{y}}$

which leads after integration to

$\displaystyle\varphi_{k}=\vartheta+(k+\frac{1}{2})\omega T-\frac{\pi}{2}$

where $\displaystyle\vartheta$ has the same definition as in the continuous
limit. This solution corresponds to the average of the continuous function
$\displaystyle\varphi(t)$ in the interval $\displaystyle[kT,(k+1)T]$. The
differential equation becomes in this case

$\displaystyle\dot{Z}=i\omega Z-e^{i\vartheta}e^{i(k+\frac{1}{2})\omega T}$

We have in the interval $\displaystyle[kT,(k+1)T]$

$\displaystyle Z_{k+1}=e^{i\omega T}Z_{k}-\frac{2}{\omega}\sin(\frac{\omega
T}{2})e^{i(k+1)\omega T}e^{i\vartheta}$

with $\displaystyle Z_{k}=Z(kT)$. For the initial state $\displaystyle Z_{0}$,
the final state is given by

$\displaystyle Z_{N}=e^{iN\omega T}Z_{0}-\frac{2N}{\omega}\sin(\frac{\omega
T}{2})e^{i\vartheta}e^{iN\omega T}$

For the control problem, we deduce that

$\displaystyle e^{-iN\omega T}=-\frac{2N}{\omega}\sin(\frac{\omega
T}{2})e^{i\vartheta}$

The parameter $\displaystyle T$ fulfills the condition

$\displaystyle T=\frac{2}{\omega}\arcsin(\frac{\omega}{2N})$

The total time is thus

$\displaystyle t_{f}=NT=\frac{2N}{\omega}\arcsin(\frac{\omega}{2N})$

The evolution of this time with respect to the minimum time in the continuous
limit ($\displaystyle t_{f}^{(c)}=1$) when $\displaystyle N\gg 1$ is then

$\displaystyle\frac{t_{f}}{t_{f}^{(c)}}-1\simeq\frac{\omega^{2}}{24N^{2}}$

We show here that the convergence in this linear approximation is polynomial
and not exponential as in the quantum case. These analytical results are
illustrated in Fig. 6. We conjecture that this speed of convergence is the
same in any linear control system. The exponential convergence observed in
quantum systems seems due to their bilinear structure, and therefore has an
intrinsic link with the quantum nature of their dynamics.

Figure 6: Evolution of the minimum time $\displaystyle t_{f}$ as a function of
the number of steps $\displaystyle N$ (Panel (a)). The solid line is just to
guide the eye. The panel (b) displays the optimal controls in the continuous
limit (red or gray curve) and in the discrete case (black curve) for
$\displaystyle N=4$. The parameter $\displaystyle\omega$ is set to 0.5 and the
minimum time is equal to 1.0008 for $\displaystyle N=4$. Quantities plotted
are dimensionless

##  6 Gradient-based optimal control algorithm

We describe in this section a formulation of gradient-based algorithm for
piecewise constant pulses. We show how this algorithm is linked to the
discrete version of the PMP.

We consider a general quantum system of state $\displaystyle|\psi(t)\rangle$
solution of the Schrödinger equation

$\displaystyle
i|\dot{\psi}\rangle=[\cos(\varphi)H_{x}+\sin(\varphi)H_{y}]|\psi\rangle$

with $\displaystyle|\psi(0)\rangle=|\psi_{0}\rangle$ the initial state, the
control being $\displaystyle\varphi(t)\in\mathbb{R}$. The optimal control
problem is defined through the figure of merit $\displaystyle J$ to maximize

$\displaystyle J=|\langle\psi_{f}|\psi(t_{f})\rangle|^{2}$

where $\displaystyle|\psi_{f}\rangle$ is the target state and the control time
$\displaystyle t_{f}$ is fixed. The control parameter $\displaystyle\varphi$
is a piecewise constant function of values $\displaystyle(\varphi_{k})$,
$\displaystyle k\in\\{0,1,\cdots,N-1\\}$, with a sampling period
$\displaystyle T$. The final state of the dynamics can be expressed as

$\displaystyle|\psi(t_{f})\rangle=|\psi_{N}\rangle=U_{N-1}U_{N-2}\cdots
U_{0}|\psi_{0}\rangle$

where $\displaystyle
U_{k}=\exp[-i(\cos(\varphi_{k})H_{x}+\sin(\varphi_{k})H_{y})T]$. We denote by
$\displaystyle|\psi_{k}\rangle=U_{k-1}U_{k-2}\cdots U_{0}|\psi_{0}\rangle$ and
$\displaystyle|\chi_{k}\rangle=U_{k}^{\dagger}U_{k+1}^{\dagger}\cdots
U_{N-1}^{\dagger}|\psi_{f}\rangle$, the state and the adjoint state at step
$\displaystyle k$. The figure of merit $\displaystyle J$ is then given by
$\displaystyle J=|\langle\chi_{k+1}|U_{k}|\psi_{k}\rangle|^{2}$. The next step
consists in computing the gradient of $\displaystyle J$ with respect to the
control $\displaystyle\varphi$, which can be found from the derivative of
$\displaystyle U_{k}$ with respect to $\displaystyle\varphi_{k}$. For that
purpose, we use the Wilcox formula for the derivative of the exponential of a
matrix [27]

$\displaystyle\frac{\partial
e^{tA}}{\partial\theta}=e^{tA}\int_{0}^{t}e^{-t^{\prime}A}\frac{\partial
A}{\partial\theta}e^{t^{\prime}A}dt^{\prime}$

In our case, we get:

$\displaystyle\displaystyle\frac{\partial}{\partial\varphi_{k}}\exp[-i\tau(\cos(\varphi_{k})H_{x}+\sin(\varphi_{k})H_{y})]=$
$\displaystyle\displaystyle
U_{k}(\tau)\int_{0}^{\tau}U_{k}^{\dagger}(t)(-i)[-\sin(\varphi_{k})H_{x}+\cos(\varphi_{k})H_{y}]U_{k}(t)dt$

which leads to

$\displaystyle\displaystyle\frac{\partial J}{\partial\varphi_{k}}=$
$\displaystyle\displaystyle-2\sin(\varphi_{k})\Im[\langle\psi_{N}|\psi_{f}\rangle\langle\chi_{k}|\int_{0}^{\tau}U_{k}^{\dagger}(t)H_{x}U_{k}(t)dt|\psi_{k}\rangle]$
$\displaystyle\displaystyle+2\cos(\varphi_{k})\Im[\langle\psi_{N}|\psi_{f}\rangle\langle\chi_{k}|\int_{0}^{\tau}U_{k}^{\dagger}(t)H_{y}U_{k}(t)dt|\psi_{k}\rangle]$

and finally to

$\displaystyle\displaystyle\frac{\partial
J}{\partial\varphi_{k}}=2\int_{0}^{\tau}\Im[\langle\psi_{N}|\psi_{f}\rangle\langle\chi_{k}(t)[-\sin(\varphi_{k})H_{x}$
$\displaystyle\displaystyle+\cos(\varphi_{k})H_{y}]|\psi_{k}(t)\rangle]dt$

or in a compact form:

$\displaystyle\frac{\partial
J}{\partial\varphi_{k}}=-\sin(\varphi_{k})I_{x}^{(k)}+\cos(\varphi_{k})I_{y}^{(k)}$

with

$\begin{split}&I_{x}^{(k)}=2\int_{0}^{\tau}\Im[\langle\psi_{N}|\psi_{f}\rangle\langle\chi_{k}(t)|H_{x}|\psi_{k}(t)\rangle]dt\\\
&I_{y}^{(k)}=2\int_{0}^{\tau}\Im[\langle\psi_{N}|\psi_{f}\rangle\langle\chi_{k}(t)|H_{y}|\psi_{k}(t)\rangle]dt\end{split}$

The maximization condition is achieved when all the gradients are zero. We
recover here a formula for the optimal control very similar as the one given
by the discrete version of the PMP

$\tan(\varphi_{k})=\frac{I_{y}^{(k)}}{I_{x}^{(k)}}$ (5)

This derivation gives a qualitative justification of the extended PMP
formulation for piecewise constant pulses. Note that this argument can be
applied to other expressions of the Hamiltonian of the system. Equation (5)
also allows to define a new version of GRAPE which has the key advantage of
giving an exact value of the gradient [7]. We now study the control problem
described in Sec. 3 with
$\displaystyle|\psi(0)\rangle=\frac{1}{\sqrt{2}}[|1\rangle+|2\rangle]$ and
$\displaystyle|\psi_{f}\rangle=\frac{1}{\sqrt{2}}[|1\rangle+i|2\rangle]$ where
$\displaystyle\\{|1\rangle,|2\rangle\\}$ is the canonical basis of the Hilbert
space. We compare in Fig. 7 this formulation to the standard one based on a
split-operator method [9] and to the recent auxiliary matrix approach [28]. We
point out that this latter method gives also an exact gradient but at the cost
of increasing by a factor 2 the size of the system, while the first only
approximates the gradient of $\displaystyle J$. In order to get a fair
comparison between the algorithms, the same parameters in the optimization
procedure have been used. We have checked that the qualitative conclusions do
not depend on a specific choice of the numerical parameters. We choose to plot
the fidelity $\displaystyle
d(t_{f})=1-|\langle\psi_{f}|\psi(t_{f})\rangle|^{2}$ as a function of the wall
time. Similar results have been achieved for the CPU time. Note that the
algorithms based on the PMP and on the auxiliary matrix method give exactly
the same performance when the figure of merit is plotted in terms of the
number of iterations. This observation was expected since the gradient is
computed without any approximation in the two cases. The slight difference
observed with the wall time is due to the computation time for calculating the
gradient. Here the PMP version gives better results because the integrals
involved in Eq. (5) can be determined exactly. However, the auxiliary matrix
approach seems to be more efficient in terms of computational time if applied
to large quantum systems. This aspect which is beyond the scope of this paper
should be an interesting point to check. A comparison with an automatic
differentiation method [29, 30, 31] has also been performed. The same final
distance to the target state is reached but with longer computation times than
the PMP or the auxilary matrix method. Another key advantage of the PMP
formulation of GRAPE lies in the fact that it is not restricted to bilinear
dynamics, as it is the case, e.g., for the auxiliary matrix approach. The same
algorithm can be used if nonlinearities are added to the model system [32, 33,
34, 35, 36].

Figure 7: Comparison of the different formulations of GRAPE for the control
problem described in Sec. 3. The figure of merit $\displaystyle d$ is defined
as $\displaystyle d(t_{f})=1-|\langle\psi_{f}|\psi(t_{f})\rangle|^{2}$. The
red, black and blue lines depict the evolution of $\displaystyle d$ as a
function of the wall time for respectively the split-operator, the auxiliary
matrix and the discrete PMP approaches. The parameter $\displaystyle N$ is
fixed to 7. Quantities plotted are dimensionless.

Finally, we propose a numerical analysis of the minimum time with the GRAPE
algorithm based on the auxiliary matrix approach. Note that the same results
in terms of precision are obtained with the PMP formulation of GRAPE. For a
fixed control time $\displaystyle t_{f}$, the goal of the optimization
procedure is to minimize the figure of merit $\displaystyle d(t_{f})$. A good
estimation of the optimal solution can be achieved by running the algorithm
with many different initializations. This operation is repeated for a range of
final times $\displaystyle t_{f}$ in order to estimate the minimum time which
is the lower value of $\displaystyle t_{f}$ for which $\displaystyle
d(t_{f})\simeq 0$. In the numerical simulations displayed in Fig. 8, we
observe that the figure of merit $\displaystyle d$ falls very quickly towards
zero around $\displaystyle t_{f}=2.753$, which is the minimum time of the
control process. The distance to the target is of the order of $\displaystyle
10^{-6}$. For the shooting algorithm, we obtain in the same conditions a
minimum time of 2.75292 and a distance of $\displaystyle 10^{-9}$.

Figure 8: Numerical results (crosses) obtained with GRAPE for $\displaystyle
N=3$ and $\displaystyle\delta T=T$ in the case of the process with two
controls. The control times range from 2.72 to 2.77. The vertical solid red
line indicates the value of the minimum time found by this optimization
procedure. Quantities plotted are dimensionless.

##  7 Conclusion

We have applied an extension of the PMP to piecewise constant pulses to the
control in minimum time of two-level quantum systems. We consider as
illustrative examples two benchmark control problems for which the control
procedure is known exactly in the continuous limit. We compare for different
state-to-state transfers the minimum times, the optimal trajectories and the
controls. We show in particular that the minimum time in the discrete case
converges very quickly towards its continuous limit, while the piecewise and
the continuous controls remain quite different. In the examples under study,
we observe that this convergence is exponential. These results can be
considered as good news from an experimental point of view since experimental
sampling limitations only slightly affect the minimum control time and
therefore the overall duration of the operation composed of several control
processes.

Finally, these results pave also the way for other studies using the same
approach, such as the generation of one-qubit gates [50, 18, 49], the control
of open quantum systems [52, 51, 19, 20], the selective [59, 58, 37, 60] and
robust [38, 53, 56, 57, 54, 55] control of two-level quantum systems or the
control of more complex quantum systems [3, 61, 62, 63]. It would be
interesting to verify in these examples that the convergence is always
exponential and does not depend on the control problem under study.

Acknowledgment
This research has been supported by the ANR project “QuCoBEC”
ANR-22-CE47-0008-02. We gratefully acknowledge useful discussions with J.
Billy, D. Guéry-Odelin and B. Peaudecerf, who suggested the idea of this
paper.

##  Appendix A The shooting method in the case of two controls

We describe in this section the solutions of the shooting algorithm in the
case of two controls. We first recall some results about the continuous limit
[7, 64]. In this limit, the optimal equations can be integrated exactly by
introducing spherical coordinates $\displaystyle(r,\vartheta,\varphi)$ such
that $\displaystyle x=r\sin\vartheta\cos\varphi$, $\displaystyle
y=r\sin\vartheta\sin\varphi$ and $\displaystyle z=r\cos\vartheta$, with by
definition $\displaystyle r=1$. Using the generating function [65]
$\displaystyle
F_{2}=p_{x}r\sin\vartheta\cos\varphi+p_{y}r\sin\vartheta\sin\varphi+p_{z}r\cos\vartheta$,
we deduce that the conjugate momenta can be expressed as

$\displaystyle\displaystyle
p_{r}=p_{x}\sin\vartheta\cos\varphi+p_{y}\sin\vartheta\sin\varphi+p_{z}\cos\vartheta$
$\displaystyle\displaystyle
p_{\vartheta}=p_{x}\cos\vartheta\cos\varphi+p_{y}r\cos\vartheta\sin\varphi-
p_{z}r\sin\vartheta$ $\displaystyle\displaystyle
p_{\varphi}=-p_{x}r\sin\vartheta\sin\varphi+p_{y}r\sin\vartheta\cos\varphi$

It is straightforward to show that $\displaystyle p_{x}(0)=p_{r}(0)$,
$\displaystyle p_{y}(0)=p_{\varphi}(0)$ and $\displaystyle
p_{z}(0)=-p_{\theta}(0)$ since $\displaystyle\vartheta(0)=\frac{\pi}{2}$ and
$\displaystyle\varphi(0)=0$. In the spherical coordinates, the true Pontryagin
Hamiltonian (obtained by replacing the controls by their values maximizing
$\displaystyle H_{P}$) can be written as

$\displaystyle
H_{P}=\sqrt{p_{\vartheta}^{2}+p_{\varphi}^{2}/\tan^{2}\vartheta}$

with the constraint $\displaystyle H_{P}=1$ at any time. The Hamiltonian’s
equations then read

$\displaystyle\displaystyle\dot{\vartheta}=p_{\vartheta};~{}\dot{\varphi}=p_{\varphi}/\tan^{2}\vartheta$
$\displaystyle\displaystyle\dot{p}_{\vartheta}=p^{2}_{\varphi}/\tan\vartheta(1+\frac{1}{\tan^{2}\vartheta});~{}\dot{p}_{\varphi}=0$

At time $\displaystyle t=0$, using $\displaystyle\vartheta(0)=\frac{\pi}{2}$
and $\displaystyle H_{P}=1$, we deduce that $\displaystyle
p_{\vartheta}(0)=\pm 1$. The extrema of the trajectory are reached when
$\displaystyle t=\frac{\pi\sqrt{3}}{4}$ and $\displaystyle\dot{\vartheta}=0$.
Using
$\displaystyle\dot{\vartheta}^{2}=1-\frac{p_{\varphi}^{2}}{\tan^{2}\vartheta}$,
we deduce that the corresponding angle $\displaystyle\vartheta_{m}$ satisfies
$\displaystyle\tan\vartheta_{m}=p_{\varphi}$. Starting from the time evolution
of $\displaystyle\vartheta$ given, when $\displaystyle\vartheta(t)$ increases,
by [64]

$\displaystyle\vartheta(t)=\frac{\pi}{2}+\textrm{asin}[\frac{1}{\sqrt{1+p_{\varphi}^{2}}}\sin(\sqrt{1+p_{\varphi}^{2}}t),$

we arrive at $\displaystyle p_{\varphi}=\frac{1}{\sqrt{3}}$. Note that the
value of $\displaystyle p_{r}$ does not play any role in the control process.
Finally, we obtain the possible values of the initial adjoint state
$\displaystyle P(0)=(p_{x}(0),\frac{1}{\sqrt{3}},\pm 1)$ given in the main
text of the paper. We choose to use a sphere to represent the different
initial adjoint states with the following spherical coordinates $\displaystyle
p_{x}=R_{p}\sin\Theta_{p}\cos\Phi_{p}$, $\displaystyle
p_{y}=R_{p}\sin\Theta_{p}\sin\Phi_{p}$ and $\displaystyle
p_{z}=R_{p}\cos\Theta_{p}$. The constant radius is equal to $\displaystyle
R_{p}=\sqrt{\frac{4}{3}+p_{x}^{2}(0)}$ and the angles
$\displaystyle\Theta_{p}$ and $\displaystyle\Phi_{p}$ can be expressed when
$\displaystyle p_{z}(0)=-1$ as

$\displaystyle\displaystyle\Theta_{p}=\textrm{acos}[-1/\sqrt{\frac{4}{3}+p_{x}^{2}(0)}]$
$\displaystyle\displaystyle\Phi_{p}=\textrm{atan}[\frac{1}{\sqrt{3}p_{x}(0)}]$

The corresponding curve is plotted in Fig. 9. The result of the control
process for the discrete case can also be represented on the same sphere. To
this aim, we consider all the possible values of $\displaystyle P(0)$ and we
integrate the trajectory up to the minimum time $\displaystyle t_{f}$. The
figure of merit $\displaystyle d$ is defined as the Euclidian distance from
the final state to the target. A contour plot of the corresponding performance
is displayed in Fig. 9. The optimal adjoint states correspond to
$\displaystyle d=0$. We observe that they also form a curve in the space
$\displaystyle(\Theta_{p},\Phi_{p})$ which is very close to the one of the
continuous case. Note that all these initial adjoint states lead to the same
controls in the discrete version of the PMP.

Figure 9: Figure of merit $\displaystyle d$ for the minimum discrete time
$\displaystyle t_{f}$ ($\displaystyle N=3$ and $\displaystyle\delta T=T$) as a
function of the coordinates $\displaystyle(\Theta_{p},\Phi_{p})$ of the
initial adjoint state (see the text for details). The solid and dashed lines
depict respectively the minimum value of $\displaystyle d$ and the continuous
limit. Quantities plotted are dimensionless.

##  Appendix B Application of the PMP in the case of one control

We show in this section how the PMP in the discrete case can be applied to the
control problem of Sec. 4. We introduce the switching function
$\displaystyle\Phi=p_{z}y-p_{y}z$. We denote by $\displaystyle\Phi_{k}(t)$ the
value of the switching function in the interval $\displaystyle[kT,(k+1)T]$.
The maximization condition can be written as follows:

$\displaystyle(\mathcal{U}-\omega_{k})\frac{1}{T}\int_{kT}^{(k+1)T}\Phi_{k}(t)dt\leq
0,~{}\forall\mathcal{U}\in[-1,1]$

where $\displaystyle\omega_{k}$ is the constant value of
$\displaystyle\omega(t)$ in this interval. We have the following cases, with
$\displaystyle\Gamma_{k}(\omega_{k})=\frac{1}{T}\int_{kT}^{(k+1)T}\Phi_{k}(t)dt$:

  * •

If $\displaystyle\min_{\omega_{k}}\Gamma_{k}>0$ then
$\displaystyle\Gamma_{k}(\omega_{k})>0$ and
$\displaystyle\mathcal{U}-\omega_{k}\leq 0$. We deduce that
$\displaystyle\omega_{k}=+1$.

  * •

If $\displaystyle\max_{\omega_{k}}\Gamma_{k}<0$ then
$\displaystyle\Gamma_{k}(\omega_{k})<0$ and
$\displaystyle\mathcal{U}-\omega_{k}\geq 0$. We deduce that
$\displaystyle\omega_{k}=-1$.

  * •

If $\displaystyle\max_{\omega_{k}}\Gamma_{k}>0$ and
$\displaystyle\min_{\omega_{k}}\Gamma_{k}<0$ then the only solution is
$\displaystyle\Gamma_{k}(\omega_{k})=0$ which leads to a value
$\displaystyle\omega_{k}$ in the interval $\displaystyle[-1,1]$.

Note that it is possible to get solutions distinct from -1 or +1 in the
discrete case while the optimal trajectory is regular in the continuous limit,
i.e. $\displaystyle|\omega(t)|=1$ except in isolated times.

The expression of $\displaystyle\Gamma$ can be derived explicitly as follows.
To simplify the notations, we consider a time interval $\displaystyle[0,T]$
such as the state and the adjoint state are respectively of coordinates
$\displaystyle(x_{0},y_{0},z_{0})$ and $\displaystyle(p_{x0},p_{y0},p_{z0})$
at $\displaystyle t=0$. We also introduce the coordinates of the angular
momentum $\displaystyle L=X\times P$

$\displaystyle\displaystyle L_{x}=yp_{z}-zp_{y}=\Phi$
$\displaystyle\displaystyle L_{y}=zp_{x}-xp_{z}$ $\displaystyle\displaystyle
L_{z}=xp_{y}-yp_{x}$

We have

$\displaystyle\displaystyle x(t)=x_{0}-\frac{\Delta
B}{\Omega_{0}}-\frac{\Delta}{\Omega_{0}}(A\sin(\Omega_{0}t)-B\cos(\Omega_{0}t))$
$\displaystyle\displaystyle y(t)=A\cos(\Omega t)+B\sin(\Omega t)$
$\displaystyle\displaystyle z(t)=z_{0}+\frac{\omega
B}{\Omega_{0}}+\frac{\omega}{\Omega_{0}}(A\sin(\Omega_{0}t)-B\cos(\Omega_{0}t))$

with $\displaystyle A=y_{0}$, $\displaystyle B=\frac{\Delta x_{0}-\omega
z_{0}}{\Omega_{0}}$ and
$\displaystyle\Omega_{0}=\sqrt{\omega^{2}+\Delta^{2}}$. For the adjoint state,
we get similar equations

$\displaystyle\displaystyle p_{x}(t)=p_{x0}-\frac{\Delta
D}{\Omega_{0}}-\frac{\Delta}{\Omega_{0}}(C\sin(\Omega_{0}t)-D\cos(\Omega_{0}t))$
$\displaystyle\displaystyle p_{y}(t)=C\cos(\Omega_{0}t)+D\sin(\Omega_{0}t)$
$\displaystyle\displaystyle p_{z}(t)=p_{z0}+\frac{\omega
D}{\Omega_{0}}+\frac{\omega}{\Omega_{0}}(C\sin(\Omega_{0}t)-D\cos(\Omega_{0}t))$

with $\displaystyle C=p_{y0}$ and $\displaystyle D=\frac{\Delta
p_{x0}-up_{z0}}{\Omega_{0}}$.

The switching function can then be expressed as

$\displaystyle\displaystyle\Phi=\frac{\omega}{\Omega_{0}^{2}}(\Delta
L_{z0}+\omega L_{x0})+\frac{\Delta}{\Omega_{0}^{2}}(\Delta L_{x0}-\omega
L_{z0})\cos(\Omega_{0}t)$
$\displaystyle\displaystyle-\frac{\Delta}{\Omega_{0}}L_{y0}\sin(\Omega_{0}t)$

The corresponding function
$\displaystyle\Gamma=\frac{1}{T}\int_{0}^{T}\Phi(t)dt$ is given by

$\displaystyle\displaystyle\Gamma=\frac{\omega}{\Omega_{0}^{2}}(\Delta
L_{z0}+\omega L_{x0})+\frac{\Delta}{\Omega_{0}^{2}}(\Delta L_{x0}-\omega
L_{z0})\textrm{sinc}(\Omega_{0}T)$
$\displaystyle\displaystyle+\frac{\Delta}{\Omega_{0}}L_{y0}\frac{\cos(\Omega_{0}T)-1}{\Omega_{0}T}$

As expected, we obtain in the limit $\displaystyle T\to 0$ that
$\displaystyle\Gamma\simeq\Phi$. In the general case, $\displaystyle\Gamma$ is
a function of $\displaystyle\omega$ in the interval $\displaystyle[-1,1]$. At
each time step, we compute numerically the maximum, the minimum and the zeros
of $\displaystyle\Gamma$ to find the corresponding value of
$\displaystyle\omega$.

## References

  * [1] S. J. Glaser, U. Boscain, T. Calarco, C. Koch, W. Kockenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny and F. Wilhelm, Eur. Phys. J. D 69, 279 (2015)
  * [2] C. Brif, R. Chakrabarti and H. Rabitz, New J. Phys. 12, 075008 (2010)
  * [3] C.P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbruggen, D. Sugny and F. K. Wilhelm, EPJ Quantum Technol. 9, 19 (2022)
  * [4] N. V. Vitanov, T. Halfmann, B. W. Shore, and K. Bergmann, Annual review of physical chemistry 52, 763 (2001).
  * [5] D. Stefanatos and E. Paspalakis, Eur. Phys. Lett. 132, 60001 (2020)
  * [6] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019)
  * [7] U. Boscain, M. Sigalotti and D. Sugny, PRX Quantum 2, 030203 (2021)
  * [8] L. S. Pontryagin et al., _The Mathematical Theory of Optimal Processes_ (John Wiley and Sons, New York, London, 1962)
  * [9] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen and S. J. Glaser, J. Magn. Reson. 172, 296 (2005)
  * [10] D. M. Reich, M. Ndong and C. P. Koch, J. Chem. Phys. 136, 104103 (2012)
  * [11] J. Werschnik and E. K. U. Gross, J. Phys. B 40, R175 (2007)
  * [12] P. Doria, T. Calarco, and S. Montangero, Phys Rev. Lett. 106, 190501 (2011).
  * [13] D. D’Alessandro, IEEE Trans. Autom. Control 46, 866 (2001)
  * [14] U. Boscain and P. Mason, J. Math. Phys. 47, 062101 (2006)
  * [15] G. C. Hegerfeldt, Phys. Rev. Lett. 111, 260501 (2013)
  * [16] N. Khaneja, R. Brockett and S. J. Glaser, Phys. Rev. A 63, 032308 (2001)
  * [17] N. Khaneja, S. J. Glaser and R. Brockett, Phys. Rev. A 65, 032301 (2002)
  * [18] A. Garon, S. J. Glaser and D. Sugny, Phys. Rev. A 88, 043422 (2013)
  * [19] M. Lapert, E. Assémat, S. J. Glaser and D. Sugny, Phys. Rev. A 88, 033407 (2013)
  * [20] M. Lapert, Y. Zhang, M. Braun, S. J. Glaser and D. Sugny, Phys. Rev. Lett. 104, 083001 (2010)
  * [21] S. Deffner and S. Campbell, J. Phys. A: Math. Theor. 50, 453001 (2017)
  * [22] M. R. Frey, Quantum Inf. Process., 15, 3919 (2016)
  * [23] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Phys. Rev. Lett. 103, 240501 (2009)
  * [24] L. Bourdin and E. Trélat, SIAM J. Control Optim. 59, 489 (2021)
  * [25] L. Bourdin and E. Trélat, Math. Control and Related fields 6, 53 (2016)
  * [26] L. Bourdin and E. Trélat, Automatica 79, 273 (2017)
  * [27] H. Tsai and K. S. Chan, Bernoulli 9, 895 (2003)
  * [28] D. L. Goodwin and I. Kuprov, J. Chem. Phys. 143, 084113 (2015)
  * [29] N. Leung, M. Abdelhafez, J. Koch and D. Schuster, Phys. Rev. A 95, 042318 (2017)
  * [30] M. Abdelhafez, D. I. Schuster, and J. Koch, Phys. Rev. A 99, 052327 (2019)
  * [31] M. H. Goerz, S. C. Carrasco and V. S. Malinovsky, Quantum 6, 871 (2022)
  * [32] M. Lapert, R. Tehini, G. Turinici, and D. Sugny, Phys. Rev. A 78, 023408 (2008)
  * [33] Y. Ohtsuki and K. Nakagami, Phys. Rev. A 77, 033414 (2008)
  * [34] G. Jäger, D. M. Reich, M. H. Goerz, C. P. Koch, and U. Hohenester, Phys. Rev. A 90, 033628 (2014).
  * [35] U. Hohenester, P. K. Rekdal, A. Borzi, and J. Schmiedmayer, Phys. Rev. A 75, 023602 (2007)
  * [36] D. Hocker, J. Yan, and H. Rabitz, Phys. Rev. A 93, 053612 (2016)
  * [37] L. Van Damme, Q. Ansel, S. J. Glaser, and D. Sugny, Phys. Rev. A 98, 043421 (2018)
  * [38] L. Van Damme, Q. Ansel, S. J. Glaser and D. Sugny, Phys. Rev. A 95, 063403 (2017)
  * [39] M. A. Bernstein, K. F. King, and X. J. Zhou, _Handbook of MRI Pulse Sequences_ (Academic, New York, 2004)
  * [40] M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, and O. Morsch, Nat. Phys. 8, 147 (2012)
  * [41] A. Zenesini, H. Lignier, G. Tayebirad, J. Radogostowicz, D. Ciampini, R. Mannella, S. Wimberger, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 103, 090403 (2009)
  * [42] G. Tayebirad, A. Zenesini, D. Ciampini, R. Mannella, O. Morsch, E. Arimondo, N. Lorch, and S. Wimberger, Phys. Rev. A 82, 013633 (2010)
  * [43] L. Garcia, J. M. Bofill, I. Moreira, and G. Albareda, Phys. Rev. Lett. 129, 180402 (2022)
  * [44] M. Poggi, F. C. Lombardo, and D. Wisniacki, Europhys. Lett. 104, 40005 (2013).
  * [45] J.-S. Li, J. Ruths, and S. J. Glaser, Nat. Commun. 8, 446 (2017)
  * [46] J.-S. Li, IEEE Trans. A. C. 56, 345 (2011)
  * [47] V. Martikyan, A. Devra, D. Guéry-Odelin, S. J. Glaser, and D. Sugny, Phys. Rev. A 102, 053104 (2020)
  * [48] K. Kobzar, T. E. Skinner, N. Khaneja, S. J. Glaser and B. Luy, J. Magn. Reson. 194, 58 (2008)
  * [49] A. D. Boozer, Phys. Rev. A 85, 012317 (2012)
  * [50] K. Kozbar, S. Ehni, T. E. Skinner, S. J. Glaser, and B. Luy, J. Magn. Reson. 225, 142 (2012)
  * [51] V. Mukherjee, A. Carlini, A. Mari, T. Caneva, S. Montangero, T. Calarco, R. Fazio, and V. Giovannetti, Phys. Rev. A 88, 062326 (2013)
  * [52] N. Khaneja, B. Luy, and S. J. Glaser, Proc. Natl. Acad. Sci. USA 100, 13162 (2003)
  * [53] D. Daems, A. Ruschhaupt, D. Sugny and S. Guérin, Phys. Rev. Lett. 111, 050404 (2013)
  * [54] J. Zeng and E. Barnes, Phys. Rev. A 98, 012301 (2018)
  * [55] J. Zeng, C. H. Yang, A. S. Dzurak, and E. Barnes, Phys. Rev. A 99, 052321 (2019)
  * [56] A. Ruschhaupt, X. Chen, D. Alonso and J. G. Muga, New J. Phys. 14, 093040 (2012)
  * [57] G. T. Genov, D. Schraft, T. Halfmann, and N. V. Vitanov, Phys. Rev. Lett. 113, 043001 (2014)
  * [58] M. S. Silver, R. I. Joseph, and D. I. Hoult, Phys. Rev. A 31, 2753 (1985)
  * [59] M. S. Silver, R. I. Joseph, C.-N. Chen, V. J. Sank, and D. I. Hoult, Nature (London) 310, 681 (1984)
  * [60] D. Basilewitsch, H. Yuan, and C. P. Koch, Phys. Rev. Research 2, 033396 (2020)
  * [61] N. Dupont, G. Chatelain, L. Gabardos, M. Arnal, J. Billy, B. Peaudecerf, D. Sugny and D. Guéry-Odelin, PRX Quantum 2, 040303 (2021)
  * [62] M. Lapert, G. Ferrini and D. Sugny, Phys. Rev. A 85 023611 (2012)
  * [63] M. Nimbalkar, R. Zeier, J. L. Neves, S. B. Elavarasi, H. Yuan, N. Khaneja, K. Dorai and S. J. Glaser, Phys. Rev. A 85, 012325 (2012)
  * [64] D. Sugny and C. Kontz, Phys. Rev. A 77, 063420 (2008)
  * [65] H. Goldstein, _Classical Mechanics_ (Addison-Wesley, Reading, MA, 1950).