# Optimal Control for Closed and Open System Quantum Optimization

Lorenzo Campos Venuti  Department of Physics & Astronomy, University of
Southern California, Los Angeles, California 90089, USA  Center for Quantum
Information Science & Technology, University of Southern California, Los
Angeles, California 90089, USA  Domenico D’Alessandro  Department of
Mathematics, Iowa State University, Ames, Iowa 50014, USA  Daniel A. Lidar
Department of Physics & Astronomy, University of Southern California, Los
Angeles, California 90089, USA  Center for Quantum Information Science &
Technology, University of Southern California, Los Angeles, California 90089,
USA  Department of Electrical & Computer Engineering, University of Southern
California, Los Angeles, CA 90089, USA  Department of Chemistry, University of
Southern California, Los Angeles, CA 90089, USA

(August 27, 2024)

###### Abstract

We provide a rigorous analysis of the quantum optimal control problem in the
setting of a linear combination $s(t)B+(1-s(t))C$ of two noncommuting
Hamiltonians $B$ and $C$. This includes both quantum annealing (QA) and the
quantum approximate optimization algorithm (QAOA). The target is to minimize
the energy of the final “problem” Hamiltonian $C$, for a time-dependent and
bounded control schedule $s(t)\in[0,1]$ and $t\in\mathcal{I}:=[0,t_{f}]$. It
was recently shown, in a purely closed system setting, that the optimal
solution to this problem is a “bang-anneal-bang” schedule, with the bangs
characterized by $s(t)=0$ and $s(t)=1$ in finite subintervals of
$\mathcal{I}$, in particular $s(0)=0$ and $s(t_{f})=1$, in contrast to the
standard prescription $s(0)=1$ and $s(t_{f})=0$ of quantum annealing.

Here we extend this result to the open system setting, where the system is
described by a density matrix rather than a pure state. This is the natural
setting for experimental realizations of QA and QAOA. For finite-dimensional
environments and without any approximations we identify sufficient conditions
ensuring that either the bang-anneal, anneal-bang, or bang-anneal-bang
schedules are optimal, and recover the optimality of $s(0)=0$ and
$s(t_{f})=1$. However, for infinite-dimensional environments and a system
described by an adiabatic Redfield master equation we do not recover the bang-
type optimal solution. In fact we can only identify conditions under which
$s(t_{f})=1$, and even this result is not recovered in the fully Markovian
limit.

The analysis, which we carry out entirely within the geometric framework of
Pontryagin Maximum Principle, simplifies using the density matrix formulation
compared to the state vector formulation. This analysis reveals that the bang-
anneal-bang optimality result requires the assumption that the optimal control
schedule is such that it is possible to lower the cost by increasing the total
time $t_{f}$. A necessary condition for this is that $t_{f}$ is smaller than
the critical time $t_{c}$ needed to reach the ground state of $C$ exactly. In
previous work this condition was believed to also be sufficient, but we give a
counterexample.

We derive a “switching equation” which describes the behavior of the optimal
schedule switches between the two types of bangs and the anneals, and use it
to identify the general features of the optimal control protocols. As an
illustration of the theory, we analyze the simple example of a single
spin-$1/2$, and prove that in this case the optimal solution in the closed
system setting is the bang-bang schedule, switching midway from $s\equiv 0$ to
$s\equiv 1$.

##  I Introduction

There is a great deal of interest in optimization algorithms that can be run
on today’s noisy intermediate scale quantum (NISQ) information processors
Preskill (2018). Two prime examples are quantum annealing (QA) Kadowaki and
Nishimori (1998) and the quantum approximate optimization algorithm (QAOA)
Farhi  _et al._ (2014). Both algorithms switch between two non-commuting
Hamiltonians: a “driver” (or “mixer”) $B$ and a “target” (or “problem”) $C$.
The latter encodes the solution to the optimization problem as its ground
state. The two algorithms can be viewed as complementary: QA switches
continuously while QAOA switches discretely; hence they are particularly well
suited for analog and gate-model devices, respectively. In addition, both
algorithms are related to the quantum adiabatic algorithm Farhi  _et al._
(2000), which is guaranteed by the adiabatic theorem Kato (1950) to converge
to the optimal solution in the limit of arbitrarily long evolution times
Jansen  _et al._ (2007); Lidar  _et al._ (2009); Ge  _et al._ (2016). QA
relaxes the strict adiabaticity condition while retaining continuity Kadowaki
and Nishimori (1998), and the adiabatic algorithm becomes an instance of QAOA
when the continuous evolution is “Trotterized” (replaced by pulsed segments)
(Farhi  _et al._ , 2014, Sec. VI). There have been numerous studies of the two
algorithms, including some that have compared them, with mixed results Bapat
and Jordan (2019); Zhou  _et al._ (2020); Streif and Leib (2019); Pagano  _et
al._ (2020).

In essence, the question of which algorithm performs best – QA or QAOA – boils
down to an optimization of the switching schedule. Various results have
already been established within the framework of the adiabatic algorithm, QA,
or QAOA. For example, it is well known that the adiabatic algorithm can
benefit from schedule optimization, even to an extent that can affect whether
it provides a quantum speedup or not, as in the case of the Grover search
problem Roland and Cerf (2002); Rezakhani  _et al._ (2010). It has also been
established that a variational approach can optimize the adiabatic switching
schedule Rezakhani  _et al._ (2009). Likewise, optimality results are known
for QA Morita and Nishimori (2008); Galindo and Kreinovich (2020) and QAOA
Szegedy (2019). A natural question is whether one can jointly treat QA and
QAOA under a single schedule optimization framework. The first such attempt
was made by Yang et al. Yang  _et al._ (2017) using the framework of the
Pontryagin Maximum Principle (PMP) of optimal control Kipka and Ledyaev (2015)
(see also Refs. Lin  _et al._ (2019); Mbeng  _et al._ (2019a, b)), whose
conclusions favoring a strict QAOA-type schedule were later shown to be overly
restrictive by Brady et al. Brady  _et al._ (2021), who showed that in general
a hybrid discrete-continuous schedule is optimal.

The results of Brady et al. were obtained in a closed system setting of purely
unitary dynamics. Here, we generalize the theory to the open system setting,
and obtain their closed system results as a special case. We proceed to first
provide the general background for the problem, after which we outline the
structure of the rest of the paper.

##  II Background

The closed system setting involves a system evolving unitarily in a
$d$-dimensional Hilbert space $\mathcal{H}$ subject to the Schrödinger
equation:

$\frac{d}{dt}|\psi(t)\rangle=-iH(t)|\psi(t)\rangle,\qquad|\psi(0)\rangle=|\psi_{0}\rangle\
.$ (1)

The protoypical quantum annealing problem concerns finding the optimal
schedule $s(t)\in[0,1]$ for the time-dependent Hamiltonian given by111In the
control literature the notation $u$ or $u(t)$ is used to denote the control
function, rather than $s$ or $s(t)$. Here we choose to use the notation that
is more familiar in the quantum computing community.

$\displaystyle H(t)$ $\displaystyle=s(t)B+(1-s(t))C\ ,\quad t\in\mathcal{I}$
(2a) $\displaystyle=C+s(t)(B-C)\ .$ (2b)

The _control interval_ is $\mathcal{I}=[0,t_{f}]$. Often the Hermitian
operator $C$ is an Ising-type Hamiltonian of the form
$\sum_{i=1}^{n}h_{i}\sigma^{z}_{i}+\sum_{i<j}^{n}J_{ij}\sigma^{z}_{i}\sigma^{z}_{j}$
(where $h_{i}$ and $J_{ij}$ are local longitudinal fields and couplings,
respectively, and $\sigma^{z}_{i}$ is the Pauli matrix acting on the $i$’th
qubit), and the Hermitian operator $B$ is a transverse field of the form
$\sum_{i=1}^{n}\sigma^{x}_{i}$ Kadowaki and Nishimori (1998). For our purposes
it only matters that $[B,C]\neq 0$.

The initial state $|\psi_{0}\rangle$ is assumed to be the ground state of $B$,
and in both QA and QAOA the target state is the ground state of $C$. A
relaxation of this, which we consider as the objective in the present work, is
to minimize the expectation value of $C$ at a given final time $t_{f}$, i.e.,

$J:=\langle\psi(t_{f})|C|\psi(t_{f})\rangle\ .$ (3)

Minimizing $J$ is equivalent to minimizing the energy of the $C$ Hamiltonian,
and if the global minimum is found then this corresponds to finding the ground
state of $C$ (i.e., solving the optimization problem defined by
$\\{h_{i},J_{ij}\\}$ when $C$ is in Ising form).

It is known in quantum control theory (see, e.g., Ref. D’Alessandro (2007))
that if $\hat{\mathfrak{L}}$ is the Lie algebra generated by $B$ and $C$ and
$e^{\hat{\mathfrak{L}}}$ the corresponding Lie group, assumed to be compact,
the set of states reachable from $|\psi_{0}\rangle$ with free final time
$t_{f}$ is

$\mathfrak{R}:=\\{X|\psi_{0}\rangle\,|\,X\in e^{\hat{\mathfrak{L}}}\\},$ (4)

so that the absolute minimum of the cost $J$ is

$J_{\min}=\min_{|\psi\rangle\in\mathfrak{R}}\langle\psi|C|\psi\rangle\ .$ (5)

If the dynamical Lie algebra $\hat{\mathfrak{L}}$ is the whole $su(d)$ then
any state in the Hilbert space can be reached (starting from any other state),
in particular the ground state of $C$, in which case the system is said to be
_controllable_ Jurdjevic and Sussmann (1972). However, requiring full
controllability may be overly restrictive, as we only need to reach a
particular state. The following is a simple generalization that provides a
sufficient condition for reaching the ground state.

######  Proposition 1.

Suppose $[B,P_{0}]=[C,P_{0}]=0$ where $P_{0}$ is an orthogonal projector. The
Hilbert space decomposes according to the block structure
$\mathcal{H}=\mathrm{Ran}(P_{0})\oplus\mathrm{Ker}(P_{0})$. The Lie algebra
generated by $B$ and $C$, $\hat{\mathfrak{L}}$, must have the same block
structure. Suppose that, according to this structure,
$\hat{\mathfrak{L}}=su(d_{0})\oplus\hat{\mathfrak{L}}_{1}$ with unspecified
$\hat{\mathfrak{L}}_{1}$ and $d_{0}=\dim(\mathrm{Ran}(P_{0}))$; then, if the
initial state belongs to $\mathrm{Ran}(P_{0})$, any state in
$\mathrm{Ran}(P_{0})$ can be reached (in finite time).

The proof is self-evident, since the full controllability result Jurdjevic and
Sussmann (1972) is now applicable in $\mathrm{Ran}(P_{0})$. This
generalization can be applied, for example, in case both $B$ and $C$ commute
with a third operator, say $M$, and one knows to which sector of $M$ the
ground state of $C$ belongs; see Appendix A for an example. In any case, it is
clear that something must be assumed in order to guarantee the reachability of
the ground state of $C$. Clearly, a necessary condition is $[B,C]\neq 0$, but
even when $[B,C]\neq 0$ it is easy to come up with examples where the ground
state of $C$ cannot be reached; see Appendix A. In the following we will
tacitly assume that conditions are such that the ground state of $C$ can be
reached.

Brady et al. Brady  _et al._ (2021) used optimal control methods to prove that
for the cost as defined as in Eq. (3), the optimal schedule is one where at
the beginning and end of the control interval $s\equiv 0$ and $s\equiv 1$,
respectively.222We use the notation $f\equiv c$ to mean that the function f is
“identically” equal to $c$ in some interval $\Theta$, i.e., $f\equiv c$ is
equivalent to $f(t)=c$ for $t\in\Theta$. From a quantum annealing perspective
this might appear as a counterintuitive result, since it means that rather
than the usual “forward” formulation of quantum annealing Albash and Lidar
(2018); Hauke  _et al._ (2020), where one interpolates smoothly from $H(0)=B$
to $H(t_{f})=C$, the optimal protocol starts from the system being in the
ground state of $B$ but the initial Hamiltonian is $C$, and the final
Hamiltonian is not $C$ but rather $B$. The result, however, can be understood
by noting that in the adiabatic approach, one interpolates so slowly from
$H(0)=B$ to $H(t_{f})=C$ that the system always remains in the ground state.
Instead, in optimal control, we optimize over the set of possible states
obtained by applying either $C$ or $B$ to the initial state, in a continuous
fashion. In this sense, applying $B$ at the beginning is a waste of time as it
does not change the initial state. Applying $C$ at the end, when the system is
supposed to be close to the ground state of $C$, is similarly wasteful. This
relaxation of the approach of strict adiabaticity is in line with other
alternatives, such as shortcuts to adiabaticity del Campo (2013); Takahashi
(2017) and diabatic quantum annealing Crosson and Lidar (2021).

More precisely, Ref. Brady  _et al._ (2021) showed, provided that a certain
condition holds (see below), that the optimal control function starts (ends)
with $s\equiv 0$ ($s\equiv 1$) in an interval of positive measure after $t=0$
(before $t_{f}$). Elsewhere the optimal control $s(t)$ is “singular”, except
for possible interruptions by a sequence of “bang” controls, where $s\equiv 0$
or $s\equiv 1$. In control theory a “singular” interval or arc, is an interval
of time where the PMP control Hamiltonian in Eq. (10) below does not depend on
the control $s$. The remaining “nonsingular” arcs give rise to the “bang”
controls. In the numerical simulations of Ref. Brady  _et al._ (2021), the
control appeared to be continuous (even smooth) on such singular arcs. Hence
the term “anneal” was used in lieu of “singular”, with the intention of
stressing the continuous (or possibly even smooth) nature of the control on
the singular arcs. They suggestively called the resulting optimal control a
“bang-anneal-bang” protocol. At present, a rigorous proof that the control
function is continuous (let alone smooth) on singular arcs is lacking, and
there is some risk of confusion in interpreting the singular arcs as always
being continuous, or even differentiable as is typically assumed in QA and
adiabatic quantum computing Jansen  _et al._ (2007); Albash and Lidar (2018).
Nonetheless, while keeping these caveats in mind, we shall adopt the same
(numerically supported) terminology as Ref. Brady  _et al._ (2021), and use
“continuous (or anneal) $=$ singular” as well as “bang $=$ nonsingular”
interchangeably.

Here, we consider the _open system_ version of the same optimal control
problem. We reformulate the problem in terms of the density matrix $\rho$,
whose dynamics is described by the following, rather general master
equation:333The form we have assumed is called a time-convolutionless master
equation. The most general master equation is in Nakajima-Zwanzig form and
includes a memory kernel superoperator $\mathcal{K}(t,t^{\prime})$ acting
jointly on the system and the environment $E$, such that (for a factorized
initial condition)
$\dot{\rho}=\operatorname{Tr}_{E}\int_{0}^{t}\mathcal{K}(t,t^{\prime})\rho(t^{\prime})\otimes\rho_{E}dt^{\prime}$,
with $\rho_{E}$ a fixed environment state and $\operatorname{Tr}_{E}$ denoting
the partial trace over the environment Breuer and Petruccione (2002).

$\dot{\rho}=\mathcal{L}\rho,\qquad\rho(0)=\rho_{0}\ ,$ (6)

where the Liouvillian $\mathcal{L}$ depends linearly on the control $s$ (and
the controlled operators $B,C$). Note that the Liouvillian is not explicitly
time-dependent (i.e. $\partial_{t}\mathcal{L}=0,\ \forall t\in\mathcal{I}$)
and depends on time only through the control schedule $s(t)$. This is an
important requirement that will play a crucial role in our ability to apply
the Pontryagin principle in the form we need, as we discuss in more detail
below. Furthermore, to be physically meaningful, $\mathcal{L}$ must preserve
hermiticity, i.e., $[\mathcal{L}(X)]^{\dagger}=\mathcal{L}(X^{\dagger}),\
\forall X$. Instead of Eq. (3), the cost $J$ takes the form

$J:=\operatorname{Tr}\left[C\rho(t_{f})\right]=\langle C,\rho(t_{f})\rangle\
,$ (7)

where we used the Hilbert-Schmidt scalar product $\langle
X,Y\rangle:=\operatorname{Tr}(X^{\dagger}Y)$ for operators $X,Y$. We shall see
that a description and treatment of the optimal control problem in the setting
of the density matrix is not only more general but also more elegant since the
cost $J$ is linear in the state rather than quadratic, as in Eq. (3).
Moreover, we obtain the closed system result as a special case. Unlike Ref.
Brady  _et al._ (2021), which used a mixture of the PMP and a variational
(Lagrange multiplier type of) argument, we use only the PMP, which
significantly simplifies the proof.

The rest of this paper is organized as follows. In Sec. III, we apply general
results from optimal control theory and the necessary conditions of the PMP to
the problem of minimizing $J$ [Eq. (7)] for a given final time $t_{f}$ and the
general dynamical system of the form of Eq. (6). In Sec. IV we specialize to
the case of closed systems, which are described by the von Neumann equation.
In particular, we confirm but also sharpen the results of Ref. Brady  _et al._
(2021). We also analyze in depth the optimal control problem of a single
spin-$1/2$, and prove that the optimal schedule is of the bang-bang type. In
Sec. V we consider the case of open systems. This includes both the most
general case of a reduced description of quantum system obtained by tracing
out the environment it is coupled to, and the case where the open quantum
system is described by adiabatic master equations, both non-Markovian and
Markovian. In Sec. VI we derive a “switching equation,” which allows us to
provide a general characterization of the switches between non-singular and
singular arcs, and derive conditions for the presence or absence of singular
arcs. We also give a heuristic derivation of the shortening of the length of
the arcs between two switches with increasing system size. We conclude in Sec.
VII. In a series of appendices we provide additional background on optimal
control theory and technical details and proofs of various results from the
main text.

##  III Statement of the Pontryagin Maximum Principle

We state the PMP as it applies to our problem of interest (see Appendices B
and C):

######  Theorem 1.

Assume that $\rho^{*}$ and $s^{*}$ are, respectively, an optimal state and
control pair for the problem defined by Eqs. (6) and (7) for a fixed final
time $t_{f}$.444We also use an asterisk to denote complex conjugation; the
meaning will always be clear by context. Then there exists a nonzero $n\times
n$ Hermitian time-dependent matrix $p=p(t)$ called the co-state that
satisfies555$\mathcal{L}^{{\dagger}}$ in Eq. (8) indicates the Hilbert-Schmidt
adjoint of $\mathcal{L}$, defined via
$\langle\mathcal{L}^{\dagger}(A),B\rangle:=\langle
A,\mathcal{L}(B)\rangle=\operatorname{Tr}[A^{\dagger}\mathcal{L}(B)],\ \forall
A,B$; see Appendix C.

$\dot{p}=-\mathcal{L}^{{\dagger}}p\ ,$ (8)

with the final condition

$p(t_{f})=-C\ .$ (9)

Furthermore, define the _PMP control Hamiltonian_ function

$\mathbb{H}\left(p,\rho,s\right):=\langle p,\mathcal{L}\rho\rangle\ .$ (10)

We then have the _maximum principle_ :

$\mathbb{H}\left(p(t),\rho^{*}(t),s^{*}(t)\right)=\max_{v\in[0,1]}\mathbb{H}\left(p(t),\rho^{*}(t),v\right)\
,$ (11)

and there exists a real constant $\lambda$ such that

$\mathbb{H}\left(p(t),\rho^{*}(t),s^{*}(t)\right)\equiv\lambda\ .$ (12)

A few remarks are in order.

  * •

The PMP control Hamiltonian _function_ (10) is, of course, different from the
Hamiltonian _operator_ (2) generating the dynamics.

  * •

Since $p$ and $\rho$ are Hermitian and $\mathcal{L}$ is Hermiticity-preserving
[$\left[\mathcal{L}\left(X\right)\right]^{\dagger}=\mathcal{L}\left(X^{\dagger}\right)$
$\forall X$], “expectation values” of the form $\langle
p,\mathcal{L}\rho\rangle$ are real, and hence so is the PMP control
Hamiltonian (10).

  * •

The condition $\partial_{t}\mathcal{L}=0$ $\forall t\in\mathcal{I}$ must be
satisfied and is implicit in Eq. (6). In other words, $\mathcal{L}$ may not
depend explicitly on time. Without this condition Eq. (12) does not hold with
a constant $\lambda$.

  * •

It is worth highlighting that at the final time the co-state becomes the
(negative of the) target Hamiltonian [Eq. (9)], a fact we use repeatedly in
our applications of the Theorem 1 below.

  * •

As discussed after Theorem 8 in Appendix B, if the optimal trajectory is such
that the constraint of the final time $t_{f}$ is “active”, i.e., a small
perturbation $t_{f}+\delta$ allows us to decrease the cost, then $\lambda>0$
in Eq. (12). This is an important point that will be further discussed in the
next section.

  * •

Given that the PMP is formulated in terms of real-valued quantities in the
optimal control literature (see Appendix B), one must first transform the
relevant equations into real-valued ones. This can easily be done since the
space of $n\times n$ Hermitian matrices is isomorphic to the space of $n^{2}$
real variables via coordinatization. We discuss this in Appendix C.

  * •

Since $p(t)$ and $\rho(t)$ are solution of differential equations, they are
continuous function of time. This implies that expressions of the form
$\langle p,\mathcal{K}\rho\rangle$ with the superoperator $\mathcal{K}$
independent of time (both explicitly and implicitly), are continuous functions
of $t$, a fact which we repeatedly and implicitly use below.

##  IV The closed system case

We first consider the closed system case. Let us define the superoperator

$\mathcal{K}_{X}:=-i\left[X,\bullet\right]\ .$ (13)

Note that $\mathcal{K}_{X}$ is linear with respect to $X$. For Hermitian $X$,
$\mathcal{K}_{X}$ is anti-Hermitian (see Appendix D):

$\mathcal{K}_{X}^{{\dagger}}=-\mathcal{K}_{X^{\dagger}}=-\mathcal{K}_{X}\ .$
(14)

The von Neumann equation corresponding to Eq. (1) is

$\dot{\rho}=\mathcal{K}_{C}\rho+s(t)\left(\mathcal{K}_{B}\rho-\mathcal{K}_{C}\rho\right),\quad\rho(0)=\rho_{0}\
,$ (15)

where henceforth we denote the initial and final conditions of operators $X$
by $X(0):=X_{0}$ and $X(t_{f}):=X_{f}$, respectively. I.e., one has Eq. (6)
with

$\mathcal{L}=\mathcal{K}_{C}+s(t)\mathcal{K}_{B-C}\ .$ (16)

Since in this case $\mathcal{L}^{{\dagger}}=-\mathcal{L}$, Eq. (8) tells us
that the co-state matrix $p$ satisfies the same equation as $\rho$:

$\dot{p}=\mathcal{K}_{C}p+s(t)\mathcal{K}_{B-C}p\ ,$ (17)

but with the final condition (9). The PMP control Hamiltonian reads

$\mathbb{H}=\langle p,\mathcal{K}_{C}\rho\rangle+s(t)\langle
p,\mathcal{K}_{B-C}\rho\rangle\ .$ (18)

###  IV.1 The “bang-anneal-bang” protocol is optimal

Applying Theorem 1 to the anti-Hermitian superoperator $\mathcal{L}$ of Eq.
(16), we obtain the following extension of the result of Ref. Brady  _et al._
(2021) to the density matrix setting:

######  Theorem 2.

(i) Assume $s^{*}\in[0,1]$ is the optimal control in an interval $[0,t_{f}]$
minimizing the cost (7) for Eq. (15). Then there exists a nonzero Hermitian
matrix solution of Eq. (17) with terminal condition (9) such that $s^{*}\equiv
0$ on intervals where $\langle p,\mathcal{K}_{B-C}\rho\rangle<0$, and
$s^{*}\equiv 1$ on intervals where $\langle p,\mathcal{K}_{B-C}\rho\rangle>0$.
On all other intervals $\langle p,\mathcal{K}_{B-C}\rho\rangle\equiv 0$ (these
are called singular arcs).

(ii) Assume furthermore that the constraint on the final time $t_{f}$ is
active (so that $\lambda>0$). Then $s^{*}(t)=1$ for
$t\in(t_{f}-\epsilon,t_{f}]$ for some $\epsilon>0$. Moreover, if the initial
condition $\rho_{0}$ commutes with the driver Hamiltonian $B$, i.e.,
$\mathcal{K}_{B}\rho_{0}=0$, one also has $s^{*}(t)=0$ for
$t\in[0,\epsilon^{\prime})$ for some $\epsilon^{\prime}>0$.

Before proving this theorem we offer a few remarks.

  * •

Part (i) implies that the optimal control is, in general, an alternation of
“bang” (nonsingular) arcs and “anneal” (singular) arcs where $\langle
p,\mathcal{K}_{B}\rho\rangle=\langle p,\mathcal{K}_{C}\rho\rangle$. Using the
PMP, this is an immediate consequence of the fact that the control enters
linearly in the equation and it is coupled to the superoperator
$\mathcal{K}_{B-C}$. _The latter is what is “special” about the quantum
annealing problem._

  * •

Part (ii) implies that under the assumption of an active time constraint and
for a particular initial condition, the optimal control starts and ends with
nonsingular arcs. In particular, it starts with an arc $s\equiv 0$ and ends
with an arc $s\equiv 1$.

  * •

In practice, whether there are additional nonsingular arcs in the middle is
problem dependent, and there is numerical evidence that such optimal scenarios
do indeed exist Brady  _et al._ (2021), but such nonsingular arcs do not exist
in the single qubit example discussed in Subsection IV.3.

  * •

Assuming that $\mathcal{K}_{B}\rho_{0}=0$ one can prove $\lambda\geq 0$ (see
also Appendix B, Proposition 1), but the condition $\lambda>0$ is more subtle
and it must be assumed independently. We return to this point in the next
subsection.

###### Proof.

Part (i): Eq. (18) states that the PMP control Hamiltonian $\mathbb{H}$
depends on the control only via the term $s(t)\langle
p,\mathcal{K}_{B-C}\rho\rangle$. If $\langle
p,\mathcal{K}_{B-C}\rho\rangle<0$, then to maximize this term as per Eq. (11)
subject to the constraint that $s(t)\in[0,1]$, clearly we must set
$s^{\ast}\equiv 0$. Likewise, if $\langle p,\mathcal{K}_{B-C}\rho\rangle>0$,
then to maximize this term subject to the same constraint requires
$s^{\ast}\equiv 1$. This is the case of nonsingular arcs. Conversely, if
$\langle p,\mathcal{K}_{B-C}\rho\rangle\equiv 0$ (a singular arc), then we
cannot conclude anything about the control from the PMP.

Part (ii): To investigate the form of the control at the end of the control
interval $[0,t_{f}]$, consider Eq. (12) with $\lambda>0$. Using Eq. (9) we
have $\langle
p_{f},\mathcal{K}_{C}\rho_{f}\rangle=\langle\mathcal{K}_{C}^{{\dagger}}p_{f},\rho_{f}\rangle=\langle\mathcal{K}_{C}C,\rho_{f}\rangle=0$.
This means that $\mathbb{H}(t_{f})=s(t_{f})\langle
p_{f},\mathcal{K}_{B}\rho_{f}\rangle=\lambda>0$, which in turn, since
$s\in[0,1]$, implies that $\langle p_{f},\mathcal{K}_{B}\rho_{f}\rangle>0$. By
continuity there must exist an interval $(t_{f}-\epsilon,t_{f}]$ (for some
$\epsilon>0$) such that $\langle p(t),\mathcal{K}_{B-C}\rho(t)\rangle>0$ for
$t\in(t_{f}-\epsilon,t_{f}]$, and in this interval we must have
$s^{\ast}(t)=1$ by (i).

The argument for the initial time is similar but instead of Eq. (9) it uses
the extra assumption $\mathcal{K}_{B}\rho_{0}=0$. Let us evaluate the control
Hamiltonian at $t=0$. Because of the assumption $\mathcal{K}_{B}\rho_{0}=0$ we
have $\mathbb{H}(t=0)=(1-s(0))\langle
p_{0},\mathcal{K}_{C}\rho_{0}\rangle=\lambda>0$. Since $s\in[0,1]$ this
implies that $\langle p_{0},\mathcal{K}_{C}\rho_{0}\rangle>0$ and $\langle
p_{0},\mathcal{K}_{-C}\rho_{0}\rangle<0$. By continuity there must exist an
$\epsilon^{\prime}>0$ such that, for $t\in[0,\epsilon^{\prime})$, $\langle
p(t),\mathcal{K}_{B-C}\rho(t)\rangle<0$ and in this interval we must have
$s^{\ast}(t)=0$ by (i). ∎

###  IV.2 The active constraint assumption and a sharpening of the results of
Ref. Brady  _et al._ (2021)

The condition $\lambda>0$ (that is, an active constraint on the final time
$t_{f}$) requires some extra discussion. It is a known fact in the geometric
theory of quantum control systems, and it follows as an application of general
results on control systems on Lie groups (see, e.g., Ref. (Jurdjevic and
Sussmann, 1972, Th. 7.2)), that there exists a _critical time_ $t_{c}$ such
that, the set of states reachable at time $t$, coincides for every $t\geq
t_{c}$. In other words, the reachable set does not grow past a certain time
$t_{c}$. Therefore, for every $t_{f}\geq t_{c}$ the time constraint is never
active. The minimum time $t_{\min}$ to reach the ground state of $C$ is $\leq
t_{c}$. If the final time $t_{f}$ is greater than or equal to $t_{\min}$, then
again the time constraint can never be active. In order to avoid this
situation, it was claimed in Ref. Brady  _et al._ (2021) that having
$t_{f}<t_{\min}$, is sufficient for having $\lambda>0$ in Eq. (12). Their
argument only uses $[\rho_{0},B]=0$. However, in Sec. IV.3 below we give an
example satisfying this assumption for which $\lambda=0$ for arbitrarily small
$t_{f}$. Thus, the assumption $t_{f}<t_{\min}$ is certainly necessary for
$\lambda>0$ but is in fact not sufficient. Rather, $\lambda>0$ is a feature of
the optimal trajectory rather than of the problem itself. This can be
explained more easily geometrically, as we now do.

The optimal cost at time $t_{f}$ is the minimum of a continuous function on
the reachable set of states $\mathfrak{R}_{t_{f}}$ (see Appendix B). It is
also known, under conditions that apply in our case, that the reachable set
$\mathfrak{R}_{t_{f}}$ varies continuously with $t_{f}$ Ayala  _et al._
(2017). We can map the space of Hermitian matrices $\rho$ diffeomorphically to
$\mathbb{R}^{n^{2}}$ (see Appendix C) and consider its reachable set there.
Since the cost function (7) is linear on this set, the minimum occurs on the
boundary. Therefore, the optimal trajectory is a curve starting from the
initial condition $\rho_{0}$ and ending on the boundary of
$\mathfrak{R}_{t_{f}}$. At the endpoint, the trajectory will have a tangent
vector which indicates its future direction. Now $\lambda>0$ if, going
(infinitesimally) in that direction combined with an increase in the size of
the reachable set $\mathfrak{R}_{t_{f}+\epsilon}$ for some small $\epsilon$,
will result in a reduced cost, and this is what we mean by the time constraint
being active. If $t_{f}$ is such that the reachable set does not increase at
$t_{f}$, for instance if $t_{f}\geq t_{c}$, then clearly this is not possible
and we must have $\lambda=0$. However, it is also possible that the reachable
set increases but not in a way to (strictly) decrease the cost, in particular
the portion of the boundary where we landed might not move at all, or it might
move but not in a direction that decreases the cost. This geometric discussion
is illustrated with figures in Appendix E.

The phenomenon that the optimal cost does not decrease with an increasing
final time $t_{f}$ may occur even though $t_{f}$ is arbitrarily small. Let us
denote by $J_{\min}(t_{f})$ the minimum cost (7) as a function of $t_{f}$. The
example we provide below (Sec. IV.3) shows that, even assuming
$[\rho_{0},B]=0$, we can have $J_{\min}(t_{f})=J_{\min}(0)$ for
$t_{f}\in[0,\epsilon)$ and some $\epsilon>0$, that is, the cost cannot be
lowered for some time, independently of the control. However, under the
additional assumption that $\rho_{0}$ is the _nondegenerate ground state of
$B$_ this does not happen, and we have the following theorem which we prove in
Appendix F:

######  Theorem 3.

Assume that $\rho_{0}$ in Eq. (6) is the nondegenerate ground state of $B$.
Then there exists an $\epsilon>0$ such that, for every $t_{f}\in(0,\epsilon)$,
$J_{\min}(t_{f})<J_{\min}(0)=\operatorname{Tr}(C\rho_{0})$.

In other words, if we start from the nondegenerate ground state of $B$ we can
always decrease the cost for sufficiently small $t_{f}$. Note that this,
however, does not prove that $\lambda>0$. As we have explained, the condition
$\lambda>0$ is a condition about the optimal trajectory, and it is an open
problem to find sufficient conditions such that every optimal trajectory
satisfies the $\lambda>0$ requirement for sufficiently small $t_{f}$.

###  IV.3 Example: optimal control of a spin-$1/2$ particle

We now give an example showing that without the assumption that $\rho_{0}$ is
the nondegenerate ground state of $B$, the cost (7) cannot be lowered even for
arbitrarily small $t_{f}$’s.

Consider a spin-$1/2$ particle (qubit) in a magnetic field. The model is given
by Eq. (15) with $C=\frac{1}{2}\sigma^{z}$ and $B=\frac{1}{2}\sigma^{x}$. As
an orthonormal, Hermitian operator basis we choose
$F_{i}=\frac{1}{\sqrt{2}}\sigma_{i}$, where we denote the standard Pauli
matrices $\sigma^{x}\equiv\sigma_{1}$ etc., i.e.:

$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\quad\sigma_{2}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),\quad\sigma_{3}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)\ .$ (19)

They satisfy the su$(2)$ commutation relations

$\left[\sigma_{1},\sigma_{2}\right]=2i\sigma_{3},\quad\left[\sigma_{3},\sigma_{1}\right]=2i\sigma_{2},\quad\left[\sigma_{2},\sigma_{3}\right]=2i\sigma_{1}\
.$ (20)

We parametrize the density matrix as $\rho=\frac{1}{2}\left(\leavevmode{\rm
1\mkern-4.8muI}+\boldsymbol{v}\cdot\boldsymbol{\sigma}\right)$, where
$\boldsymbol{v}\in\mathbb{R}^{3}$ is the Bloch vector ($\|\boldsymbol{v}\|\leq
1$) and $\boldsymbol{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})^{T}$. The
Bloch vector satisfies Eq. (15) where, using
$\boldsymbol{\mathcal{K}}_{ij}=\operatorname{Tr}[F_{i}\boldsymbol{\mathcal{K}}(F_{j})]$
(see Appendix C), we have

$\boldsymbol{\mathcal{K}}_{B}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&-1\\\
0&1&0\end{array}\right),\quad\boldsymbol{\mathcal{K}}_{C}=\left(\begin{array}[]{ccc}0&-1&0\\\
1&0&0\\\ 0&0&0\end{array}\right)\ .$ (21)

Equivalently, the dynamics are given by the Bloch equation

$\displaystyle\dot{\boldsymbol{v}}$
$\displaystyle=\boldsymbol{\mathcal{M}}(s)\boldsymbol{v}$ (22a)
$\displaystyle\boldsymbol{\mathcal{M}}(s)$
$\displaystyle=\left(\begin{array}[]{ccc}0&-(1-s)&0\\\ 1-s&0&-s\\\
0&s&0\end{array}\right)\ .$ (22e)

Geometrically, $\boldsymbol{\mathcal{K}}_{C}$ is the infinitesimal generator
of a counterclockwise rotation about the $v_{3}$ axis, while
$\boldsymbol{\mathcal{K}}_{B}$ is the infinitesimal generator of a
counterclockwise rotation about the $v_{1}$ axis. For $s\in(0,1)$,
$\boldsymbol{\mathcal{M}}(s)=(1-s)\boldsymbol{\mathcal{K}}_{C}+s\boldsymbol{\mathcal{K}}_{B}$
generates a counterclockwise rotation about an intermediate axis in the
$(v_{1},v_{3})$ plane. The cost (7) in this case becomes $J=\frac{1}{2}v_{3}$,
i.e., it corresponds to the minimization of the $v_{3}$ component.
Furthermore, let us assume for simplicity that the initial state $\rho_{0}$ is
pure ($\|\boldsymbol{v}\|=1$). There are only two such states compatible with
the condition $[B,\rho_{0}]=0$ (equivalently:
$\boldsymbol{\mathcal{K}}_{B}\boldsymbol{v}_{0}=\boldsymbol{0}$): the
$\sigma_{1}$ eigenstates, i.e., $\boldsymbol{v}_{0}=(\pm 1,0,0)^{T}$.

Now, if the initial state is $\boldsymbol{v}_{0}=(1,0,0)^{T}$, i.e., the
excited state of $B$, then for sufficiently small $t$ we have $v_{3}(t)\geq 0$
independently of the control $s\in[0,1]$ (see Appendix D). Therefore, an
optimal control in $[0,t_{f}]$ for $t_{f}$ small will be $s\equiv 0$ (which
will keep the value of $v_{3}$ at zero). The value of the minimum cost
$J_{\min}(t_{f})$ is equal to $J(0)$ for any arbitrarily small $t_{f}$. The
constraint on the final time is not active here, _even for arbitrarily small
$t_{f}$_. As a consequence, in this case we cannot draw the conclusions of
Theorem 2 following from the assumption $\lambda>0$. On the other hand, for
$\boldsymbol{v}_{0}=(-1,0,0)^{T}$, which corresponds to the (nondegenerate)
_ground state_ , with sufficiently small $t_{f}$ we can lower the cost
according to Theorem 3. We prove in Appendix H that the optimal control in
this case is a simple bang-bang protocol:

######  Theorem 4.

The optimal control for the system of one spin-$1/2$ particle considered
above, starting from the ground state and minimizing the cost
$\operatorname{Tr}(\sigma_{3}\rho)$ in time $t_{f}<\pi$, is the sequence
$s^{\ast}\equiv 0$ for time $\frac{t_{f}}{2}$ followed by $s^{\ast}\equiv 1$
for time $\frac{t_{f}}{2}$ (see Fig. 1).

Here $t_{c}=\pi$, i.e., if $t_{f}\geq\pi$ one trivially finds the ground state
exactly (by a $\pi/2$ rotation from the $-1$ eigenstate of $\sigma^{x}$ to the
$-1$ eigenstate of $\sigma^{y}$, followed by another $\pi/2$ rotation to the
$-1$ eigenstate of $\sigma^{z}$) and one cannot do better by increasing
$t_{f}$. This optimal bang-bang schedule result for a single spin-$1/2$ joins
previous such results for systems as diverse as pairs of one-dimensional
quasicondensates Rahmani  _et al._ (2013) or “gmon” qubits Bao  _et al._
(2018), as well as braiding of Majorana zero modes Karzig  _et al._ (2015).

Figure 1: (Color online) Single qubit case. Optimal trajectory on the Bloch
sphere with initial condition $\boldsymbol{v}_{0}=(-1,0,0)^{T}$. Here
$t_{f}=0.95\pi$ so the ground state of $C$, corresponding to the point
$(0,0,-1)$, is not reached exactly.

##  V The open system case

In this section we generalize the results for closed systems to the open
system setting. We consider two different approaches: an approximation-free
treatment of a system + environment where both are finite-dimensional, and a
master equation approach subject to a Markovian approximation, which applies
for infinite-dimensional environments E.B. Davies (1976); Alicki and Lendi
(2007); Breuer and Petruccione (2002); Rivas and Huelga (2012); Lidar (2019).
We show that under a number of additional assumptions, we can (partially)
recover the results from the closed system setting, but that the bangs
characterizing the latter are not a particularly robust feature in the open
system setting.

###  V.1 Optimal control for the Liouville-von Neumann equation

One approach for extending the closed system results of the previous section
to open systems is to consider the _full dynamics_ of a jointly evolving
system + environment. In this case $\rho$ in Eq. (6) is the density matrix of
the system and the environment, with an initial condition which is usually
taken to be of the factorized form $\rho(0):=\rho_{0}\otimes\rho_{E}$, where
now $\rho_{0}$ refers to the initial state of the system only. The Liouville-
von Neumann equation is [extending Eq. (15)]:

$\dot{\rho}=\mathcal{K}_{H_{\mathrm{tot}}}\rho\
,\quad\rho(0)=\rho_{0}\otimes\rho_{E}\ ,$ (23)

where $\mathcal{K}$ is defined in Eq. (13), with the total Hamiltonian

$H_{\mathrm{tot}}=H_{S}\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}+H_{I}+\leavevmode{\rm 1\mkern-4.8muI}_{S}\otimes H_{E}\ .$
(24)

Here $H_{I}$ is the interaction between the system and environment, $H_{E}$
generates the dynamics of the environment, while the system Hamiltonian, as
before, contains the controllable part:

$H_{S}(t)=C+s(t)\left(B-C\right)\ .$ (25)

Finally the cost is given by

$J=\operatorname{Tr}[\left(C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}\right)\rho]\ .$ (26)

Theorem 1 holds with $\mathcal{L}=\mathcal{K}_{H_{\mathrm{tot}}}$ and the PMP
control Hamiltonian has the form

$\mathbb{H}=\langle p,\mathcal{K}_{H_{S}\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho\rangle+\langle
p,\mathcal{K}_{H_{I}}\rho\rangle+\langle p,\mathcal{K}_{\leavevmode{\rm
1\mkern-4.8muI}_{S}\otimes H_{\mathrm{E}}}\rho\rangle\ .$ (27)

The treatment of Sec. IV applies, mutatis mutandis. In particular, condition
(9) is replaced by

$p(t_{f})=-C\otimes\leavevmode{\rm 1\mkern-4.8muI}_{E}\ .$ (28)

Remarkably, no additional modifications of the statement of the PMP Theorem 1
are needed. Moreover, it is clear from Eqs. (25) and (28) that once again the
control enters $\mathbb{H}$ only via the term $s(t)\langle
p,\mathcal{K}_{(B-C)\otimes\leavevmode{\rm 1\mkern-4.8muI}_{E}}\rho\rangle$,
so that the proof of Part (i) of Theorem 2 applies without any change. This
shows that:

######  Corollary 1.

In the general open system setting of Eq. (23) with the cost (26), the optimal
control $s(t)$ is an alternation of bang arcs where $s\equiv 0$ (when $\langle
p,\mathcal{K}_{(B-C)\otimes\leavevmode{\rm 1\mkern-4.8muI}}\rho\rangle<0$) or
$s\equiv 1$ (when $\langle p,\mathcal{K}_{(B-C)\otimes\leavevmode{\rm
1\mkern-4.8muI}}\rho\rangle>0$), and singular arcs where $\langle
p,\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}}\rho\rangle\equiv\langle
p,\mathcal{K}_{B\otimes\leavevmode{\rm 1\mkern-4.8muI}}\rho\rangle$.

Let us consider the generalization of Part (ii) of Theorem 2, which addresses
the characterization of the control function at the beginning and at the end.
We first consider the final arc. We have the following:

######  Theorem 5.

Assume $s^{*}\in[0,1]$ is the optimal control in an interval $[0,t_{f}]$
minimizing the cost (26) for Eq. (23). Assume furthermore that
$\left[H_{I},C\otimes\leavevmode{\rm 1\mkern-4.8muI}_{E}\right]=0$ and that
the final time constraint is active, i.e., $\lambda>0$. Then $s^{*}(t)=1$ for
$t\in(t_{f}-\epsilon,t_{f}]$ for some $\epsilon>0$.

Note that the assumption $\left[H_{I},C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}\right]=0$ implies that the choice of control $s\equiv 0$
leaves the cost $J=\langle C\otimes\leavevmode{\rm
1\mkern-4.8muI},\rho\rangle$ unchanged since in this case $H_{\mathrm{tot}}$
commutes with $C\otimes\leavevmode{\rm 1\mkern-4.8muI}_{E}$.

###### Proof.

Let us compute the PMP Hamiltonian at $t=t_{f}$. Note first that it follows
from Eq. (28) that

$\displaystyle\langle p_{f},\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho_{f}\rangle$
$\displaystyle=\langle\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}^{{\dagger}}p_{f},\rho_{f}\rangle$
$\displaystyle=\langle\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E},\rho_{f}\rangle=0\ .$ (29)

The other two terms of the PMP control Hamiltonian Eq. (27) also vanish at
$t=t_{f}$: using the same calculation as in Eq. (29) the second term vanishes
because of the assumption $\left[H_{I},C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}\right]=0$, and the third term does as well because
$[\leavevmode{\rm 1\mkern-4.8muI}_{S}\otimes H_{E},C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}]=0$. So we obtain $\mathbb{H}=s(t_{f})\langle
p_{f},\mathcal{K}_{B\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho_{f}\rangle=\lambda>0$. This implies that $\langle
p_{f},\mathcal{K}_{B\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho_{f}\rangle>0$ and by continuity there must exist an
interval $(t_{f}-\epsilon,t_{f}]$ for some $\epsilon>0$ such that $\langle
p(t),\mathcal{K}_{(B-C)\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho(t)\rangle>0$ for $t\in(t_{f}-\epsilon,t_{f}]$.
Finally we must have $s^{\ast}\equiv 1$ in this interval by Corollary 1. ∎

Regarding the arc at the beginning we have instead:

######  Theorem 6.

Assume $s\in[0,1]$ is the optimal control in an interval $[0,t_{f}]$
minimizing the cost Eq. (26) for Eq. (23). Assume that
$\left[\rho_{E},{H}_{E}\right]=0$ and that
$\left[H_{I},\rho_{0}\otimes\rho_{E}\right]=0$. Assume furthermore that
$\left[B,\rho_{0}\right]=0$ and that the final time constraint is active
$\lambda>0$. Then the control satisfies $s^{\ast}(t)=0$ for $t\in[0,\epsilon)$
for some $\epsilon>0$.

Note that the assumption $\left[H_{I},\rho_{0}\otimes\rho_{E}\right]=0$
implies that the interaction alone does not modify the initial state of the
system.

###### Proof.

We abbreviate the proof since it is very similar to the ones we presented
above in more detail. Using the assumption
$\left[H_{I},\rho_{0}\otimes\rho_{E}\right]=0$ and
$\left[B,\rho_{0}\right]=0$, evaluating $\mathbb{H}(t=0)$ we obtain
$\mathbb{H}(0)=\left(1-s(0)\right)\langle
p(0),\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho(0)\rangle=\lambda>0$. This implies that $\langle
p(0),\mathcal{K}_{C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho(0)\rangle>0$ or equivalently that $\langle
p(0),\mathcal{K}_{-C\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho(0)\rangle<0$. By continuity this in turn implies that
there exist an interval $[0,\epsilon)$ for some $\epsilon>0$ such that
$\langle p(t),\mathcal{K}_{(B-C)\otimes\leavevmode{\rm
1\mkern-4.8muI}_{E}}\rho(t)\rangle<0$ for $t\in[0,\epsilon)$. Finally we must
have $s^{\ast}\equiv 0$ in this interval by Corollary 1. ∎

We comment on the implications of the additional assumptions used in these
theorems in Sec. VII.

###  V.2 Optimal control for quantum master equation dynamics

The treatment of the open system case in the previous subsection did not
involve any approximations. On the other hand, we tacitly assumed that the
environment is _finite dimensional_. This was helpful since all the results on
optimal control which we have elaborated upon in Sec. IV and used so far, are
classically stated and proved for finite dimensional systems. Extending such
results, in particular concerning the PMP and the topology and continuity of
the reachable sets for infinite dimensional systems, is possible and is a
current area of research in control theory (see, e.g., Refs. Fiacca  _et al._
(1998); Kipka and Ledyaev (2015)), although the results in this area become
considerably more technical. An alternative we discuss in this subsection is
to replace the Liouville-von Neumann equation (23) with an approximate quantum
master equation. This can be viewed as an investigation of the result of Sec.
V.1 when the environment dimension is sent to infinity.

Without loss of generality we write $H_{I}=\sum_{\alpha}S_{\alpha}\otimes
E_{\alpha}$, where $S_{\alpha}=S_{\alpha}^{\dagger}$ and
$E_{\alpha}=E_{\alpha}^{\dagger}$ $\forall\alpha$. The goal is now to find a
master equation for the dynamics of the system density matrix in the case of a
time-dependent system Hamiltonian. After the Born approximation and tracing
out the environment, one arrives at a time dependent Redfield master equation
(see, e.g., the Schrödinger picture Redfield master equation (SPRME) derived
in Ref. Campos Venuti and Lidar (2018)). From this point there are multiple
ways to proceed, e.g., by introducing an additional adiabatic approximation or
an additional Markovian approximation, or both. These different paths, and
exactly how they are taken, lead to a plethora of different master equations
Childs  _et al._ (2001); Alicki  _et al._ (2006); Albash  _et al._ (2012);
Campos Venuti and Lidar (2018); Mozgunov and Lidar (2020); Yamaguchi  _et al._
(2017); Dann  _et al._ (2018); Nathan and Rudner (2020); Davidović (2020);
Winczewski  _et al._ (2021). We next focus on two representative cases of
master equations derived from first principles.

####  V.2.1 Adiabatic Redfield Master Equation

The Adiabatic Redfield Master Equation (ARME) is derived in Ref. Campos Venuti
and Lidar (2018). It results from assuming that $t_{f}\gg\tau_{B}$, where
$\tau_{B}$ is the environment time scale, and the adiabatic approximation
$\mathcal{T}\exp\left[-i\int_{t-r}^{t}H_{S}(t^{\prime})dt^{\prime}\right]\approx
e^{-irH_{S}(t)}$, dropping a correction of $O\left((r/t_{f})^{2}\right)$. The
ARME has the form of Eq. (6) with a time-dependent Redfield generator
$\mathcal{L}$ given by

$\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{K}_{H_{S}}+\mathcal{D}$
(30a) $\displaystyle\mathcal{D}\rho$
$\displaystyle=\sum_{\alpha\beta}\int_{0}^{t_{\max}}dr\
G_{\alpha\beta}(r)\,\left[S_{\beta}(-r)\rho,S_{\alpha}\right]+\mathrm{h.c.},$
(30b)

where the system Hamiltonian is as in Eq. (25). $G_{\alpha\beta}(t)$ is the
environment correlation function

$G_{\alpha\beta}(t)=\langle
E_{\alpha}(t)E_{\beta}(0)\rangle=\operatorname{Tr}[E_{\alpha}(t)E_{\beta}(0)\rho_{E}]\
,$ (31)

where $\left\langle X\right\rangle$ denotes the environmental thermal average
of $X$. When $G_{\alpha\beta}(r)$ decays exponentially, the relative error of
the resulting dynamics due to the adiabatic approximation above is
$O\left((\tau_{B}/t_{f})^{2}\right)$. Finally,

$S_{\beta}(-r)=e^{-irH_{S}(s)}S_{\beta}e^{irH_{S}(s)}\ .$ (32)

The parameter $t_{\max}$ can either be set to $t_{f}$ or infinity on account
of the fact that the environment correlation function decays very rapidly. The
ARME is not in Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form Gorini  _et
al._ (1976); Lindblad (1976); Chruściński and Pascazio (2017), hence does not
generate a completely positive map. However, it generates non-Markovian
dynamics, hence has a wider range of applicability than Markovian master
equations, within its range of applicability Mozgunov and Lidar (2020).

Crucially, the generator in Eq. (30) with Eqs. (25), (31) and (32) _depends on
time only through the control function $s(t)$_. This implies that the PMP
control Hamiltonian $\mathbb{H}=\langle p,\mathcal{L}\rho\rangle$ is constant
and hence the PMP in the form of Theorem 1 is directly applicable.666It can be
shown that the conditions of Filippov’s theorem (see Ref. (Fleming and Rishel,
1975, Th. 2.1)) for the existence of the optimal control solution are
satisfied also in this case. However, the control now enters nonlinearly in
$\mathcal{D}$, in particular in an exponential through Eq. (32). As a
consequence it is not possible to derive the form of the control on the
nonsingular arcs, or even to determine simple equations for the appearance of
singular arcs. One can ask, however, what remains of the results of the
previous subsection. We do not have an analog of Theorem 6 for the initial
arc. However, if we again make the assumption of Theorem 5 that
$\left[S_{\alpha},C\right]=0$ $\forall\alpha$, then the analog of this theorem
for the final arc holds, even when the environment is infinite-dimensional,
and under the approximations used to derive Eq. (30). However, instead of an
arc, we obtain a bang only at a point:

######  Theorem 7.

Assume $s^{*}\in[0,1]$ is the optimal control in an interval $[0,t_{f}]$
minimizing the cost Eq. (7) for Eq. (6) with $\mathcal{L}$ given by Eq. (30).
Assume further that $\left[S_{\alpha},C\right]=0$ $\forall\alpha$ and that the
final time constraint is active ($\lambda>0$). Then the optimal control
satisfies $s^{*}(t_{f})=1$.

###### Proof.

It is convenient to write the Redfield dissipator as

$\displaystyle\mathcal{D}\rho$
$\displaystyle=\sum_{\alpha\beta}\left(\big{[}W_{\alpha\beta}\rho,S_{\alpha}\big{]}+\big{[}S_{\alpha},\rho
W_{\alpha\beta}^{\dagger}\big{]}\right)$ (33a) $\displaystyle
W_{\alpha\beta}(t)$ $\displaystyle=\int_{0}^{t_{\max}}dr\
G_{\alpha\beta}(r)\,S_{\beta}(-r)\ .$ (33b)

Using Eq. (33a) one obtains, for the adjoint of $\mathcal{D}$:

$\mathcal{D}^{{\dagger}}X=\sum_{\alpha\beta}\left(W_{\alpha\beta}^{\dagger}\left[X,S_{\alpha}\right]+\left[S_{\alpha},X\right]W_{\alpha\beta}\right)$
(34)

(see Appendix D). From the above expression and the assumption
$\left[S_{\alpha},C\right]=0,\,\forall\alpha$ we obtain
$\mathcal{D}^{{\dagger}}C=0$. Using $p_{f}=-C$ we have $\langle
p_{f},\mathcal{D}\rho_{f}\rangle=-\langle\mathcal{D}^{{\dagger}}C,\rho_{f}\rangle=0$.
Evaluating the PMP control Hamiltonian at the final time we obtain
$\mathbb{H}(t_{f})=s(t_{f})\langle
p_{f},\mathcal{K}_{B}\rho_{f}\rangle=\lambda>0$. This implies that $\langle
p_{f},\mathcal{K}_{B}\rho_{f}\rangle>0$, and so $s(t_{f})=1$ from the maximum
principle. ∎

Since, as argued in Sec. VI.2, the size of the bang intervals is generically
expected to shrink when the total system size increases, this result can be
seen as a generalization of Theorem 5 when the environment dimension is sent
to infinity and the bang interval at the end shrinks to a point.

The implications of the additional assumptions used in Theorem 7 are discussed
in Sec. VII.

####  V.2.2 Markovian, completely positive master equations

The most significant drawback of the ARME is the violation of complete
positivity, which means that the density matrix can develop unphysical,
negative eigenvalues. Hence we also consider Markovian, completely positive
master equations. There are a variety of such master equations derived from
first principles under different assumptions. However, in most cases the
generator $\mathcal{L}$ is explicitly time-dependent (e.g., the coarse-grained
master equation (CGME) (Mozgunov and Lidar, 2020, Eq. (22)), the master
equation of Ref. (Yamaguchi  _et al._ , 2017, Eq. (21)), the non-adiabatic
master equation (NAME) (Dann  _et al._ , 2018, Eq. (16)), and the universal
Lindblad equation (ULE) (Nathan and Rudner, 2020, Eq. (27))) and hence we
cannot apply Theorem 1.

In this subsection we give an example of a Markovian master equation derived
from first principles where, like in the ARME case, the generator
$\mathcal{L}$ depends on time only through the schedule $s(t)$. In this case
the PMP can be applied in the simplified form described in Theorem 1.

Consider the “geometric-arithmetic master equation” (GAME) (Davidović, 2020,
Eq. (46)), which is claimed there to have a higher degree of accuracy than all
the previous Markovian master equations. In the adiabatic limit it has the
Schrödinger picture form

$\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{K}_{H_{S}}+\mathcal{D}$
(35a) $\displaystyle\mathcal{D}\rho$
$\displaystyle=\sum_{j}\left([L_{\alpha}(s)\rho,L_{\alpha}^{\dagger}(s)]+[L_{\alpha}(s),\rho
L_{\alpha}^{\dagger}(s)]\right)\ ,$ (35b)

where $L_{\alpha}(s)=S_{\alpha}\circ\sqrt{\gamma(s)}$, the circle denotes the
Hadamard (element-wise) product, $\gamma$ is the spectral density matrix
[Fourier transform of the environment correlation function (31)] whose
elements $\gamma_{nm}:=\gamma(\omega_{nm}(s))$ depend on the instantaneous
Bohr frequencies $\omega_{nm}(s)=E_{n}(s)-E_{m}(s)$, where
$H_{S}(s)|n(s)\rangle=E_{n}(s)|n(s)\rangle$, and the dependence on time is
only through the schedule $s(t)$.777In writing these expressions we have
adapted the results of Ref. Davidović (2020) to the adiabatic limit by using
the instantaneous eigenbasis of $H_{S}$, and also performed the
“$\int^{t}\to\int^{\infty}$ approximation” (otherwise $\gamma$ would have a
$t$ dependence, which would prevent us from being able to apply Theorem 1);
see Ref. Davidović (2020) for complete details. The adjoint dissipator is now:

$\displaystyle\mathcal{D}^{{\dagger}}X$
$\displaystyle=\sum_{\alpha}\big{(}(S_{\alpha}\circ\sqrt{\gamma})^{\dagger}\left[X,S_{\alpha}\circ\sqrt{\gamma}\right]$
$\displaystyle\qquad+[(S_{\alpha}\circ\sqrt{\gamma})^{\dagger},X]S_{\alpha}\circ\sqrt{\gamma}\big{)}\
.$ (36)

Unfortunately, since
$[S_{\alpha},C]=0\not\Rightarrow[S_{\alpha}\circ\sqrt{\gamma},C]=0$, this
means that even if $[S_{\alpha},C]=0$, we do _not_ obtain
$\mathcal{D}^{\dagger}C=0$ as in the ARME case, and hence the proof of Theorem
7 does not carry through.888Of course $[S_{\alpha},C]=0\
\Rightarrow[S_{\alpha}\circ\sqrt{\gamma},C]=0$ when $\gamma_{nm}=c$ $\forall
n,m$, but this is a highly nongeneric scenario.

While these arguments are not a proof that in general Markovian dynamics do
not admit $s(t_{f})=1$ as an optimal control solution, we conjecture that in
fact, they do not. It thus appears that the “counterintuitive” appearance of
the driver Hamiltonian at the end of the control interval is not a feature of
the optimal schedule in the Markovian limit of open quantum systems. We
revisit this point in Sec. VII.

##  VI Switching operator and analysis of the optimal control

In order to study the qualitative behavior of the optimal control law, in
particular its switching properties and the existence and nature of the
singular arcs, it is convenient to introduce one more operator, besides the
state $\rho$ and the co-state $p$, which we call the _switching operator_. The
switching operator determines the behavior of the optimal control, i.e., the
points where there is a switch between $s=0$ and $s=1$, and where there is a
singular arc. For clarity we focus on the closed system case of Sec. IV but
our definitions and treatment naturally extend with a change of notation to
the open system case of Sec. V.1.

###  VI.1 Switching equation

The switching operator $S$ is the Hermitian operator defined as

$S:=i\left[p,\rho\right]\ .$ (37)

In the closed system case the Liouvillian has the form
$\mathcal{L}=\mathcal{K}_{H}$, with the system Hamiltonian of Eq. (2). One has
the following property:

$\mathcal{K}_{H}\left(\left[X,Y\right]\right)=\left[\mathcal{K}_{H}(X),Y\right]+\left[X,\mathcal{K}_{H}(Y)\right]$
(38)

valid for any operators $X,Y,H$. Now, differentiating Eq. (37), we obtain

$\displaystyle\dot{S}$
$\displaystyle=i\left[\dot{p},\rho\right]+i\left[p,\dot{\rho}\right]$ (39a)
$\displaystyle=i\left[\mathcal{K}_{H}p,\rho\right]+i\left[p,\mathcal{K}_{H}\rho\right]$
(39b)
$\displaystyle=i\mathcal{K}_{H}\left(\left[p,\rho\right]\right)=\mathcal{K}_{H}S\
,$ (39c)

where in the third equality we used Eq. (38). Thus, $S$ satisfies the same
equation as $\rho$ and $p$. To understand why $S$ determines the optimal
control switching times, let us define

$\displaystyle x_{C}$ $\displaystyle:=\langle
p,\mathcal{K}_{C}\rho\rangle=-i\operatorname{Tr}\left(p\left[C,\rho\right]\right)$
(40a)
$\displaystyle=-i\operatorname{Tr}\left(C\left[\rho,p\right]\right)=\operatorname{Tr}\left(CS\right)=\langle
C,S\rangle\ ,$ (40b)

and similarly

$x_{B}:=\langle p,\mathcal{K}_{B}\rho\rangle=\langle B,S\rangle\ ,$ (41)

so that the PMP control Hamiltonian Eq. (18) can be written in a form closely
resembling the Hamiltonian (2b):

$\mathbb{H}=x_{C}+s\left(x_{B}-x_{C}\right)\ .$ (42)

The quantity $x_{B}-x_{C}=\langle B-C,S\rangle$, that is, the orthogonal
component of $S$ along $B-C$, regulates the switches of the candidate optimal
control. By the same PMP argument we have used repeatedly in our proofs, when
$x_{B}-x_{C}<0$ we have $s\equiv 0$, and when $x_{B}-x_{C}>0$ we have $s\equiv
1$. The switch occurs when $x_{B}-x_{C}=0$, while a singular arc occurs when
$x_{B}-x_{C}\equiv 0$ for an interval of positive measure.

The initial condition $S_{0}$ of the switching operator determines the optimal
control candidate $s$ uniquely. This initial condition is not completely
arbitrary. In particular, under the assumptions of Theorem 2 the following
holds. Since $s(0)=0$, we have $x_{C}(0)=\langle C,S_{0}\rangle=\lambda$.
Furthermore, since $[B,\rho_{0}]=0$, we have

$\displaystyle x_{B}(0)$ $\displaystyle=\langle
B,S_{0}\rangle=i\operatorname{Tr}(B[p_{0},\rho_{0}])$
$\displaystyle=i\operatorname{Tr}\left(p_{0}[\rho_{0},B]\right)=0\ .$ (43)

At the final time $t_{f}$, since $s(t_{f})=1$, we have from Eq. (42)
$x_{B}(t_{f})=\lambda$, while using Eq. (9) we have

$\displaystyle x_{C}(t_{f})$ $\displaystyle=\langle
C,S_{f}\rangle=i\operatorname{Tr}\left(C[p_{f},\rho_{f}]\right)$
$\displaystyle=i\operatorname{Tr}\left(\rho_{f}[C,p_{f}]\right)=0\ .$ (44)

Thus, at $t=0$ we have $x_{B}(0)=0$ and in the initial arc $s\equiv 0$ and
$x_{C}\equiv\lambda$ by Eq. (42). The next arc can be either nonsingular with
$s\equiv 1$ ($x_{B}>x_{C}$) or a singular arc with $x_{B}\equiv x_{C}$. Either
way, at the switching point we must have $x_{B}=\lambda$, hence we reach the
point $(x_{C},x_{B})=(\lambda,\lambda)$ at the end of the first arc. When
$(x_{C},x_{B})=(\lambda,\lambda)$ there is either a switch to an arc with
$s\equiv 1$, or a return to $s\equiv 0$, or a singular arc where
$(x_{C},x_{B})\equiv(\lambda,\lambda)$. On this arc $s$ is unspecified, but
nonetheless certain equations need to be satisfied and they can be used to
obtain information on the dynamics on such singular arcs (see Appendix G).
Note that every switching event, whether from a bang arc to a singular arc or
v.v.., or from a bang arc to another bang arc, happens at
$(x_{C},x_{B})=(\lambda,\lambda)$. When $s\equiv 1$, from Eq. (42) we have
that $x_{B}$ is constant, while $x_{C}$ is allowed to change. Therefore the
optimal control can be described schematically as in Fig. 2 in the
$(x_{C},x_{B})$ plane.

Figure 2: Switching diagram for optimal control candidates. Starting from the
point $(\lambda,0)$ the dynamics stay on the line $x_{C}=\lambda$. In
principle even negative values of $x_{B}$ can be attained, but eventually the
point $(\lambda,\lambda)$ must be reached, followed by an alternation of
vertical and horizontal lines (where in principle one can also extend to
negative values of $x_{C}$) through $(\lambda,\lambda)$ with time intervals
where the dynamics do not move from $(\lambda,\lambda)$, representing the
singular arcs. In the last interval, the dynamics follow a horizontal line
(with $s\equiv 1$) reaching the point $(x_{C},x_{B})=(0,\lambda)$.

In the case where the final time is not active, i.e., when we do not have the
guarantee $\lambda>0$ in Eq. (42), the above reasoning can still be applied to
conclude $x_{B}(0)=0$ and $x_{C}(t_{f})=0$. Furthermore, if $\lambda=0$, from
Eq. (42), we obtain $(1-s(0))x_{C}(0)=0$ and $s(t_{f})x_{B}(t_{f})=0$, in
addition to the maximization condition Eq. (11).

###  VI.2 Shortening of nonsingular arcs

####  VI.2.1 Dependence on $n$

Here we give a heuristic argument that explains why the terminal arcs
($s\equiv 0$ and $s\equiv 1$) for the optimal control become shorter as the
number $n$ of spins (or qubits) increases. In fact the heuristic holds also
for intermediate arcs taking place anywhere along the optimal trajectory.

Consider a bang arc where $s(t)=0$ for $t\in\left[t_{0},t_{1}\right]$.
Equation (39) for the switching operator in this region is
$\dot{S}=\mathcal{K}_{C}S$ with the initial condition $S(t_{0})=S_{0}$ for
some $S_{0}$. The solution in this interval is
$S(t)=e^{-i(t-t_{0})C}S_{0}e^{i(t-t_{0})C}.$ The coordinate $x_{C}$ equals
$\lambda$ in the interval: $x_{C}=\lambda=\langle C,S(t)\rangle=\langle
C,S_{0}\rangle$. The switching happens when $x_{B}=\lambda$, i.e., at the
first solution $t_{1}$ of $\langle B,S(t_{1})\rangle=\langle C,S_{0}\rangle$.
More explicitly, the interval of the bang arc $\Delta t:=t_{1}-t_{0}$ is given
by the first solution of

$x_{B}(\Delta t):=\operatorname{Tr}\left(Be^{-i\Delta tC}S_{0}e^{i\Delta
tC}\right)=\operatorname{Tr}\left(CS_{0}\right)\ .$ (45)

Using the spectral resolution $C=\sum_{k}E_{k}\Pi_{k}$ (with eigenvalues
$E_{k}$ and eigenprojectors $\Pi_{k}$), the left-hand hand side of Eq. (45)
can be written as

$x_{B}(\Delta t)=\sum_{kl}M_{kl}e^{-i\Delta t\omega_{kl}}\ ,$ (46)

with amplitudes $M_{kl}:=\operatorname{Tr}\left(B\Pi_{k}S_{0}\Pi_{l}\right)$
and Bohr frequencies $\omega_{kl}=E_{k}-E_{l}$. The function $x_{B}(\Delta t)$
is a real trigonometric polynomial with $O(d^{2})$ terms (where $d$ is the
Hilbert space dimension) starting from $\operatorname{Tr}(BS_{0})$ at $\Delta
t=0$, and reaching $\lambda>0$ at $\Delta t$. Now, as the number of qubits $n$
increases, both $d=2^{n}$ and the frequencies increase. Both of these facts
contribute to making $x_{B}(t)$ oscillate faster. As a consequence, the
solution $\Delta t$ of Eq. (45) tends to decrease with $n$. The same
considerations hold for the case of an $s\equiv 1$ arc with $B$ and $C$
interchanged. See Appendix I for additional comments. Note that this heuristic
applies both to the initial and final bang arcs, as well as to possible
intermediate bang arcs if they are present.

####  VI.2.2 Dependence on $t_{f}$

It was concluded in Ref. (Brady  _et al._ , 2021, Sec. S3) that “these bangs
should become smaller and smaller as $t_{f}$ is increased. Eventually in the
true $t_{f}\to\infty$ adiabatic limit, these bangs disappear recovering the
standard form expected for quantum adiabatic computing.” However, there is in
fact no guarantee that the optimal control coincides with the adiabatic path
even in this limit, since the adiabatic theorem provides a sufficient, but not
a necessary condition for convergence to the minimum of the cost function.
Indeed, it is easy to construct a counterexample, as we now do. First note
that, as we have seen, when $t_{f}<t_{c}$, the optimal schedule always starts
with a bang $s\equiv 0$ and ends with a bang $s\equiv 1$ (provided the initial
state commutes with $B$). We have not addressed the question of uniqueness of
this optimal schedule, which we leave for future work. However, when
$t_{f}>t_{c}$, the optimal schedule is certainly not unique, as one has the
possibility of “wasting time” by adding a bang $s\equiv 0$ at the end (thus
applying $C$ there), or by adding a bang $s\equiv 1$ at the beginning (thus
applying $B$ there), or both. The resulting schedules do not resemble the
smooth adiabatic schedule interpolating slowly from $s(0)=1$ to $s(t_{f})=0$.

##  VII Summary and Discussion

The quest to discover the optimal schedule for quantum optimization algorithms
such as QA and QAOA naturally leads to the use of optimal control theory via
Pontryagin’s principle. Previous work concluded that QAOA is optimal Yang  _et
al._ (2017), but a more careful analysis showed that in fact a hybrid bang-
anneal-bang protocol is generally optimal for closed systems when not enough
time is allowed for the desired state to be reached perfectly Brady  _et al._
(2021). Here we confirmed this result using a density matrix approach, which
both generalizes the analysis to mixed states and simplifies it since it makes
the cost-function linear in the state. We also showed that the assumption that
$t_{f}$ is smaller than the critical time $t_{c}$ needed to reach the ground
state of $C$ exactly, is necessary but not sufficient for the result of Ref.
Brady  _et al._ (2021), by giving a counterexample to the latter.

We introduced a switching operator and found its equation of motion, which
characterizes the points at which the optimal schedule switches between the
two different types of bang arcs and the anneal arc. In Theorem 4 we gave the
explicit optimal schedule for the example of a single spin-$1/2$ particle,
which consists of two bangs of equal duration.

Using the density matrix formulation we extended the theory to the open system
setting, both for the exact reduced system dynamics in the case of a finite-
dimensional environment, and under the approximation of dynamics governed by a
master equation due to coupling to an infinite-dimensional environment. We
proved that in the first setting, depending on additional assumptions
concerning the initial states of the system and the environment and their
interaction, either an anneal-bang (Theorem 5) or bang-anneal (Theorem 6)
schedule is optimal.

In the second setting (infinite-dimensional environment) we considered both an
adiabatic Redfield equation accounting for non-Markovian dynamics but without
a complete positivity guarantee, and a completely positive Markovian master
equation. In the former (Redfield) case we could only prove that the optimal
schedule terminates with the driver Hamiltonian, i.e., $s(t_{f})=1$ (Theorem
7). One could interpret this result as a manifestation of the phenomenon of
the shortening of the nonsingular arcs as the _total_ system (i.e., the
subsystem plus its environment) size increases. Indeed, in the Redfield case
the environment Hilbert space dimension is infinite, which is consistent with
the bang arc having shrunk down to a point. In the fully Markovian case, even
this last remnant of the bang-arc was not recovered, as we found no evidence
of natural conditions under which $s(t_{f})=1$ holds.

Let us now comment on the differences between these theorems and their closed
system counterpart, Part (ii) of Theorem 2. Regarding Theorem 5 concerning the
final arc, the main change is the addition of the assumption that the
interaction Hamiltonian commutes with the cost function, i.e.,
$\left[H_{I},C\otimes\leavevmode{\rm 1\mkern-4.8muI}_{E}\right]=0$. Writing
the interaction in the general form $H_{I}=\sum_{\alpha}S_{\alpha}\otimes
E_{\alpha}$, where $S_{\alpha}$ and $E_{\alpha}$ are system and environment
operators, respectively, the assumption is equivalent to $[S_{\alpha},C]=0$
$\forall\alpha$, which is the same assumption as in Theorem 7. Thus $C$ must
belong to the commutant of the algebra generated by the set
$\\{S_{\alpha}\\}$, i.e., $C\in\text{Alg}\\{S_{\alpha}\\}^{\prime}$.999The
commutant of an algebra $\mathcal{A}=\text{Alg}\\{A\\}$ is defined as the set
$\mathcal{A}^{\prime}:=\\{X|[X,A]=0\ \forall A\in\mathcal{A}\\}$. For example,
if $S_{\alpha}=\sum_{i=1}^{n}\sigma_{i}^{\alpha}$ for a system of $n$ qubits,
where $\alpha\in\\{x,y,z\\}$, i.e., the collective decoherence case Zanardi
(1998); Lidar  _et al._ (1998), then, if $C$ is at most a two-body
interaction, it follows that it must be of the Heisenberg interaction form:
$C=\sum_{ij}J_{ij}\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}$, where
$J_{ij}$ are constants Bacon  _et al._ (2000). Or, for a classical target
Hamiltonian arising in optimization such as the Ising-type Hamiltonian
mentioned in Sec. II, this means that the interaction must be of the pure-
dephasing type, i.e., $S_{\alpha}\propto\sigma^{z}$ (or products of
$\sigma^{z}$ over different qubits). This is a realistic model, e.g., for
superconducting qubits undergoing flux noise Kjaergaard  _et al._ (2020).

Regarding Theorem 6 concerning the initial arc, the main change is the
addition of the two assumptions that (i) the environment Hamiltonian commutes
with the environment’s initial state ($\left[{H}_{E},\rho_{E}\right]=0$), and
(ii) the interaction Hamiltonian commutes with the joint system-environment
initial state ($\left[H_{I},\rho_{0}\otimes\rho_{E}\right]=0$). The first of
these is natural and is known as the stationary environment assumption Breuer
and Petruccione (2002). It is satisfied, e.g., if the environment is in
thermal equilibrium, i.e., in the Gibbs state: $\rho_{E}\propto e^{-\beta
H_{E}}$, where $\beta$ is the inverse temperature. The second assumption means
that $\rho_{0}\in\text{Alg}\\{S_{\alpha}\\}^{\prime}$ and
$\rho_{E}\in\text{Alg}\\{E_{\alpha}\\}^{\prime}$. This assumption is the least
natural of the ones we have encountered so far. E.g., it is clearly violated
in the standard quantum annealing setting where $\rho_{0}$ is the ground state
of a transverse field $-\sum_{i}\sigma_{i}^{x}$ and
$S_{\alpha}\propto\sigma^{z}$. Even the condition
$\rho_{E}\in\text{Alg}\\{E_{\alpha}\\}^{\prime}$ is not very natural. For
example, for a bosonic environment one typically has $E_{\alpha}$ as the
position operator of an oscillator, while $H_{E}$ might be the number
operator, in which case the Gibbs state $\rho_{E}$ would not commute with
$E_{\alpha}$.

We note that while the conditions given in Theorem 6 are sufficient, we do not
know if they are necessary, which we thus leave as an open problem. We
conjecture that the initial bang does not appear as a feature of optimal
schedules for open systems coupled to an infinite-dimensional environment.
This state of affairs would be reminiscent of the existence of an arrow of
time for open systems, which breaks the symmetry between the initial and final
times (see, e.g., Ref. Campos Venuti and Lidar (2018) for a similar effect in
the pure QA setting). On the other hand, given the naturalness of the
sufficient conditions under which a final bang arc (Theorem 5) or a schedule
terminating with the driver Hamiltonian (Theorem 7) are optimal, such
schedules may find utility in the design of quantum algorithms for
optimization problems in the setting of open quantum systems. This is true in
particular for systems that are well described by the adiabatic Redfield
master equation, e.g., superconducting flux qubits used for quantum annealing
Harris  _et al._ (2010); Yan  _et al._ (2016); Smirnov and Amin (2018); Khezri
_et al._ (2021).

However, the conditions under which the adiabatic Redfield master equations
hold need not apply in general, e.g., for Hamiltonians that arise naturally in
systems such as Rydberg atoms or transmons, which have been used to
demonstrate QAOA Zhou  _et al._ (2020); Harrigan  _et al._ (2021). Especially
for quantum optical systems such as Rydberg atoms, the Markovian limit may be
more appropriate, and we have not found evidence of the optimality of an
initial or final bang arc in this limit, or even the optimality of
$s(t_{f})=1$.

While our analysis does not strictly rule out bang-type schedules for open
systems coupled to infinite-dimensional environments, we conjecture that they
are indeed not a feature of optimal schedules in this case, primarily due to
the shortening of arcs in the open system setting. If this could be confirmed,
it would mean that after all, continuous annealing-type schedules are optimal
for optimization purposes when using open quantum systems, which would have
implications for all NISQ-era optimization algorithms.

Finally, we remark that throughout this work we have used a simplified form of
the PMP as described in Theorem 1. A more general PMP for time-dependent
dynamics is described in Ref. Fleming and Rishel (1975). The main difference
is that $\mathbb{H}\equiv\lambda$ [Eq. (12)] is no longer valid in the given
form. This equation is, in fact, a special case of another one which contains
the derivative of the dynamics with respect to $t$, and this term vanishes
when the dynamics are not explicitly time-dependent, as in our case, where we
considered adiabatic time-dependent master equations. The time-dependence of
the dynamics in typical quantum master equations Childs  _et al._ (2001);
Alicki  _et al._ (2006); Albash  _et al._ (2012); Campos Venuti and Lidar
(2018); Mozgunov and Lidar (2020); Yamaguchi  _et al._ (2017); Dann  _et al._
(2018); Nathan and Rudner (2020); Davidović (2020); Winczewski  _et al._
(2021) is, however, different from the one of models usually encountered in
classical control theory, and therefore further study is required before the
PMP can be applied to a broader class of quantum master equations.

###### Acknowledgements.

LCV’s and DL’s research is based upon work (partially) supported by the Office
of the Director of National Intelligence (ODNI), Intelligence Advanced
Research Projects Activity (IARPA) and the Defense Advanced Research Projects
Agency (DARPA), via the U.S. Army Research Office contract W911NF-17-C-0050.
DD’s research was supported by the NSF under Grant ECCS 1710558\. DL’s
research was also sponsored by the Army Research Office and was accomplished
under Grant Number W911NF-20-1-0075. The views and conclusions contained
herein are those of the authors and should not be interpreted as necessarily
representing the official policies or endorsements, either expressed or
implied, of the ODNI, IARPA, DARPA, or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Governmental
purposes notwithstanding any copyright annotation thereon.

##  Appendix A Examples of reachability/unreachability of the ground state of
$C$

Here we provide two examples, one where Prop. 1 guarantees reachability of the
ground state of $C$ despite the algebra generated by $B$ and $C$ being smaller
than the full $su(d)$, another illustrating that $[B,C]\neq 0$ is not a
sufficient condition.

###  A.1 First example

Let $B=\sigma^{z}\otimes\leavevmode{\rm 1\mkern-4.8muI}$,
$C=\sigma^{x}\otimes\sigma^{x}+\sigma^{y}\otimes\sigma^{y}+\sigma^{z}\otimes\sigma^{z}$.
Both $B$ and $C$ commute with $M=\sigma^{z}\otimes\leavevmode{\rm
1\mkern-4.8muI}+\leavevmode{\rm 1\mkern-4.8muI}\otimes\sigma^{z}$. We know
that the ground state of $C$ is a singlet
$\frac{1}{\sqrt{2}}(\ket{\downarrow\uparrow}-\ket{\uparrow\downarrow})$ with
total spin zero and hence belongs to the sector $M=0$. This statement follows
from a theorem of Marshall Marshall (1955); Lieb and Mattis (1962) and holds
for general antiferromagnetic Heisenberg models defined on a bipartite
lattice. In particular, it does not require knowledge of the ground state. The
projector $P_{0}$ is given by

$P_{0}=|{\downarrow\uparrow}\rangle\\!\langle\downarrow\uparrow|+|{\uparrow\downarrow}\rangle\\!\langle\uparrow\downarrow|\
.$ (47)

In $\mathrm{Ran}(P_{0})$ one has:

$\displaystyle P_{0}BP_{0}$ $\displaystyle=\sigma^{z}$ (48a) $\displaystyle
P_{0}CP_{0}$ $\displaystyle=-\leavevmode{\rm 1\mkern-4.8muI}+2\sigma^{x}\ .$
(48b)

Therefore in $\mathrm{Ran}(P_{0})$, $B$ and $C$ generate the full $su(2)$
algebra and any state in $\mathrm{Ran}(P_{0})$ can be reached starting from
any state in $\mathrm{Ran}(P_{0})$, in particular the ground state of $C$.

###  A.2 Second example

It is straightforward to find examples where the ground state of $C$ cannot be
reached even when $[B,C]\neq 0$. Consider, e.g.,
$B=\sigma^{x}\otimes\leavevmode{\rm 1\mkern-4.8muI}$ and
$C=\sigma^{z}\otimes\leavevmode{\rm 1\mkern-4.8muI}+\leavevmode{\rm
1\mkern-4.8muI}\otimes\sigma^{z}$. In this case the Lie algebra generated by
$B$ and $C$ is $su(2)\otimes u(1)\neq su(4)$ and only the first qubit can be
fully steered anywhere on the Bloch sphere. As a consequence, the ground state
$\ket{\downarrow\downarrow}$ of $C$ cannot be reached unless one starts with a
state of the form
$\rho_{0}=\tilde{\rho}\otimes|{\downarrow}\rangle\\!\langle\downarrow|$, with
$\tilde{\rho}$ being an arbitrary single qubit state.

##  Appendix B General results on optimal control; the Pontryagin Maximum
Principle

We review here standard results in optimal control theory emphasizing a
geometric viewpoint and the results needed for the applications in the main
body of the paper.

###  B.1 Setup

In optimal control theory (see, e.g., Ref. Fleming and Rishel (1975)) one
considers a general control system

$\dot{x}=f\left(x,s\right),\qquad x(0)=x_{0},$ (49)

with $x\in\mathbb{R}^{N}$,101010In our case $x$ is either a wavefunction
$\psi$ ($N=n$) or a density matrix $\rho$ ($N=n^{2}$), after the
coordinatization described in Appendix C. the control $s$ with values from a
compact subset $S\subseteq\mathbb{R}^{M}$, $f$ a smooth map
$\mathbb{R}^{N}\times\mathbb{R}^{M}\rightarrow\mathbb{R}^{N}$ that does not
depend explicitly on time. The terminal cost (of _Mayer type_ 111111In control
theory, problems with a cost that depends only on the final values $t_{f}$ and
$x(t_{f})$, are called of Meyer type. In our case the cost does not depend
explicitly on time so we simply consider cost of the form $\phi(x(t_{f}))$ to
simplify the notation.)

$J:=\phi\left(x(t_{f})\right)\ ,$ (50)

is to be minimized at the terminal (final) time $t_{f}$, where $\phi$ is a
smooth function. In particular, and this is the case that interests us, one
can fix $t_{f}$ so that $J$ only depends on the final state $x(t_{f})$.

The geometric approach to the necessary condition of optimal control (see,
e.g., Ref. Agrachev and Sachkov (2004)) is based on the concept of a reachable
set (or attainable set) for Eq. (49) with values of the control in $S$, which
we denote by $\mathfrak{R}_{t}$. The set $\mathfrak{R}_{t}$ is the set of
values for the state $x$ that can be reached (from $x_{0}$) at time _exactly_
$t$ for control functions with values in the set $S$. With this definition,
the minimum of the cost $J$ in Eq. (50) is the minimum of the function $\phi$
over $\mathfrak{R}_{t_{f}}$. One also defines the reachable set
$\mathfrak{R}_{\leq t_{f}}:=\cup_{0\leq t\leq t_{f}}\mathfrak{R}_{t}$, and if
$\mathfrak{R}_{t}$ is non decreasing with $t$,
$\mathfrak{R}_{t_{f}}=\mathfrak{R}_{\leq t_{f}}$. This is the case for Eq.
(15) if one assumes, as we do, that $\rho_{0}$ commutes with $B$ (in this case
the system can remain in the state $\rho_{0}$ for an arbitrary length of time
with the choice $s\equiv 1$). If the set of admissible controls $\tilde{S}$ is
compact (as assumed in the main text) and under general conditions on the map
$f$ in Eq. (49), which are also satisfied in our cases, Filippov’s theorem
(see, e.g., (Agrachev and Sachkov, 2004, Th. 10.1)) states that the reachable
sets are _compact_ and this implies the existence of the minimum of the
function $\phi$ and therefore of the optimal control.121212Filippov’s
existence theorem can also be applied in the version of Theorem 2.1 in
(Fleming and Rishel, 1975, Ch. III) if we assume [as for Eq. (15)] that $f$ in
Eq. (49) is linear in both the control $s$ and state $x$, the cost $J$ in Eq.
(7) is smooth, and the set $\tilde{S}$ is compact. Application of this theorem
also takes into account that $\rho$ in Eq. (7) is bounded. The only condition
that is not directly satisfied for Theorem 2.1 in (Fleming and Rishel, 1975,
Ch. III) is condition (c) for which we can, however, use the alternative (c’)
of Corollary 2.2 in (Fleming and Rishel, 1975, Ch. III), i.e., the existence
of a constant $\mu_{1}$ such that
$\\{x_{1}\,|\,J=\phi(x_{1})\leq\mu_{1}\\}\subseteq\mathbb{R}^{n}$ is compact.
The introduction of the concept of a reachable set effectively reduces the
optimal control problem to a _static_ optimization problem for the function
$\phi$, where the set of possible dynamics is described by the reachable set
$\mathfrak{R}_{t}$, that roughly separates the minimization problem from the
analysis of the dynamics.

###  B.2 The Pontryagin Maximum Principle (PMP)

The basic necessary conditions of optimality are given by the Pontryagin
Maximum Principle (PMP), which we restate below in a more general formulation
than in the main text, but in a context relevant for the problem of interest
to us, i.e., a fixed final time and a free final state. Assume that
$s^{*}=s^{*}(t)$ is the optimal control function and $x^{*}=x^{*}(t)$ the
optimal trajectory. We shall refer to $(x^{*},s^{*})$ as an _optimal pair_. We
have the following.

######  Theorem 8.

Assume that $(x^{*},s^{*})$ is an optimal pair. Then there exists a nonzero
vector of functions $p=p(t)\in\mathbb{R}^{n}$ called the co-state that
satisfies the terminal problem

$\dot{p}^{T}=-p^{T}\frac{\partial f}{\partial
x}\left(x^{*}(t),s^{*}(t)\right),\quad
p^{T}\left(t_{f}\right)=-\frac{\partial\phi}{\partial
x}\left(x^{*}(t_{f})\right),\ $ (51)

with $f$ defined in Eq. (49) and $\phi$ defined in Eq. (50). Furthermore,
define the Hamiltonian function

$\mathbb{H}\left(p,x,s\right):=p^{T}f\left(x,s\right)\ .$ (52)

Then we have (maximum principle):

$\mathbb{H}\left(p,x^{*},s^{*}\right)=\max_{v\in\tilde{S}}\mathbb{H}\left(p,x^{*},v\right)\
,$ (53)

(where $\tilde{S}$ is the set of admissible controls) and

$\mathbb{H}\left(p(t),x^{*}(t),s^{*}(t)\right)\equiv\lambda,$ (54)

for a constant $\lambda$.

The constant $\lambda$ describes the dependence of the optimal cost on the
terminal time $t_{f}$. To see this, given the optimal control $s^{*}$ defined
in $[0,t_{f}]$, calculate the variation of the cost with this control at
$t=t_{f}$,

$\displaystyle\left.\frac{d}{dt}\right|_{t=t_{f}}\phi(x(t))$
$\displaystyle=\frac{\partial\phi}{\partial
x}(x(t_{f}))f\left(x(t_{f}),s(t_{f})\right)$
$\displaystyle=-p^{T}\left(t_{f}\right)f\left(x(t_{f}),s(t_{f})\right)$
$\displaystyle=-\lambda\ ,$ (55)

where we used Eqs. (49), (51) and (54). In particular, $\lambda>0$ indicates
that it is possible to lower the cost by increasing the time or, in other
words, the constraint $t=t_{f}$ is _active_. If $\lambda=0$ the constraint on
the final time is not active. If $\lambda<0$, the above calculation shows that
the cost is actually _increasing_ with $t$ at $t=t_{f}$.

If there exists a value $s_{0}$ in the admissible control set $\tilde{S}$ such
that $f(x_{0},s_{0})=0$ in Eq. (49), we can show that we must have
$\lambda\geq 0$. The argument is as follows. Assume $\lambda<0$ in Eq. (55).
Since $\phi(x(t))$ is increasing in $t$, there exists an $\epsilon>0$ such
that $\phi(x(t_{f}-\epsilon))<\phi(x(t_{f}))=J$. Now construct the following
control function

$s_{1}(t)=\begin{cases}s_{0}&t\in[0,\epsilon]\\\
s^{*}(t-\epsilon)&t\in(\epsilon,t_{f}]\ ,\end{cases},$ (56)

i.e., $s_{1}$ leaves the cost unchanged in the first interval of time
$[0,\epsilon]$ and then follows the optimal schedule shifted by $\epsilon$.
Let us denote by $\phi_{1}(t)$ the cost function at time $t$ obtained with
control $s_{1}$. Then, by construction,
$\phi_{1}(t_{f})=\phi(x(t_{f}-\epsilon))<\phi(x(t_{f}))$, which contradicts
the fact that $s^{*}$ is optimal. Therefore $\lambda$ cannot be negative.
Summarizing, we have:

######  Proposition 1.

Assume there exists a value $s_{0}\in\tilde{S}$ such that $f(x_{0},s_{0})=0$
in Eq. (49). Then $\lambda\geq 0$ in Eq. (54), and $\lambda>0$ if and only if
the constraint on the final time $t_{f}$ is active.

For the problem of interest in the main body of the paper the above assumption
on the existence of the value $s_{0}\in\tilde{S}$ is valid in Theorems 2 and
6, since we assume that the initial condition commutes with the Hamiltonian
$B$ (i.e., we can choose $s_{0}=1$). Furthermore, we note, concerning
reachable sets, that (i) because of the existence of such a value
$s_{0}\in\tilde{S}$, we have ${\cal R}_{\leq t}={\cal R}_{t}$; (ii) from
standard result in control theory, (e.g., Ref. (Ayala  _et al._ , 2017, Th. 1)
and references therein) we know that for the bilinear class of models [such as
Eq. (15)], the sets ${\cal R}_{\leq t}={\cal R}_{t}$ are compact and
continuous with $t$ with respect to the Hausdorff metric.

Since the optimal cost $J$ is the minimum of the continuous function
$\phi=\phi(x)$ on $\mathfrak{R}_{t_{f}}$, it depends continuously on $t_{f}$,
and since $\mathfrak{R}_{t_{f}}$ is nondecreasing with $t_{f}$, the optimal
cost is _nonincreasing_ with $t_{f}$. That the constraint on the final time is
active means that the cost is actually strictly _decreasing_. Notice in
particular that the function $\phi=\operatorname{Tr}(C\rho)$ giving the cost
in Eq. (7) is _linear_ in the state ($\rho$); hence the minimum is necessarily
achieved on the _boundary_ of the reachable set.131313For a linear function
$\phi(x)=k^{T}x$, if the minimum is achieved at an interior point setting the
derivatives equal to zero would imply $k=0$. Therefore $\lambda>0$ in Eq. (54)
implies that at the optimal final point, the boundary of the reachable set
$\mathfrak{R}_{t_{f}}$ “moves” in such a way so as to make the cost decrease.

##  Appendix C Coordinatization in terms of an orthonormal real matrix basis

The application to quantum systems of the PMP, which is typically formulated
over real vector spaces as in Appendix B, has to account for the fact that in
the quantum case the equations are complex-valued. To show how this can be
done we start with some basic preliminaries.

We denote the Hilbert-Schmidt scalar product between operators $A$ and $B$
acting on an $n$-dimensional Hilbert space $\mathcal{H}$ by $\langle
A,B\rangle:=\operatorname{Tr}(A^{\dagger}B)$. For superoperators
$\mathcal{L}$, we denote the Hilbert-Schmidt adjoint of $\mathcal{L}$ by
$\mathcal{L}^{\dagger}$, which is defined via

$\langle\mathcal{L}^{\dagger}(A),B\rangle:=\langle A,\mathcal{L}(B)\rangle\ \
\forall A,B\ .$ (57)

We now choose an orthonormal basis $\\{F_{j}\\}$ for the _real_ vector space
of Hermitian $n\times n$ matrices with Hilbert-Schmidt inner product $\langle
A,B\rangle=\operatorname{Tr}(AB)$. Let us “coordinatize”
$X=\sum_{j}\boldsymbol{X}_{j}F_{j}$ in this basis, where henceforth we use the
notation $\boldsymbol{X}=\\{\boldsymbol{X}_{j}\\}$ for the vector of real-
valued coordinates of the operator $X$. Then:

$\displaystyle\langle A,B\rangle$
$\displaystyle=\sum_{jk}\boldsymbol{A}_{j}\boldsymbol{B}_{k}\operatorname{Tr}(F_{j}F_{k})=\sum_{jk}\boldsymbol{A}_{j}\boldsymbol{B}_{k}\delta_{jk}$
$\displaystyle=\boldsymbol{A}^{T}\boldsymbol{B}\ ,$ (58)

so that in these coordinates the inner product $\langle A,B\rangle$
corresponds to the standard inner product in $\mathbb{R}^{n^{2}}$. In
particular, the Hermitian density operator
$\rho=\sum_{i}\boldsymbol{\rho}_{i}F_{i}$, is now represented by a real,
$n^{2}$-dimensional vector $\boldsymbol{\rho}$.141414This representation is
closely related to the familiar “coherence vector” $\boldsymbol{v}$, wherein
one represents $\rho$ as $\rho=\frac{1}{n}(\leavevmode{\rm
1\mkern-4.8muI}+\sum_{j=1}^{n^{2}-1}\boldsymbol{v}_{j}F_{j})$, where
$F_{n^{2}}:=\leavevmode{\rm 1\mkern-4.8muI}$ and the other orthonormal basis
elements $\\{F_{j}\\}_{j=1}^{n^{2}-1}$ are in addition traceless; see, e.g.,
Ref. M.S. Byrd and N. Khaneja (2003). Since the Liouvillian is Hermitian
preserving, i.e.,
$\left[\mathcal{L}\left(X\right)\right]^{\dagger}=\mathcal{L}\left(X^{\dagger}\right)\forall
X$, after coordinatization the operator $\mathcal{L}$ can be seen as an
operator $\mathbb{R}^{n^{2}}\mapsto\mathbb{R}^{n^{2}}$. Indeed, denoting the
corresponding matrix by $\boldsymbol{\mathcal{L}}$ in the chosen basis, i.e.,
$\boldsymbol{\mathcal{L}}_{ij}=\operatorname{Tr}\left(F_{i}\mathcal{L}\left(F_{j}\right)\right)$,
one has

$\displaystyle\boldsymbol{\mathcal{L}}_{ij}^{*}$
$\displaystyle=\operatorname{Tr}\left[(F_{i}\mathcal{L}\left(F_{j}\right))^{\dagger}\right]=\operatorname{Tr}\left(\left[\mathcal{L}\left(F_{j}\right)\right]^{\dagger}F_{i}^{\dagger}\right)$
$\displaystyle=\operatorname{Tr}\left(\mathcal{L}\left(F_{j}\right)F_{i}\right)=\boldsymbol{\mathcal{L}}_{ij}\
.$ (59)

i.e., the matrix $\boldsymbol{\mathcal{L}}$ is real and defines an operator
$\mathbb{R}^{n^{2}}\mapsto\mathbb{R}^{n^{2}}$. Accordingly, Eq. (6) is
transformed into a real-valued equation:

$\dot{\boldsymbol{\rho}}=\boldsymbol{\mathcal{L}}\boldsymbol{\rho},\qquad\boldsymbol{\rho}(0)=\boldsymbol{\rho}_{0}\
.$ (60)

The cost (7) takes the form

$J=\boldsymbol{C}^{T}\boldsymbol{\rho}(t_{f})\ .$ (61)

We are now ready to state the PMP in the standard setting of real-valued
functions, in the form needed for our purposes. Namely:

######  Theorem 9.

Assume $(\boldsymbol{\rho}^{*},s^{*})$ is an optimal pair for the problem
defined by Eqs. (60) and (61) for a fixed final time $t_{f}$.151515Recall that
here the asterisk denotes optimality, not complex conjugation. Then there
exists a co-state vector $\boldsymbol{p}$ that satisfies161616In the notation
of Appendix B the co-state vector satisfies the adjoint equation (51), and in
the present case $\frac{\partial f}{\partial x}=\boldsymbol{\mathcal{L}}$.

$\dot{\boldsymbol{p}}^{T}=-\boldsymbol{p}^{T}\boldsymbol{\mathcal{L}},\quad\boldsymbol{p}(t_{f})=-\boldsymbol{C}\
.$ (62)

Furthermore, define the PMP control Hamiltonian function

$\mathbb{H}\left(\boldsymbol{p},\boldsymbol{\rho},s\right)=\boldsymbol{p}^{T}\mathcal{L}\boldsymbol{\rho}\
.$ (63)

We then have the _maximum principle_ :

$\mathbb{H}\left(\boldsymbol{p}(t),\boldsymbol{\rho}^{*}(t),s^{*}(t)\right)=\max_{v\in[0,1]}\mathbb{H}\left(\boldsymbol{p}(t),\boldsymbol{\rho}^{*}(t),v\right)\
,$ (64)

and there exists a real non-negative constant $\lambda$ such that

$\mathbb{H}\left(\boldsymbol{p}(t),\boldsymbol{\rho}^{*}(t),s^{*}(t)\right)=\lambda\
.$ (65)

Theorem 1 is readily obtained applying coordinatization in reverse. To see
that Eq. (62) corresponds to Eq. (8), let us write Eq. (62) explicitly as
$\dot{\boldsymbol{p}_{i}}=-\sum_{j}\boldsymbol{\mathcal{L}}_{ji}\boldsymbol{p}_{j}$.
Next, since $\boldsymbol{\mathcal{L}}$ is real we have
$\boldsymbol{\mathcal{L}}_{ji}=\langle F_{j},\mathcal{L}F_{i}\rangle=\langle
F_{j},\mathcal{L}F_{i}\rangle^{*}=\langle\mathcal{L}F_{i},F_{j}\rangle=\langle
F_{i},\mathcal{L}^{{\dagger}}F_{j}\rangle$. Thus

$\displaystyle\dot{p}$
$\displaystyle=\sum_{i}\dot{\boldsymbol{p}_{i}}F_{i}=-\sum_{ij}\boldsymbol{\mathcal{L}}_{ji}\boldsymbol{p}_{j}F_{i}$
(66a) $\displaystyle=-\sum_{ij}\langle
F_{i},\mathcal{L}^{{\dagger}}F_{j}\rangle\boldsymbol{p}_{j}F_{i}$ (66b)
$\displaystyle=-\sum_{i}\langle F_{i},\mathcal{L}^{{\dagger}}p\rangle
F_{i}=-\mathcal{L}^{{\dagger}}p\ .$ (66c)

Finally, the correspondence between Eqs. (63) and (10) is a direct consequence
of Eq. (58): $\langle
p,\mathcal{L}\rho\rangle=\boldsymbol{p}^{T}\mathcal{L}\boldsymbol{\rho}$.

##  Appendix D Proof of various formulas

###  D.1 Proof of Eq. (14)

The proof is, for for arbitrary operators $A,B,X$:

$\displaystyle\langle A,\mathcal{K}_{X}^{{\dagger}}(B)\rangle$
$\displaystyle=\langle\mathcal{K}_{X}(A),B\rangle$ (67a)
$\displaystyle=\langle-i[X,A],B\rangle=i\operatorname{Tr}([X,A]^{\dagger}B)$
(67b)
$\displaystyle=i\operatorname{Tr}(A^{\dagger}X^{\dagger}B-A^{\dagger}BX^{\dagger})$
(67c) $\displaystyle=i\operatorname{Tr}(A^{\dagger}[X^{\dagger},B])=i\langle
A,[X^{\dagger},B]\rangle$ (67d) $\displaystyle=\langle
A,-\mathcal{K}_{X^{\dagger}}(B)\rangle\ ,$ (67e)

where we used $\langle X,Y\rangle=\operatorname{Tr}[X^{\dagger}Y]$.

###  D.2 Proof of Eq. (34)

Let $\mathcal{D}\rho:=[W\rho,V^{\dagger}]+[W,\rho V^{\dagger}]$. Then:

$\displaystyle\langle X,\mathcal{D}\rho\rangle$
$\displaystyle=\operatorname{Tr}\big{(}X^{\dagger}[W\rho,V^{\dagger}_{\alpha}]+X^{\dagger}[V,\rho
W^{\dagger}]\big{)}$
$\displaystyle=\operatorname{Tr}\Big{[}\big{(}V^{\dagger}X^{\dagger}W-X^{\dagger}V^{\dagger}W$
$\displaystyle\qquad+W^{\dagger}X^{\dagger}V-W^{\dagger}VX^{\dagger}\big{)}\rho\Big{]}$
(68a)
$\displaystyle=\operatorname{Tr}\Big{[}\big{(}W^{\dagger}XV-W^{\dagger}VX$
$\displaystyle\qquad+V^{\dagger}XW-XV^{\dagger}W\big{)}^{\dagger}\rho\Big{]}$
(68b)
$\displaystyle=\operatorname{Tr}\Big{[}\big{(}W^{\dagger}[X,V]+[V^{\dagger},X]W\big{)}^{\dagger}\rho\Big{]}$
(68c) $\displaystyle=\langle\mathcal{D}^{\dagger}X,\rho\rangle\ ,$ (68d)

which yields Eq. (34) when we replace $W$ by $W_{\alpha\beta}$, $V$ by
$S_{\alpha}=S_{\alpha}^{\dagger}$, and sum over $\alpha,\beta$.

###  D.3 Proof of $v_{3}(t)\geq 0$ for sufficiently small $t$

Let us compute the Dyson series solution of Eq. (22a) to second order:

$\displaystyle\boldsymbol{v}(t)$ $\displaystyle=\Big{[}\leavevmode{\rm
1\mkern-4.8muI}+\int_{0}^{t}dt_{1}\ \boldsymbol{\mathcal{M}}(s(t_{1}))+$ (69)
$\displaystyle\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\
\boldsymbol{\mathcal{M}}(s(t_{1}))\boldsymbol{\mathcal{M}}(s(t_{2}))+O(t^{3})\big{]}\boldsymbol{v}_{0}\
.$

Using Eq. (22e), for the initial condition $\boldsymbol{v}_{0}=(1,0,0)^{T}$
there is no contribution from the first (and in fact also the third) order,
while the second order contributes via
$\left[\boldsymbol{\mathcal{M}}(s(t_{1}))\boldsymbol{\mathcal{M}}(s(t_{2}))\right]_{31}=s(t_{1})(1-s(t_{2}))\geq
0$. Hence, using $s(t)\in[0,1]$,
$v_{3}(t)=\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\
s(t_{1})(1-s(t_{2}))+O(t^{4})\geq 0$ for sufficiently small $t$. In contrast,
for the initial condition $\boldsymbol{v}_{0}=(-1,0,0)^{T}$ we have
$v_{3}(t)\leq 0$ by the same argument.

(a)

(b)

Figure 3: (a) Behavior of the optimal cost in the regular case where the
reachable set increases with the final time at the optimal point
$\boldsymbol{\rho}(t_{f})$. In this case we expect $\lambda>0$. (b) Behavior
of the optimal cost in the case where the reachable set increases with the
final time but not at the optimal final point $\boldsymbol{\rho}(t_{f})$. In
this case we expect $\lambda=0$.

##  Appendix E Optimal control and the geometry of the reachable set

The cost considered in this work is linear in the state $\rho$, which after
coordinatization we have identified with a point $\boldsymbol{\rho}$ in
$\mathbb{R}^{n^{2}}$, that is, $J:={\bf C}^{T}\boldsymbol{\rho}$. In
$\mathbb{R}^{n^{2}}$, we also consider the reachable sets $\mathfrak{R}_{t}$
at various times $t$. It is of interest to consider the level lines
(hyperplanes) ${\bf C}^{T}\boldsymbol{\rho}=k$ for various $k$’s. If ${k}$ is
the minimum cost at the final time $t_{f}$, the level line ${\bf
C}^{T}\boldsymbol{\rho}=k$ intersects the boundary of the reachable set
$\mathfrak{R}_{t_{f}}$ at the point $\boldsymbol{\rho}(t_{f})$. In our case,
the reachable sets $\mathfrak{R}_{t}$ are always nondecreasing with $t$.
Figure 3(a) describes the regular situation of an active time constraint
($\lambda>0$ in the main text). The intersection occurs at a point where the
reachable set is increasing with time. Therefore an increase (decrease) of the
final time $t_{f}$ results in a decrease (increase) of the optimal cost.
However, in principle, a different situation may occur which is described in
Fig. 3(b). In this case, the reachable set $\mathfrak{R}_{t}$ increases with
$t$ but not at the point where the optimum occurs. In this case the final time
constraint is not active ($\lambda=0$). Notice that by continuity of the
reachable set $\mathfrak{R}_{t}$ with $t$, the point
$\boldsymbol{\rho}(t_{f})$ where the optimum was achieved with the final time
$t_{f}$ will give the optimal for $t_{f}+\epsilon$ for sufficiently small
$\epsilon>0$. The corresponding control will be a zero control which keeps the
state at the initial value for time $\epsilon$ followed by the same control
applied to reach $\boldsymbol{\rho}(t_{f})$.

Since we have not claimed uniqueness of the optimal control, the two
situations may occur simultaneously for two different optimal trajectories.
For one of them the final time constraint is active, while for the other one
it is not. Given the additional structure of our problem, it should be
possible to say more about the geometry of the reachable sets for the systems
of interest here beside what is known from, for instance, Ref. Ayala  _et al._
(2017). However, this is beyond the scope of this work.

##  Appendix F Proof of Theorem 3

###### Proof.

Let us denote by $\tilde{J}(t_{f})$ the cost function obtained with control
$\tilde{s}:=\tilde{s}(t)$. By definition, the optimal control satisfies
$J_{\min}(t_{f})\leq\tilde{J}(t_{f})$. Hence, to prove
$J_{\min}(t_{f})<J_{\min}(0)=\langle C,\rho_{0}\rangle$, it suffices to show
that there exists a control $\tilde{s}$ for Eq. (6) such that the
corresponding cost $\tilde{J}(t_{f})$ satisfies
$\tilde{J}(t_{f})<J_{\min}(0)$. This is the proof strategy we employ here.

Specifically, we consider a bang-bang control schedule $\tilde{s}$ in the
interval $[0,2t]$ and denote the corresponding cost starting from $\rho_{0}$
by $\tilde{J}(2t)$. We will show that $\tilde{J}(2t)<\tilde{J}(0)$ for
sufficiently small $t$, which gives

$J_{\min}(2t)\leq\tilde{J}(2t)<\tilde{J}(0)=J_{\min}(0)=\operatorname{Tr}(C\rho_{0})\
,$ (70)

and this proves the theorem with $t_{f}=2t$.

The class of controls we consider is $\tilde{s}\equiv 1$ [corresponding to
$\mathcal{K}_{C}$ in Eq. (15)] for an interval of length $t$, followed by
$\tilde{s}\equiv 0$ [corresponding to $\mathcal{K}_{B}$ in Eq. (15)] for a
second interval of length $t$. This gives for ${J}(t)$ [Eq. (7)]:

$\tilde{J}(2t)=\operatorname{Tr}\left(Ce^{-iBt}e^{-iCt}\rho_{0}e^{iCt}e^{iBt}\right)\
.$ (71)

We work in a basis where $B$ is diagonal with eigenvalues in decreasing order:
$B=\texttt{diag}(\lambda_{n},\lambda_{n-1},\dots,\lambda_{1})$ and
$\lambda_{1}<\lambda_{j}$, for each $j=2,3,\dots,n$ (nondegeneracy). In this
basis $\rho_{0}=\texttt{diag}(0,\dots,0,1)$.

In Eq. (71), the Baker-Campbell-Hausdorff (BCH) formula yields:

$\displaystyle\tilde{J}(2t)=\operatorname{Tr}\Big{[}e^{iBt}Ce^{-iBt}\times$
$\displaystyle\qquad\Big{(}\rho_{0}-i[C,\rho_{0}]t-[C,[C,\rho_{0}]]\frac{t^{2}}{2}+O(t^{3})\Big{)}\Big{]}$
(72)

Applying the BCH formula again, this time to $e^{iBt}Ce^{-iBt}$, we obtain:

$\tilde{J}(2t)=\operatorname{Tr}\Bigg{[}\left(C+i[B,C]t-[B,[B,C]]\frac{t^{2}}{2}\right)\times\\\
\left(\rho_{0}-i[C,\rho_{0}]t-[C,[C,\rho_{0}]]\frac{t^{2}}{2}\right)\Bigg{]}+O(t^{3}).$
(73)

Expanding, we obtain:

$\displaystyle\tilde{J}(2t)$
$\displaystyle=\operatorname{Tr}(C\rho_{0})+i\operatorname{Tr}([B,C]\rho_{0})t-\operatorname{Tr}([B,[B,C]]\rho_{0})\frac{t^{2}}{2}$
$\displaystyle\quad-i\operatorname{Tr}(C[C,\rho_{0}])t+\operatorname{Tr}\left([B,C][C,\rho_{0}]\right)t^{2}$
$\displaystyle\quad-\operatorname{Tr}\left(C[C,[C,\rho_{0}]]\right)\frac{t^{2}}{2}+O(t^{3}).$
(74)

Several of the terms in the above equation vanish. In particular,

$\operatorname{Tr}\left([B,C]\rho_{0}\right)=\operatorname{Tr}\left([\rho_{0},B]C\right)=0\
,$ (75)

since $B$ and $\rho_{0}$ commute. $\operatorname{Tr}([B,[B,C]]\rho_{0})=0$ for
the same reason.
$-i\operatorname{Tr}(C[C,\rho_{0}])=-i\operatorname{Tr}(\rho_{0}[C,C])=0$, and
$\operatorname{Tr}\left(C[C,[C,\rho_{0}]]\right)=\operatorname{Tr}\left([C,\rho_{0}][C,C]\right)=0$.
Therefore we have

$\tilde{J}(2t)-\operatorname{Tr}(C\rho_{0})=\operatorname{Tr}\left([B,C][C,\rho_{0}]\right)t^{2}+O(t^{3})\
.$ (76)

Write

$B=\begin{pmatrix}\Lambda&0\\\ 0&\lambda_{1}\end{pmatrix},\qquad
C=\begin{pmatrix}C_{1}&a\\\ a^{\dagger}&c\end{pmatrix}\ ,$ (77)

with $\Lambda=\texttt{diag}(\lambda_{n},\dots,\lambda_{2})$, $C_{1}$ an
$(n-1)\times(n-1)$ Hermitian matrix, $c$ a real number and $a$ an $(n-1)$-th
dimensional complex vector. With these notations, we have

$\displaystyle[B,C]$
$\displaystyle=\begin{pmatrix}[\Lambda,C_{1}]&(\Lambda-\lambda_{1}\leavevmode{\rm
1\mkern-4.8muI})a\\\ a^{\dagger}(\lambda_{1}\leavevmode{\rm
1\mkern-4.8muI}-\Lambda)&0\end{pmatrix}$ (78) $\displaystyle{}[C,\rho_{0}]$
$\displaystyle=\begin{pmatrix}0&a\\\ -a^{\dagger}&0\end{pmatrix}.$ (79)

From this we obtain:

$\operatorname{Tr}\left([B,C][C,\rho_{0}]\right)=2a^{\dagger}(\lambda_{1}\leavevmode{\rm
1\mkern-4.8muI}-\Lambda)a\ ,$ (80)

so that we have, from Eq. (76):

$\tilde{J}(2t)-\operatorname{Tr}(C\rho_{0})=-t^{2}\left(\sum_{j=2}^{n}(\lambda_{j}-\lambda_{1})|a_{j}|^{2}\right)+O(t^{3})\
,$ (81)

where $a_{j}$ are the components of $a$. Since $\lambda_{j}>\lambda_{1}$ for
each $j$, we have for sufficiently small $t$:

$\tilde{J}(2t)-\operatorname{Tr}(C\rho_{0})=\tilde{J}(t)-J_{\min}(0)<0\ ,$
(82)

as required. We have assumed here that at least one of the components of $a$
is nonzero. If that were not the case then $[\rho_{0},C]=0$ and $\rho_{0}$
would be fixed not just under $B$ but also under $C$. There would then be no
dynamics, which is a case that is naturally excluded. ∎

##  Appendix G Singular arcs

Along singular arcs we have $x_{C}\equiv x_{B}$, i.e.,

$\langle B,S\rangle\equiv\langle C,S\rangle\ .$ (83)

Differentiating Eq. (83), using Eq. (39) we find $\langle
B,\mathcal{K}_{H}S\rangle=\langle C,\mathcal{K}_{H}S\rangle$. Using
$H=sB+(1-s)C$ and the antihermiticity of $\mathcal{K}_{H}$ we thus obtain:

$\displaystyle\langle\mathcal{K}_{H}^{\dagger}(B),S\rangle$
$\displaystyle=\langle\mathcal{K}_{H}^{\dagger}(C),S\rangle\Longrightarrow$
(84a) $\displaystyle\langle(1-s)[C,B]^{\dagger},S\rangle$
$\displaystyle=\langle s[B,C]^{\dagger},S\rangle\Longrightarrow$ (84b)
$\displaystyle-(1-s)\langle[C,B],S\rangle$
$\displaystyle=s\langle[C,B],S\rangle\ \Longrightarrow\ $ (84c)
$\displaystyle\langle[C,B],S\rangle$ $\displaystyle=0\ .$ (84d)

Analogously, differentiating Eq. (84) and using the antihermiticity of
$\mathcal{K}_{H}$ again, setting $D:=[C,B]$, we have:

$\displaystyle 0$ $\displaystyle=\langle D,\dot{S}\rangle=\langle
D,\mathcal{K}_{H}(S)\rangle=\langle\mathcal{K}_{H}^{\dagger}(D),S\rangle$
(85a)
$\displaystyle=-(1-s)\langle[C,D]^{\dagger},S\rangle-s\langle[B,D]^{\dagger},S\rangle$
(85b) $\displaystyle=(1-s)\langle[C,D],S\rangle+s\langle[B,D],S\rangle\ .$
(85c)

Conditions (83)-(85) have to hold along a singular arc. In an algorithm to
calculate the _dynamical Lie algebra_ for the controllability of Eq. (15)
D’Alessandro (2007), the matrices $B$ and $C$ are the matrices of “depth” zero
in the calculation via iterated Lie brackets. The matrix $[C,B]$ is of depth
one and the matrices $[[C,B],C]$ and $[[C,B],B]$ are of depth two. Now, one
can have either (i) $\langle[C,D],S\rangle\neq\langle[B,D],S\rangle$ or (ii)
$\langle[C,D],S\rangle=\langle[B,D],S\rangle$. In case (i) it follows from Eq.
(85) that

$s=\frac{\langle[C,D],S\rangle}{\langle[C,D],S\rangle-\langle[B,D],S\rangle}$
(86)

(compare with Ref. (Brady  _et al._ , 2021, Eq. (12))). This shows the
continuity of $s$ in the corresponding open set(s). Case (ii) implies
$\langle[C,D],S\rangle\equiv\langle[B,D],S\rangle\equiv 0$ (in some closed
set). One can further differentiate one of these equations and obtain an
analog of Eq. (85) at a higher order, at which point similar reasoning can be
applied. In principle $s$ can be defined in different intervals by equations
such as Eq. (86) or its higher order generalizations. In each interval $s$ is
continuous because of the continuity of $S$. We leave a more general proof of
continuity of $s$ on the _entire_ singular arc as an open problem. In any
case, equations (84)-(85) provide information on the dynamics along singular
arc intervals. They are used in the example discussed in Appendix H.

##  Appendix H Optimal control protocol for the spin-$\frac{1}{2}$ model

Here we analyze in detail the optimal control problem for the
spin-$\frac{1}{2}$ model treated in Sec. IV.3. Our goal is to give a simple
but explicit example to show how the results developed in this paper can be
used to find the optimal control. We shall use the same notation as in the
example of Sec. IV.3. To avoid the situation of an inactive terminal time
constraint described in the example, we assume that the initial state is the
ground state $x_{0}=(-1,0,0)^{T}$, so that we can apply Theorem 3.

###  H.1 The global minimum is found using two non-singular arcs in time
$t_{f}\geq\pi$

Recall that $C=\sigma^{z}/2$ and $B=\sigma^{x}/2$. The global minimum of the
cost (7) is $J_{\min}=\operatorname{Tr}[C\rho_{f}]=-1/2$, achieved when
$\rho_{f}=(\leavevmode{\rm 1\mkern-4.8muI}-\sigma^{z})/2$, the ground state of
$C=\sigma^{z}/2$. Given that our initial condition is
$\rho_{0}=(\leavevmode{\rm 1\mkern-4.8muI}-\sigma^{x})/2$ (the ground state of
$B$, corresponding to $\boldsymbol{v}=(-1,0,0)^{T}$), we can trivially reach
$\rho_{f}=(\leavevmode{\rm 1\mkern-4.8muI}-\sigma^{z})/2$ by applying two
consecutive bangs (i.e., unitary single-qubit gates): first $e^{-i(\pi/2)C}$
(rotation to $(\leavevmode{\rm 1\mkern-4.8muI}-\sigma^{y})/2$) with $s\equiv
0$, then $e^{-i(\pi/2)B}$ (rotation to $(\leavevmode{\rm
1\mkern-4.8muI}-\sigma^{z})/2$) with $s\equiv 1$. Each bang lasts for a time
$\frac{\pi}{2}$, therefore the total bang-bang sequence last for a total time
of $\pi$. _This sequence presents no singular arcs_. For any $t_{f}>\pi$ the
constraint on $t_{f}$ becomes inactive, i.e., increasing $t_{f}$ cannot
further lower the value of $J_{\min}$. Since we assume that the global ground
state is not reached (recall the discussion in Sec. IV.2), henceforth we
assume that $t_{f}<\pi$. In principle, this setting could still allow for the
appearance of singular arcs. However, we shall show that this is not the case.

###  H.2 Conditions on the singular arcs for the spin-$1/2$ model

Let us derive the conditions on the singular arcs in the present problem,
which are a special case of the computations carried out in Appendix G. Using
Eq. (20), we obtain $\left[C,B\right]=i\sigma^{y}/2$,
$\left[[C,B],C\right]=-\sigma^{x}/2$, $\left[[C,B],B\right]=\sigma^{z}/2$.
Using these, Eqs. (40b), (41), and (83) for the switching operator $S$ become:

$2x_{C}=\operatorname{Tr}\left(S\sigma^{z}\right)\equiv\operatorname{Tr}\left(S\sigma^{x}\right)=2x_{B}\
.$ (87)

Condition (84) becomes:

$\operatorname{Tr}\left(S\sigma^{y}\right)\equiv 0\ ,$ (88)

and condition (85) becomes
$(1-s)\operatorname{Tr}(\sigma^{x}S)-s\operatorname{Tr}(\sigma^{z}S)\equiv 0$,
which using Eq. (87) gives:

$\left(1-2s\right)\operatorname{Tr}\left(\sigma^{x}S\right)\equiv 0\ .$ (89)

Thus, either $s=1/2$ or $\operatorname{Tr}(\sigma^{x}S)\equiv 0$. Let assume
the latter. From Eqs. (87) and (88) we obtain
$\operatorname{Tr}(\sigma^{z}S)=\operatorname{Tr}(\sigma^{y}S)\equiv 0$. Since
we can expand
$S=\frac{1}{2}\sum_{i=1}^{3}\operatorname{Tr}(S\sigma_{i})\sigma_{i}$ ($S$ is
traceless since it is defined as a commutator), this would then imply that
$S\equiv 0$ on a singular interval. However, since $S$ satisfies the linear
equation (39), this would imply $S\equiv 0$ on the whole interval $[0,t_{f}]$
and, in particular, $[p_{0},\rho_{0}]=[p_{f},\rho_{f}]=-[C,\rho_{f}]=0$. This
would imply that $\rho_{f}$ is a linear combination of eigenprojectors of $C$,
but since $\rho_{f}$ is a pure state it must in fact be equal to a single
eigenprojector. Moreover, this must be the ground state of $C$ since $J(t)$ is
minimized at $t=t_{f}$. But, since $J(t_{f})\leq J(0)$, the only possibility
is that $J=\operatorname{Tr}(C\rho)$ reaches its global minimum at $t_{f}$,
which contradicts our assumption that $t_{f}$ is smaller than a value that
would allow the global minimum to be reached. Hence we conclude that
$\operatorname{Tr}(\sigma^{x}S)\neq 0$ in Eq. (89), _which yields $s\equiv
1/2$ on the singular arcs_. However, we shall see in Proposition 3 that
singular arcs are in fact not possible in this case.

###  H.3 Candidate optimal controls with a singular arc

Using conditions (87) and (88) along with Eq. (42) equated to $\lambda$, we
have that in the time interval of a singular arc

$S\equiv\lambda\left(\sigma^{x}+\sigma^{z}\right)\ .$ (90)

$\lambda=0$ is impossible because according to the argument at the end of the
previous subsection $S\equiv 0$ is to be excluded. Since $\lambda\neq 0$, we
can apply all the conclusions of Theorem 2 and affirm that the optimal control
starts with an $s\equiv 0$ bang arc and ends with an $s\equiv 1$ bang arc.
Therefore, preceding or following a singular arc we must have $s\equiv 0$ or
$s\equiv 1$, respectively. Let us show that after a singular arc we cannot go
to a switching point, i.e., where Eq. (87) holds [$(\lambda,\lambda)$ in Fig.
2]. (Analogously, changing the sign of time, we can show that a singular arc
cannot be preceded by a switching point.) Assume that after the singular arc
we have $s\equiv 0$. The switching operator $S=S(t)$, with $t=0$ at the end of
this singular arc, is then the solution of Eq. (39) with the initial condition
(90), i.e., $S(t)=\lambda
e^{-it\sigma^{z}/2}\sigma^{x}e^{it\sigma^{z}/2}+\lambda\sigma^{z}$. The
minimum time needed for it to return to a switching point is $t=2\pi$. This
contradicts the fact that $t_{f}<\pi$ and therefore is impossible. Similar
reasoning shows that we cannot go back to a switching point with $s\equiv 1$.
Therefore, we have learned the following fact about optimal control in the
single-qubit case:

######  Proposition 2.

The optimal control has at most one singular arc and, if it does, the optimal
control is the sequence $s\equiv 0$, $s\equiv\frac{1}{2}$, $s\equiv 1$.

Consider now the initial switching operator $S_{0}$, which together with the
differential equation Eq. (39) determines the control sequence. Since $S$ is
traceless, we can write $S_{0}=\boldsymbol{r}_{0}\cdot\boldsymbol{\sigma}$,
but using Eq. (37) and $\rho_{0}=\frac{1}{2}({\bf\leavevmode{\rm
1\mkern-4.8muI}}-\sigma^{x})$ we see by expanding $p$ in the Pauli matrix
basis that $S_{0}$ cannot contain $\sigma^{x}$, i.e., we find that $S_{0}$ has
the form

$S_{0}=r_{0y}\sigma^{y}+r_{0z}\sigma^{z}\ .$ (91)

In the first interval $s\equiv 0$ and therefore, from Eq. (39):

$\displaystyle S(t)$
$\displaystyle=r_{0z}\sigma^{z}+r_{0y}e^{-it\sigma^{z}/2}\sigma^{y}e^{it\sigma^{z}/2}$
(92a)
$\displaystyle=r_{0z}\sigma^{z}+r_{0y}\cos(t)\sigma^{y}-r_{0y}\sin(t)\sigma^{x}\
.$ (92b)

If there is a singular arc and therefore $S$ takes the form (90), then we must
have $t=\frac{\pi}{2}$ and $r_{0y}=-r_{0z}$ or $t=\frac{3\pi}{2}$ and
$r_{0y}=r_{0z}$. The second case is to be excluded since $t_{f}<\pi$. After
the singular arc we would have $s\equiv 1$, which, using Eq. (39) again would
give:

$\displaystyle S(t)$
$\displaystyle=r_{0z}\sigma^{x}+r_{0z}e^{-it\sigma^{x}/2}\sigma^{z}e^{it\sigma^{x}/2}$
(93a)
$\displaystyle=r_{0z}\sigma^{x}+r_{0z}(\cos(t)\sigma^{z}-\sin(t)\sigma^{y})\
.$ (93b)

Since we have to reach the point $(x_{C},x_{B})=(0,\lambda)$ in Fig. 2, we
must have $\cos\left(t\right)=0$, i.e., $t=\frac{\pi}{2}$ or
$t=\frac{3\pi}{2}$, which has to be added to the time used before the last
interval. Therefore the total time is greater than or equal to $\pi$, which we
have excluded.

In conclusion we have:

######  Proposition 3.

No singular arc exists in the optimal control for the spin-$1/2$ example with
$t_{f}<\pi$.

###  H.4 Candidate optimal controls without singular arcs

Now we consider the optimal control candidates knowing that they must be free
of singular arcs, i.e., they can consist only of bangs. Since $S(t)$ is
traceless we again use the parametrization
$S(t)=\boldsymbol{r}(t)\cdot\boldsymbol{\sigma}$. We already know [Eq. (91)]
that:

$\boldsymbol{r}(0)=\boldsymbol{r}_{0}=(0,r_{0y},r_{0z})^{T}\ .$ (94)

We know from Eq. (39) that the vector $\boldsymbol{r}$ evolves according to
Eqs. (22a)-(22e). More explicitly, let
$\boldsymbol{X}:=\boldsymbol{\mathcal{K}}_{B}$ and
$\boldsymbol{Z}:=\boldsymbol{\mathcal{K}}_{C}$. When $s\equiv 0$, $H=C$ and
$\boldsymbol{r}$ evolves according to $e^{t\boldsymbol{Z}}$, and likewise when
$s\equiv 1$ it evolves according to $e^{t\boldsymbol{X}}$, where, using Eq.
(21):

$\displaystyle e^{t\boldsymbol{X}}$
$\displaystyle=\left(\begin{array}[]{ccc}1&0&0\\\ 0&\cos(t)&-\sin(t)\\\
0&\sin(t)&\cos(t)\end{array}\right)$ (95d) $\displaystyle e^{t\boldsymbol{Z}}$
$\displaystyle=\left(\begin{array}[]{ccc}\cos(t)&-\sin(t)&0\\\
\sin(t)&\cos(t)&0\\\ 0&0&1\end{array}\right)\ .$ (95h)

Since we have shown that $\lambda>0$, from Theorem 2, the control law will
start with an $s\equiv 0$ bang arc and end with an $s\equiv 1$ bang arc. The
control law is determined by a sequence of intervals of lengths
$\\{\tau_{1},\tau_{2},\dots\\}$ where for $k$ odd (even) $\tau_{k}$ marks the
switch from $s\equiv 0$ to $s\equiv 1$ ($s\equiv 1$ to $s\equiv 0$), that is
$\boldsymbol{Z}\rightarrow\boldsymbol{X}$
($\boldsymbol{X}\rightarrow\boldsymbol{Z}$). That is,

$\boldsymbol{r}(t_{k})=\left\\{\begin{array}[]{c}e^{\tau_{k}\boldsymbol{Z}}\boldsymbol{r}(t_{k-1})\quad
k\text{ odd}\\\ e^{\tau_{k}\boldsymbol{X}}\boldsymbol{r}(t_{k-1})\quad k\text{
even}\end{array}\right.\ ,$ (96)

where $t_{k}=\sum_{i=1}^{k}\tau_{i}$ is the total time after $k$ intervals.
Note that in principle there is no guarantee that such a switching sequence is
finite, even if the total control interval is finite; this is known in the
control theory literature as the _Fuller phenomenon_ (see, e.g., Ref. Borisov
(2000)). We shall see in Remark 1 that this does not happen in our case and we
have a _finite_ sequence of intervals of lengths
$\\{\tau_{1},\tau_{2},\dots,\tau_{N}\\}$ with $N$ even (according to Theorem
2). Given our definitions, $t_{f}=\sum_{i=1}^{N}\tau_{i}$, where $\tau_{1}$
and $\tau_{N}$ are the lengths of the initial ($\boldsymbol{Z}$) and final
($\boldsymbol{X}$) arcs, respectively.

###  H.5 Characterization of the switching times

Note that the vector $\boldsymbol{r}=(r_{x},r_{y},r_{z})^{T}$ consists of the
components of $S$ along the Pauli basis, and that $r_{x}=x_{B}$ and
$r_{z}=x_{C}$ [Eq. (87)]. Recall also that, as argued in Sec. VI.1 (see Fig.
2), $(x_{C},x_{B})=(\lambda,\lambda)$ at every switching point. Hence
$r_{x}=r_{z}=\lambda$ at every switching point between nonsingular arcs in our
discussion below.

The optimal candidate control law is characterized by a sequence of intervals
of length $\tau_{1}$, $\tau_{2}$, etc. Define the sequence $\\{\Delta_{k}\\}$
recursively from the sequence $\\{\tau_{k}\\}$ via $\Delta_{0}=0$,
$\Delta_{k}=\tau_{k}-\Delta_{k-1}$, for $k=1,2,\dots$. Then:

######  Lemma 1.

For $n=1,2,\dots$, except for the $n$ corresponding to the last control
interval

$\left(-1\right)^{n}\sin\left(\Delta_{n}\right)r_{0y}=r_{0z},$ (97)

[cf. Eq. (91)] and

$\boldsymbol{r}\left(t_{n}\right)=\left(r_{0z},\cos(\Delta_{n})r_{0y},r_{0z}\right)^{T}\
.$ (98)

###### Proof.

The proof is by induction on $n$. For $n=1$ we have
$\Delta_{1}=t_{1}=\tau_{1}$. At the end of the first arc we must reach the
point $(x_{C},x_{B})=(\lambda,\lambda)$, which means that
$\boldsymbol{r}(\tau_{1})=(r_{0z},*,r_{0z})^{T}$. Thus, using Eq. (96) and
$e^{\tau_{1}\boldsymbol{Z}}$ in Eq. (95h), we obtain Eqs. (97) and (98).

Now assume Eq. (97) and Eq. (98) hold for $n-1$. If $n$ is even we have

$\displaystyle\boldsymbol{r}\left(t_{n}\right)=e^{\tau_{n}\boldsymbol{X}}\boldsymbol{r}\left(t_{n-1}\right)\
.$ (99)

Next, use Eq. (95d), and impose $\boldsymbol{r}(t_{n})=(r_{0z},*,r_{0z})^{T}$,
since $(x_{C},x_{B})=(\lambda,\lambda)$ at every switching point. Equality of
the $z$ component then gives
$\sin(\tau_{n})\cos(\Delta_{n-1})r_{0y}+\cos(\tau_{n})r_{0z}=r_{0z}$, and
using Eq. (97) with $n$ replaced by $n-1$, we obtain
$r_{0z}=\sin(\tau_{n})\cos(\Delta_{n-1})r_{0y}-\cos(\tau_{n})\sin(\Delta_{n-1})r_{0y}=\sin(\tau_{n}-\Delta_{n-1})r_{0y}=\sin(\Delta_{n})r_{0y}$.
Calculating the $y$ component of $\boldsymbol{r}(t_{n})$, we obtain
$\cos(\tau_{n})\cos(\Delta_{n-1})r_{0y}-\sin(\tau_{n})r_{0z}=\cos(\tau_{n})\cos(\Delta_{n-1})r_{0y}+\sin(\tau_{n})\sin(\Delta_{n-1})r_{0y}=\cos(\Delta_{n})r_{0y}$,
using again the inductive assumption Eq. (97). A similar calculation with
$\boldsymbol{Z}$ replacing $\boldsymbol{X}$ in Eq. (99) gives the result when
$n$ is odd. ∎

###  H.6 Determination of the Optimal Control

We now use the formulas in the above Lemma to determine the optimal control.
Define $\mu=\arcsin\left(\frac{r_{0z}}{r_{0y}}\right)$ and notice that from
Eq. (97) for $n=1$ and $\Delta_{1}=\tau_{1}$ we have
$\mu=\arcsin\left(-\sin\left(\tau_{1}\right)\right)$. Since $0<\tau_{1}<\pi$,
we have $\mu=-\tau_{1}$ for $\tau_{1}\in(0,\pi/2]$ and $\mu=-\pi+\tau_{1}$ for
$\tau_{1}\in[\pi/2,\pi)$, and in particular $\mu<0$.

Let us consider first the possibility that $0<\tau_{1}\leq\frac{\pi}{2}$. We
also have $\tau_{1}=\Delta_{1}=-\mu$. If there is more than one switch (i.e.,
$\tau_{2}>0$), then we can derive $\Delta_{2}$ from Eq. (97). We have either
$\Delta_{2}=\pi-\mu+2l\pi$ or $\Delta_{2}=\mu+2l\pi$ for integer $l$.
Recalling that $\tau_{2}=\tau_{1}+\Delta_{2}$, in the first case we have
$\tau_{2}=\pi-2\mu+2l\pi$, and in the second case $\tau_{2}=2l\pi$. The second
case is not possible because $\tau_{2}$ must be in $(0,\pi)$. The first case
is not possible either because $l\geq 0$ would contradict that the total time
must be less than $\pi$ while $l<0$ would give a negative or zero interval
$\tau_{2}$. Therefore, in the case $0<\tau_{1}\leq\frac{\pi}{2}$ there exists
only one switch and the control is simply the sequence of two bangs, one
corresponding to $s\equiv 0$ followed by one corresponding to $s\equiv 1$.
Before determining where the switch must occur, let us consider the case
$\frac{\pi}{2}<\tau_{1}<\pi$.

######  Lemma 2.

Assume that $\tau_{1}\in(\pi/2,\pi)$ and define
$\mu=\arcsin\left(\frac{r_{0z}}{r_{0y}}\right)<0$. Then $\Delta_{n}=\pi+\mu$
for $n$ odd and $\Delta_{n}=\mu$ for $n$ even.

###### Proof.

The claim follows by induction from Eq. (97). Applying it for $n=1$, since
$\Delta_{1}=\tau_{1}$, we have $\Delta_{1}=\pi+\mu$. Now, assume that the
claim is true for $n$ even. From Eq. (97) applied for $n$ even, we obtain
$\Delta_{n}=\mu+2l\pi$ or $\Delta_{n}=\pi-\mu+2l\pi$, for integer $l$. Since
$\tau_{n}=\Delta_{n-1}+\Delta_{n}$, using the inductive assumption, we obtain
in the two cases, $\tau_{n}=\pi+2\mu+2l\pi$ and $\tau_{n}=2\pi+2l\pi$,
respectively. The latter case is impossible because it would mean that
$\tau_{n}>\pi$. The first case is only possible with $l=0$, because $l>0$
would imply $\tau_{n}>2\pi$ while $l<0$ would give a negative time interval.
Since we saw above that $\Delta_{n}=\mu+2l\pi$ for $n$ even, this gives
$\Delta_{n}=\mu$ for such $n$.

Let us now prove that $\Delta_{n}=\pi+\mu$ for $n>1$ and odd. Again using Eq.
(97) we obtain either $\Delta_{n}=-\mu+2l\pi$ or $\Delta_{n}=\pi+\mu+2l\pi$.
The first case is impossible because it would mean
$\tau_{n}=\Delta_{n-1}+\Delta_{n}=2l\pi$ (using the inductive assumption). The
second case would give $\tau_{n}=\Delta_{n-1}+\Delta_{n}=\pi+2\mu+2l\pi$ which
is only possible for $l=0$. This gives $\Delta_{n}=\pi+\mu$. ∎

######  Remark 1.

One of the consequences of the above lemma is that the switching sequence is
finite, i.e., we do not have intervals between two switches which become
arbitrarily small and hence the Fuller phenomenon Borisov (2000) is ruled out
in our case. In particular if $0<\tau_{1}<\frac{\pi}{2}$ there is only one
switch, as we have seen, while if $\frac{\pi}{2}<\tau_{1}<\pi$ then we have
multiple switches with $\tau_{k}=\Delta_{k}+\Delta_{k-1}=\pi+2\mu$, which is a
constant independent of $k$.

In order to learn more about the optimal control, and rule out the second case
of $\frac{\pi}{2}<\tau_{1}<\pi$, we examine the final arc, which is of the
form $e^{\tau_{N}\boldsymbol{X}}$ ($s\equiv 1$). Recall that with the final
arc we have to reach the point $(x_{C},x_{B})=(0,\lambda)$, which imposes that
the final switching operator is of the form
$\boldsymbol{r}(t_{N})=(r_{0z},*,0)^{T}$ [see Eq. (93b) and the discussion
immediately below it]. Thus, using Eq. (99) with $n=N$ and Eq. (95d), we
obtain $\sin(\tau_{N})\cos(\Delta_{N-1})r_{0y}+\cos(t_{N})r_{0z}=0$. Using Eq.
(97) and $r_{0y}\neq 0$, we obtain
$\sin(\tau_{N})\cos(\Delta_{N-1})-\sin(\Delta_{N-1})\cos(\tau_{N})=\sin(\Delta_{N})=0$,
where $\Delta_{N}=\tau_{N}-\Delta_{N-1}$. Therefore
$\tau_{N}=\Delta_{N-1}+l\pi$ for $l$ integer. Now there are two cases:
Multiple switches or only one switch. In the case of multiple switches, we are
in the situation described in Lemma 2. We have
$\Delta_{N-1}=\Delta_{1}=\pi+\mu$. Therefore $\tau_{N}=\pi+\mu+l\pi$. The
integer $l$ must be zero because if it is positive we have $\tau_{N}>\pi$ and
if it is negative, we have a negative interval $\tau_{N}$. Therefore
$\tau_{N}=\pi+\mu$. However
$t_{f}\geq\tau_{1}+\tau_{N}=2\Delta_{1}=2(\pi+\mu)>\pi$ which is impossible.
Therefore the situation $\frac{\pi}{2}<\tau_{1}<\pi$ cannot occur. The only
possibility is the situation with $0<\tau_{1}<\frac{\pi}{2}$ with one switch
only. In this case, as above we have $\tau_{2}=\Delta_{1}+l\pi=\tau_{1}+l\pi$
with $l=0$ since again $l<0$ will give a negative time interval and $l$
positive will give total time greater than $\pi$. So $\tau_{2}=\tau_{1}$ and
the optimal control is the simplest one. This completes the proof of
Proposition 4.

##  Appendix I Additional considerations regarding the shortening of the
initial and final arcs

In principle one can be more quantitative about the shortening of the arcs
discussed in Sec. VI.2 by repeating the analysis of Ref. Venuti (2015) for the
quantity $x_{B}(\Delta t)$ [Eq. (45)]. The latter work provided a detailed
analysis of the return probability
$\mathcal{F}(t):=\left|\langle\psi|e^{-itH}|\psi\rangle\right|^{2}$, where $H$
is a time-independent Hamiltonian and $|\psi\rangle$ an initial state. One of
the results of Ref. Venuti (2015) was an explicit form for the average number
of zeroes of the equation $\mathcal{F}(t)=v$. Indeed, $\mathcal{F}$ is a
particular case of the left-hand side of Eq. (45) and can be written in that
form with $B=S_{0}=|{\psi}\rangle\\!\langle\psi|$, $C=H$. Let $N_{x_{B}}(v)$
be the average number of solutions of the equation $x_{B}(t)=v$. The number of
zeroes $N$ in a large interval of length $T$ turns out to be proportional to
$T$: $N_{x_{B}}(v)=TD_{x_{B}}(v)$ where $D_{x_{B}}$ can be computed using the
methods of Ref. Venuti (2015). Then $\Delta t\simeq
T/N_{x_{B}}(\lambda)=1/D_{x_{B}}\left(\lambda\right)$. Using (Venuti, 2015,
Eq. (6)), which applies for the special case mentioned above, we then obtain:

$\Delta t(\lambda)\simeq\frac{\sqrt{\pi}}{2}\frac{1}{\Delta
E}\sqrt{\frac{\langle x_{B}\rangle}{\lambda}}e^{\lambda/\langle x_{B}\rangle}\
,$ (100)

where $\langle
x_{B}\rangle:=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}x_{B}(t)dt$ and $\Delta
E$ is the standard deviation of the energies $\\{E_{k}\\}$ with respect to the
distribution
$\\{p_{k}=\left\lVert\Pi_{k}|\psi\rangle\right\rVert^{4}/\sum_{n}\left\lVert\Pi_{n}|\psi\rangle\right\rVert^{4}\\}$
(recall that $C=\sum_{k}E_{k}\Pi_{k}$ is the spectral resolution of $C$). The
derivation of Eq. (100) is subtle and is carried out in Ref. Venuti (2015).
Since it only represents a special case in our context, we do not pursue it
further here, and present Eq. (100) mainly to stimulate the interest of the
reader and establish a possible entry point towards a rigorous treatment of
the shortening of the initial and final arcs.

As a final comment, note that when $C$ is a classical Ising Hamiltonian of the
form $C=J\sum_{i,j=1}^{n}\sigma_{i}^{z}\sigma_{j}^{z}$, the energies are
integer multiples of $J$ and highly degenerate (therefore the spectrum is
commensurate: $\left\\{E_{k}\right\\}\subset\left\\{0,\pm J,\pm
2J,\ldots\right\\}$). This implies that the function $x_{B}(t)$ in Eq. (46) is
periodic with period not larger than $2\pi/J$ (as opposed to being almost
periodic) and this has implications for $D_{x_{B}}(\lambda)$, but the general
considerations we have outlined still hold.

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Kadowaki and Nishimori (1998) T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,”  (2014), arXiv:1411.4028 [quant-ph] .
  * Farhi  _et al._ (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv:quant-ph/0001106 (2000).
  * Kato (1950) T. Kato, J. Phys. Soc. Jap. 5, 435 (1950).
  * Jansen  _et al._ (2007) S. Jansen, M.-B. Ruskai, and R. Seiler, J. Math. Phys. 48, 102111 (2007).
  * Lidar  _et al._ (2009) D. A. Lidar, A. T. Rezakhani, and A. Hamma, J. Math. Phys. 50, 102106 (2009).
  * Ge  _et al._ (2016) Y. Ge, A. Molnár, and J. I. Cirac, Physical Review Letters 116, 080503 (2016).
  * Bapat and Jordan (2019) A. Bapat and S. Jordan, Quantum Information & Computation 19, 424 (2019).
  * Zhou  _et al._ (2020) L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Physical Review X 10, 021067 (2020).
  * Streif and Leib (2019) M. Streif and M. Leib, arXiv preprint arXiv:1901.01903 (2019).
  * Pagano  _et al._ (2020) G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, L. T. Brady, A. Deshpande, F. Liu, S. Jordan, A. V. Gorshkov, and C. Monroe, Proceedings of the National Academy of Sciences 117, 25396 (2020).
  * Roland and Cerf (2002) J. Roland and N. J. Cerf, Phys. Rev. A 65, 042308 (2002).
  * Rezakhani  _et al._ (2010) A. T. Rezakhani, A. K. Pimachev, and D. A. Lidar, Phys. Rev. A 82, 052305 (2010).
  * Rezakhani  _et al._ (2009) A. T. Rezakhani, W. J. Kuo, A. Hamma, D. A. Lidar, and P. Zanardi, Phys. Rev. Lett. 103, 080502 (2009).
  * Morita and Nishimori (2008) S. Morita and H. Nishimori, J. Math. Phys. 49, 125210 (2008).
  * Galindo and Kreinovich (2020) O. Galindo and V. Kreinovich, in  _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_ (2020) pp. 963–967.
  * Szegedy (2019) M. Szegedy, “What do QAOA energies reveal about graphs?”  (2019), arXiv:1912.12277 [quant-ph] .
  * Yang  _et al._ (2017) Z.-C. Yang, A. Rahmani, A. Shabani, H. Neven, and C. Chamon, Physical Review X 7, 021027 (2017).
  * Kipka and Ledyaev (2015) R. J. Kipka and Y. S. Ledyaev, Set-Valued Var. Anal 23, 133 (2015).
  * Lin  _et al._ (2019) C. Lin, Y. Wang, G. Kolesov, and U. Kalabić, Physical Review A 100, 022327 (2019).
  * Mbeng  _et al._ (2019a) G. B. Mbeng, R. Fazio, and G. E. Santoro, “Optimal quantum control with digitized quantum annealing,”  (2019a), arXiv:1911.12259 [quant-ph] .
  * Mbeng  _et al._ (2019b) G. B. Mbeng, R. Fazio, and G. Santoro, “Quantum annealing: a journey through digitalization, control, and hybrid quantum variational schemes,”  (2019b), arXiv:1906.08948 [quant-ph] .
  * Brady  _et al._ (2021) L. T. Brady, C. L. Baldwin, A. Bapat, Y. Kharkov, and A. V. Gorshkov, Phys. Rev. Lett. 126, 070505 (2021).
  * D’Alessandro (2007) D. D’Alessandro,  _Introduction to Quantum Control and Dynamics_ , 1st ed. (Chapman and Hall/CRC, Boca Raton, 2007).
  * Jurdjevic and Sussmann (1972) V. Jurdjevic and H. J. Sussmann, Journal of Differential Equations 12, 313 (1972).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Reviews of Modern Physics 90, 015002 (2018).
  * Hauke  _et al._ (2020) P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Reports on Progress in Physics (2020).
  * del Campo (2013) A. del Campo, Physical Review Letters 111, 100502 (2013).
  * Takahashi (2017) K. Takahashi, Physical Review A 95, 012309 (2017).
  * Crosson and Lidar (2021) E. J. Crosson and D. A. Lidar, Nature Reviews Physics (2021).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).