# Revisiting Quantum Optimal Control Theory: New Insights for the Canonical
Solutions

Katherine Castro  Departamento de Geometría y Topología, Facultad de Ciencias,
Universidad de Granada, 18071 Granada, Spain  Ignacio R. Solá  Departamento de
Química Física, Universidad Complutense de Madrid, 28040 Madrid, Spain  Juan
J. Omiste  Departamento de Química Física, Universidad Complutense de Madrid,
28040 Madrid, Spain<EMAIL_ADDRESS>
###### Abstract

In this study, we present a revision of the Quantum Optimal Control Theory
(QOCT) originally proposed by Rabitz et al [1], which has broad applications
in physical and chemical physics. First, we identify the QOCT equations as the
Euler-Lagrange equations of the functional associated to the control scheme.
In this framework we prove that the extremal functions found by Rabitz are not
continuous, as it was claimed in previous works. Indeed, we show that the
costate is discontinuous and vanishes after the measurement time. In contrast,
we demonstrate that the driving field is continuous. We also identify a new
set of continuous solutions to the QOCT. Overall, our work provides a
significant contribution to the QOCT theory, promoting a better understanding
of the mathematical solutions and offering potential new directions for
optimal control strategies.

##  1 Introduction

Minimization and optimization are standard tools in mathematics that we
encounter in many areas of science [2, 3, 4] ranging from health and biology
to space exploration or even social behaviour [5]. Furthermore, most of the
fundamental equations in physics may be derived from a variational principle,
which aims to minimize a quantity, for instance, the action or the energy [6,
4, 7]. Prominent examples are Lagrangian mechanics, general relativity or
geometric optics. In atomic and molecular physics, let us highlight the use of
the so-called Dirac-Frenkel variational principle to propagate quantum systems
with many degrees of freedom [8], which leads to the Multiconfigurational
Time-Dependent Hartree and related methods [9, 10, 11, 12, 13].

In the literature, we find Optimal Control as the methodology designed to
obtain a driving field such that a given observable is optimal, i. e., a
maximum or minimum, for which a functional is designed and its extremal
solutions are explored, ensuring that certain constraints are fulfilled [14,
15, 16, 17, 18]. It is said to be Quantum Optimal Control if the Schrödinger
equation drives the dynamics [1]. This methodology has been widely used to
control molecular orientation [19, 20, 21], state transitions and population
transfer [22, 23, 24, 25, 26], ultracold atoms systems [27, 28], ionization
[29] or even entanglement [30, 31]. It has also been used to drive Bose-
Einstein condensates [32] using the classical [33] and Gross-Pitaevskii
approximation [34] or to drive dissipative systems [35, 36, 37, 38].

The most used algorithms are based on the gradient of the functional with
respect to the field [39, 40, 41], which can be used to ensure a fast
monotonic increase in the objective. This is the case of Krotov’s method [42,
43, 44, 35, 45, 46], as well as of many algorithms developed by Rabitz and
coworkers [47, 48, 49, 50, 51, 52, 53] There is continuously ongoing research
and improvements on Krotov’s method [54, 55, 56], being of particular interest
the algorithms which provide accelerated convergence of quantum control [49,
54] or strategies to smooth the control parameters of the driving field [34].
In this context, there are also studies on the landscape, i. e., topological
and geometric structure defined by the functional [57, 15, 58].

The Quantum Optimal Control Theory (QOCT) is a well tested method that has
proven to be efficient and accurate to control the dynamics of a quantum
system [43]. In this work we revisit the underlying theory to rigorously
obtain the control variational equations. To do so, we also derive the Euler-
Lagrange equations explaining in detail and pedagogically the theory of
variations. The latter can be of interests for many physicist or chemist as
well as grade students who do not fully understand the theory of variations
and can benefit of it to, for instance, develop numerical methods.

This paper is structured as follows: In Section 2, we apply the Euler-Lagrange
equations to obtain a set of QOCT equations, also known as equations of
motion, which differ in a term with the equations obtained in the literature
[1, 43, 39]. In Section 3, we present the necessary conditions to obtain
continuous solutions and demonstrate that the previously labeled ”canonical
solutions” are discontinuous in the costate, but not in the wavefucntion or
the driving field under certain conditions. Our concluding remarks are
summarized in Section 4. Additionally, in the appendix, we provide a brief
introduction to theory of variations, derivations of the Euler-Lagrange
equations, and explain how to impose constraints using a Lagrange multipliers.
We also include some proofs relevant to QOCT mentioned in the primary text.

##  2 Quantum Optimal control

Quantum optimal control is a family of methods to optimize a specific
observable for a quantum system. This involves imposing a cost, which is a
constraint or condition on the field driving the system. In this section, we
derive the QOCT equations by utilizing the Euler-Lagrange equations (ELE). The
ELE are derived using the theory of variations, and their solutions correspond
to the stationary functions. Specifically, the ELE corresponding to the
functional

$\int_{\Omega}\mathcal{F}[y,y_{x_{1}},y_{x_{2}},\ldots,y_{x_{n}},x_{1},x_{2},\ldots,x_{n},\Omega^{\prime}]\mathrm{d}\Omega^{\prime}$
(1)

are given by

$\dfrac{\partial}{\partial
y}\mathcal{F}+\sum_{k=1}^{m}(-1)^{k}\sum\limits_{j=1}^{n}\dfrac{\mathrm{d^{k}}}{\mathrm{d}x^{k}_{j}}\dfrac{\partial}{\partial
y_{x^{k}_{j}}}\mathcal{F}=0,$ (2)

being
$y_{x_{j}^{m}}(x_{1},x_{2},\ldots,x_{n})\equiv\frac{\partial^{m}}{\partial{x_{j}}^{m}}y$.
The ELE are derived in A and B.

Our objective is to optimize the expectation value of the observable $O$ while
minimizing a cost functional and satisfying the TDSE. To achieve this, we
define a functional $J$ and then its extremal points. The functional $J$ may
be divided into three terms:

$\displaystyle
J[\psi,\psi^{*},\chi,\chi^{*},\varepsilon]=J_{\text{opt}}[\psi,\psi^{*}]+J_{\text{cost}}[\varepsilon]+J_{\text{TDSE}}[\psi,\psi^{*},\chi,\chi^{*},\varepsilon].$
(3)

First, the functional corresponding to optimizing the expectation value of the
operator $O$

$J_{\text{opt}}[\psi,\psi^{*}]=\int_{\Omega}\psi^{*}(\Omega,T)O\psi(\Omega,T)\mathrm{d}\Omega,$
(4)

where $\psi(\Omega,t)$ is the wavefunction which describes the quantum system
and $\Omega$ the degrees of freedom (Euler angles, cartesian coordinates of
one or many-body, …). Then, we set the _cost_ , i. e., a condition on the
field, as [39]

$J_{\text{cost}}[\varepsilon]=-\alpha\int\limits_{0}^{T}\left[\varepsilon(t)-\varepsilon_{\text{ref}}(t)\right]^{2}\mathrm{d}t,$
(5)

where $\alpha$ is a penalty factor. It is worth noting that
$J_{\text{cost}}[\varepsilon]$ may be more complex by including a mask
function [39]. However, for the sake of simplicity and without loss of
generality, we will consider this simpler case. Finally, the functional which
ensures that the TDSE is fulfilled reads as

$J_{\text{TDSE}}[\psi,\psi^{*},\chi,\chi^{*},\varepsilon]=-2\operatorname{Im}\int\limits_{0}^{\widehat{T}}\left[\int_{\Omega}\chi^{*}(\Omega,t)\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)\mathrm{d}\Omega\right]\mathrm{d}t,$
(6)

where $\chi(\Omega,t)$ is the Lagrange multiplier. In this work, we set the
Hamiltonian to $H=-\frac{1}{2}\nabla^{2}+V(\Omega,\varepsilon(t))$, where
$V(\Omega,\varepsilon(t))$ is the internal potential plus the interaction with
a driving field, $\varepsilon(t)$. In contrast to the literature, we take the
upper limit of the time interval to $\widehat{T}$, a value greater than the
time at which the function to optimize is evaluated, $T$, denoted as the
_measurement time_ [43]. Specifically, we choose $\widehat{T}>T$, as shown in
Fig.1, to effectively manage the discontinuity caused by the Dirac delta in
the QOCT equations, as described below.

Figure 1: Sketch of a time interval illustrating the definition of $T\text{
and }\widehat{T}$.

To obtain the extremal points we derive the corresponding set of Euler-
Lagrange equations (2). In the case of $\varepsilon(t)$, we find

$\displaystyle
2\operatorname{Im}\int_{\Omega}\chi^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\psi(\Omega,t)\mathrm{d}\Omega-2\alpha\left[(\varepsilon(t)-\varepsilon_{\textup{ref}}(t)\right]=0\Rightarrow$
(7)
$\displaystyle\varepsilon(t)=\varepsilon_{\textup{ref}}(t)+\frac{1}{\alpha}\operatorname{Im}\int_{\Omega}\chi^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}V(\vec{r},\varepsilon)\psi(\Omega,t)\mathrm{d}\Omega.$
(8)

Now, by taking variations on the wavefunction $\psi(\Omega,t)$, we determine
the Lagrange multiplier $\chi^{*}(\Omega,t)$. Note that only the part of the
functional (6) that has a direct dependence on $\psi(\Omega,t)$ is relevant.
As $\psi(\Omega,t)$ and $\psi^{*}(\Omega,t)$ are independent of each other
(for further information, refer to C), the relevant term reads as

$j_{\psi}[\chi^{*},\varepsilon,\psi,\Omega,t]=i\chi^{*}(\Omega,t)\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)+\psi^{*}(\Omega,t)O\psi(\Omega,t)\delta(t-T).$
(9)

Then, the terms in the ELE are

$\displaystyle\dfrac{\partial j_{\psi}}{\partial\psi}$ $\displaystyle=$
$\displaystyle
i\chi^{*}(\Omega,t)\left[-V(\Omega,\varepsilon)\right]+O\psi^{*}(\Omega,t)\delta(t-T),$
(10) $\displaystyle\dfrac{\partial j_{\psi}}{\partial\dot{\psi}}$
$\displaystyle=$
$\displaystyle-\chi^{*}(\Omega,t),\quad\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\partial
j_{\psi}}{\partial\dot{\psi}}=-\dot{\chi}^{*}(\Omega,t),$ (11)
$\displaystyle\dfrac{\partial j_{\psi}}{\partial\psi_{x_{j}}}$
$\displaystyle=$ $\displaystyle 0,\quad\dfrac{\partial
j_{\psi}}{\partial\psi_{x_{j}^{2}}}=\cfrac{i}{2}\chi^{*}(\Omega,t),\quad\dfrac{\mathrm{d}^{2}}{\mathrm{d}x_{j}^{2}}\dfrac{\partial
j_{\psi}}{\partial\psi_{x_{j}^{2}}}=\cfrac{i}{2}\chi_{x_{j}^{2}}^{*}(\Omega,t)~{}\text{for
all}~{}x_{j},$ (12)

where we use the notation $\dot{\psi}(\Omega,t)\equiv\frac{\partial}{\partial
t}\psi(\Omega,t)$ and
${\psi}_{x_{j}^{n}}(\Omega,t)\equiv\frac{\partial^{n}}{\partial
x_{j}^{n}}\psi(\Omega,t)$. Note that we have used the identity
$\int_{\Omega}\psi^{*}(\Omega,T)O\psi(\Omega,T)\mathrm{d}\Omega=\int_{\Omega}\int_{0}^{\widehat{T}}\psi^{*}(\Omega,t)O\psi(\Omega,t)\delta(t-T)\mathrm{d}\Omega\mathrm{d}t$
and Eq. (2) to derive the contributions to the ELE with high order
derivatives. After summing all the terms, we obtain

$\left(i\dfrac{\partial}{\partial
t}+H(\Omega,\varepsilon)\right)\chi^{*}(\Omega,t)=-iO\psi^{*}(\Omega,t)\delta(t-T).$
(13)

To obtain the ELE for $\chi(\Omega,t)$ we can either take variations on
$\psi^{*}(\Omega,t)$ or simply take the complex conjugate of Eq. (13),

$\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\chi(\Omega,t)=-iO\psi(\Omega,t)\delta(t-T).$
(14)

The TDSE of $\psi^{*}(\Omega,t)$ is obtained by taking variations of
$\chi(\Omega,t)$. The relevant term reads

$j_{\chi}[\psi^{*},\chi,\varepsilon,\Omega,t]=-i\chi(\Omega,t)\left(i\frac{\partial}{\partial
t}+H(\Omega,\varepsilon)\right)\psi^{*}(\Omega,t)$ (15)

Each term in the ELE is

$\displaystyle\dfrac{\partial j_{\chi}}{\partial\dot{\chi}}$ $\displaystyle=$
$\displaystyle 0,\qquad\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\partial
j_{\chi}}{\partial\dot{\chi}}=0$ (16) $\displaystyle\dfrac{\partial
j_{\chi}}{\partial\chi}$ $\displaystyle=$
$\displaystyle-i\left(i\frac{\partial}{\partial
t}+H(\Omega,\varepsilon(t))\right)\psi^{*}(\Omega,t).$ (17)

By combining all the terms we derive the EOM

$\left(i\frac{\partial}{\partial
t}+H(\Omega,\varepsilon)\right)\psi^{*}(\Omega,t)=0$ (18)

Finally, we take variations with respect to the Lagrange multiplier
$\chi^{*}(\Omega,t)$. As done above, we work with the relevant term

$\displaystyle j_{\chi^{*}}[\psi,\chi^{*},\varepsilon,\Omega,t]$
$\displaystyle=$ $\displaystyle
i\chi^{*}(\Omega,t)\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t).$ (19)

The terms in the associated ELE are

$\displaystyle\dfrac{\partial j_{\chi}}{\partial\dot{\chi}^{*}}$
$\displaystyle=$ $\displaystyle
0,\qquad\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\partial
j_{\chi^{*}}}{\partial\dot{\chi}^{*}}=0,$ (20) $\displaystyle\dfrac{\partial
j_{\chi^{*}}}{\partial\chi^{*}}$ $\displaystyle=$ $\displaystyle
i\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)\ ,$ (21)

from which we obtain the TDSE for $\psi(\Omega,t)$

$\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)=0.$ (22)

As expected, Eq. (22) is the complex conjugate of Eq. (18), which ensures that
the functional $J[\psi,\psi^{*},\chi,\chi^{*},\varepsilon]$ is appropriate.
Summing up, we derived the set of equations

$\displaystyle\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)=0,$ (23)
$\displaystyle\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\chi(\Omega,t)=-iO\psi(\Omega,t)\delta(t-T),$
(24)
$\displaystyle\varepsilon(t)=\varepsilon_{\textup{ref}}(t)+\frac{1}{\alpha}\operatorname{Im}\int_{\Omega}\chi^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\psi(\Omega,t)\mathrm{d}\Omega,$
(25)

and their complex conjugates.

In the following sections we obtain two different types of solutions, namely,
continuous and discontinuous solutions. The latter are related to the
_canonical_ QOCT equations.

##  3 Extremal trajectories: Conditions and continuity

###  3.1 Continuous solutions

In Eq.(23), the continuity of $\psi(\Omega,t)$ is guaranteed, while the
continuity of $\chi(\Omega,t)$ and thus $\varepsilon(t)$ cannot be assured by
Eq.(24). Here, we show that $\chi(\Omega,t)$ is continuous if it fulfills
$\chi(\Omega,t)=\frac{i}{2\pi n}O\psi(\Omega,t)$, where $n\in\mathds{Z}-{0}$.

###### Proof.

Let us assume the ansatz
$\chi(\Omega,t)=\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}$, where
$\widetilde{\chi}(\Omega,t)$ is continuous and fullfils the TDSE,
$\phi\in\mathds{C}$ and $\Theta(t-T)$ is the Heaviside function. We substitute
the ansatz in Eq. (24)

$\displaystyle i\frac{\partial}{\partial
t}\left(\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}\right)-\hat{H}(\Omega,t)\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}=-iO\psi(\Omega,t)\delta(t-T)$
$\displaystyle ie^{i\phi\Theta(t-T)}\frac{\partial}{\partial
t}\widetilde{\chi}(\Omega,t)+i\widetilde{\chi}(\Omega,t)\frac{\partial}{\partial
t}e^{i\phi\Theta(t-T)}-\hat{H}(\Omega,t)\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}=-iO\psi(\Omega,t)\delta(t-T)$
$\displaystyle e^{i\phi\Theta(t-T)}\left(i\frac{\partial}{\partial
t}\widetilde{\chi}(\Omega,t)-\hat{H}(\Omega,t)\widetilde{\chi}(\Omega,t)\right)=\left(\phi\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}-iO\psi(\Omega,t)\right)\delta(t-T)$

To remove the discontinuity, we may set

$\lim\limits_{t\rightarrow
T}\phi\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}=iO\psi(\Omega,T)\Rightarrow\lim\limits_{t\rightarrow
T}\chi(\Omega,t)=\frac{i}{\phi}O\psi(\Omega,T)$ (26)

Finally, $\lim\limits_{t\rightarrow
T}\widetilde{\chi}(\Omega,t)e^{i\phi\Theta(t-T)}$ is defined if
$e^{i\phi\Theta(t-T)}$ is continuous, i. e., $\phi=2\pi n$ with
$n\in\mathds{Z}-\\{0\\}$. Therefore, $\chi(\Omega,t)$ is continuous if and
only if $\chi(\Omega,T)=\frac{i}{2\pi n}O\psi(\Omega,T)$. Note that it is easy
to prove that $\epsilon(t)$ is also continuous since $\psi(\Omega,t)$ and
$\chi(\Omega,t)$ are continuous. ∎

It should be noted that the condition of continuity is imposed at a specific
time $T$, which poses a challenge as $T$ is not a boundary of the time
interval (as illustrated in Fig. 2). Therefore, such solutions may not be
appropriate from a formal standpoint, especially when extending the system to
a finite _measurement time_. Let us also remark that the solution here
obtained differs from the costate provided in the literature [43, 39, 46]. In
the following section, we will provide evidence to demonstrate that the
costate reported in earlier studies corresponds to a solution that is
discontinuous.

###  3.2 Relationship to _Canonical_ quantum optimal control equations

In this section, we introduce a less restrictive condition for
$\chi(\Omega,t)$ at $t=T$, allowing for discontinuity, and impose an initial
condition such that $\chi(\Omega,T^{\prime})=0$ for
$T^{\prime}\textgreater{T}$.

###### Proof.

Let us define $\chi(\Omega,t)=\widetilde{\chi}(\Omega,t)f(t)$, where
$\widetilde{\chi}(\Omega,t)$ satisfies the TDSE
$\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\widetilde{\chi}(\Omega,t)=0$ and
$\widetilde{\chi}(\Omega,t)=\chi(\Omega,t)$ for $t\in[0,T)$. We plug it in Eq.
(14)

$\displaystyle\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\left(\widetilde{\chi}(\Omega,t)f(t)\right)=-iO\psi(\Omega,t)\delta(t-T)$
$\displaystyle i\widetilde{\chi}(\Omega,t)\dfrac{\partial}{\partial
t}f(t)+f(t)\left(i\dfrac{\partial}{\partial
t}\widetilde{\chi}(\Omega,t)-H(\Omega,\varepsilon)\widetilde{\chi}(\Omega,t)\right)=-iO\psi(\Omega,t)\delta(t-T)$
$\displaystyle\widetilde{\chi}(\Omega,t)\dfrac{\partial}{\partial
t}f(t)=-O\psi(\Omega,t)\delta(t-T)$ (27)

Next, we integrate over time in the interval $[T-\delta,T+\delta]$ with
$\delta>0$ and take the limit $\delta\rightarrow 0$

$\lim_{\delta\rightarrow 0}f(\Omega,T+\delta)-\lim_{\delta\rightarrow
0}f(\Omega,T-\delta)=-\cfrac{O\psi(\Omega,T)}{\widetilde{\chi}(\Omega,T)}.$
(28)

Setting $\chi(\Omega,t)=0$ for $t>T$ and using that
$\chi(\Omega,t)=\widetilde{\chi}(\Omega,t)$ for $t<T$, that is to say,
$\lim_{\delta\rightarrow 0}f(T+\delta)=0$ and $\lim_{\delta\rightarrow
0}f(T-\delta)=1$, we obtain

$-1=-\cfrac{O\psi(\Omega,T)}{\widetilde{\chi}(\Omega,T)}\Rightarrow\widetilde{\chi}(\Omega,T)=O\psi(\Omega,T),$
(29)

therefore, we obtain the boundary condition for the costate

$\lim_{\delta\rightarrow 0}\chi(\Omega,T-\delta)=O\psi(\Omega,T).$ (30)

Next, we plug Eq. (30) in Eq. (25)

$\displaystyle\lim_{\delta\rightarrow 0}\varepsilon(T-\delta)$
$\displaystyle=$
$\displaystyle\varepsilon_{\textup{ref}}(T)+\frac{1}{\alpha}\operatorname{Im}\int_{\Omega}O\psi^{*}(\Omega,T)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\psi(\Omega,T)\mathrm{d}\Omega$
(31) $\displaystyle=$
$\displaystyle\varepsilon_{\textup{ref}}(T)+\frac{1}{\alpha}\operatorname{Im}\left\langle\psi\left|O\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\right|\psi\right\rangle$

Note that the integral in Eq. (31) is real if
$O\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)$ is Hermitian,
that is,
$\left[O,\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\right]=0$.
As a consequence, we have $\lim\limits_{\delta\rightarrow
0}\varepsilon(T-\delta)=\varepsilon_{\text{ref}}(T)$. Moreover, even though
$\chi(\Omega,t)$ is not continuous at $t=T$, we have that
$\lim_{\delta\rightarrow 0}\varepsilon(T+\delta)=\varepsilon_{\text{ref}}(T)$,
which ensures the continuity of $\varepsilon(t)$ at $t=T$. ∎

##  4 Conclusions

In this study, we have presented a comprehensive explanation of the theory of
variations that is essential to derive the Quantum Optimal Control Theory
(QOCT) equations. Specifically, we have derived the Euler-Lagrange equations
for multiple variables and demonstrated a systematic approach to incorporate
constraints in any multidimensional functional by utilizing the Lagrange
multipliers.

Moreover, we have shown that the _canonical_ solutions of the QOCT equations,
which are commonly reported in the literature, exhibit discontinuity at the
measurement time, as it is commonly assumed. In contrast, we have identified a
new set of continuous solutions of QOCT equations. However, the lack of well-
defined boundary conditions limits their practical applicability.

Finally, let us highlight that our work presents a foundation for developing
novel optimal control strategies using our methodology.

J.J.O. gratefully acknowledges the funding from the Madrid Government
(Comunidad de Madrid Spain) under the Multiannual Agreement with Universidad
Complutense de Madrid in the line Research Incentive for Young PhDs, in the
context of the V PRICIT (Regional Programme of Research and Technological
Innovation) (Grant: PR27/21-010), Project PID2019-106732GB-I00 (MINECO) and
the Social European Fund and Juan de la Cierva-Incorporación granted by the
Ministerio de Ciencia e Innovación. I.R.S. acknowledges support from MINECO
PID2021-122796NB-I00. K.C. acknowledges support form MINECO
PID2020-118180GB-I00, and Junta de Andalucía (Spain) grants A-FQM-441-UGR18
and P20-00164.

##  Appendix A Theory of variations

As we have sketched above, the QOCT is a widely used procedure to maximize or
minimize an observable of a quantum system by setting several constraints.
Most of the theory developed in the literature addresses physicists and
chemists, and lacks a rigorous mathematical framework, which demands clear and
correct proofs. In particular, this is due to a misunderstanding in the
concept of  _variation_ , which is defined in mathematics as  _any function
$\eta(\Omega)$ that perturbs a given function $f(\Omega)$, fulfilling that
$\eta(\Omega)$ vanishes in the border of the domain of $f(\Omega)$_. Let us
also note that this vanishing condition may also apply for the derivatives of
$\eta(\Omega)$, as we explain in detail in B.

In the following sections we show and derive the mathematical tools to
correctly derive the EOM of QOC. We first derive the Euler-Lagrange equations
(ELE), which are the equations fulfilled by the functions which optimize a
given functional. By applying the ELE, we do not need to make explicit use of
the variations, facilitating and clarifying the derivation of our set of QOCT
equations. Next, we explain in detail how to add constraints to the
functional, such as the total energy of the driving field.

###  A.1 Derivation of Euler-Lagrange equations

The minimization principle consists on seeking for the minimum of the
functional

$\int_{\Omega}\mathcal{F}[\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},\ldots,\widehat{y}_{x_{n}},\Omega^{\prime}]\mathrm{d}\Omega^{\prime},$
(32)

where $\mathcal{F}$ is a functional of $\widehat{y}(x_{1},x_{2},\ldots,x_{n})$
defined in a $n$ dimensional domain $\Omega$, its first derivatives and
$\Omega^{\prime}=(x_{1},x_{2},\ldots,x_{n})$. Note that we use the notation
$y_{x_{i}}=\frac{\partial y}{\partial x_{i}}$. The Euler-Lagrange equations
(ELE) are the conditions of the extremal functions, that is, the conditions
that have to be satisfied by $\\{y_{x_{j}}\\}_{j=1}^{n}$, that optimize the
functional (32). They read as

$\displaystyle\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\ldots,\widehat{y}_{x_{n}},x_{1},\ldots,x_{n})-$
$\displaystyle\sum\limits_{j=1}^{n}\dfrac{\mathrm{d}}{\mathrm{d}x_{j}}\dfrac{\partial}{\partial\widehat{y}_{x_{j}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\ldots,\widehat{y}_{x_{n}},x_{1},\ldots,x_{n})=0.$
(33)

###### Proof.

To seek for these extremal functions we follow the same strategy than Lagrange
[2]. We take

$\widehat{y}(x_{1},x_{2},\ldots,x_{n},\epsilon)=y(x_{1},x_{2},\ldots,x_{n})+\xi\eta(x_{1},x_{2},\ldots,x_{n}),$
(34)

where $\eta(x_{1},x_{2},\ldots,x_{n})\in\mathds{C}$ is such that vanishes on
$\partial\Omega$, i. e., the border of the domain $\Omega$ and
$\epsilon\in\mathds{R}$. Let us note that we may find in the literature that
the  _variation_ is defined as $\delta
y(x_{1},x_{2},\ldots,x_{n})\equiv\xi\eta(x_{1},x_{2},\ldots,x_{n})$. Thus,
$\widehat{y}(x_{1},x_{2},\ldots,x_{n},\epsilon)$ and
$y(x_{1},x_{2},\ldots,x_{n})$ have the same value in the $\partial\Omega$ for
all $\eta(x_{1},x_{2},\ldots,x_{n})$. In Fig. 2, we illustrate how the
variations affect a curve in a plane with the calculation of the shortest
distance between two points. For the sake of simplicity we restrict our
analysis to the 2 dimensional case, which can be easily generalized. The
stationary function is a geodesic which corresponds to the straight line
linking two points. Therefore, the problem is reduced to finding
$y(x_{1},x_{2},\ldots,x_{n})$ so that

$\left.\dfrac{\partial}{\partial\xi}\int_{\Omega}\mathcal{F}[\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},\Omega^{\prime}]\mathrm{d}\Omega^{\prime}\right|_{\xi=0}=0.$
(35)

Figure 2:  Paths linking $y_{\text{initial}}$ and $y_{\text{end}}$. The
function $y(t)$ (solid red line) represents the extremal function and the
other curves correspond to this extremal solution plus a variation.

We assume that $y(x_{1},x_{2})$ is a stationary point of
$\int_{\Omega}\mathcal{F}[y,y_{x_{1}},y_{x_{2}},x_{1},x_{2}]\mathrm{d}x_{1}\mathrm{d}x_{2}$.
Then, we study the functional close to this point by adding a  _variation_
$\delta y=\xi\eta(x_{1},x_{2})$, where $\eta(x_{1},x_{2})$ is any function
which vanishes at the border of the domain of $y(x_{1},x_{2})$ (see Fig. 2).
The variation ensures that for $\xi=0$,
$\int_{\Omega}\mathcal{F}[y+\xi\eta,y_{x_{1}}+\xi\eta_{x_{1}},y_{x_{2}}+\xi\eta_{x_{2}},x_{1},x_{2}]\mathrm{d}x_{1}\mathrm{d}x_{2}$
has an extremal, i. e., it fulfills the condition

$\left.\dfrac{\partial}{\partial\xi}\int_{\Omega}\mathcal{F}[\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2}]\mathrm{d}x_{1}\mathrm{d}x_{2}\right|_{\xi=0}=0.$
(36)

Expanding the derivative, we obtain

$\displaystyle\int\limits_{\Omega}\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\dfrac{\partial\widehat{y}}{\partial\xi}+\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\dfrac{\partial\widehat{y}_{x_{1}}}{\partial\xi}+\right.$
(37)
$\displaystyle\left.\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\dfrac{\partial\widehat{y}_{x_{2}}}{\partial\xi}\right]\mathrm{d}x_{1}\mathrm{d}x_{2}=$
(38)
$\displaystyle\int\limits_{\Omega}\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})+\right.$
(39)
$\displaystyle\left.\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta_{x_{1}}(x_{1},x_{2})+\right.$
(40)
$\displaystyle\left.\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta_{x_{2}}(x_{1},x_{2})\right]\mathrm{d}x_{1}\mathrm{d}x_{2}=$
(41)
$\displaystyle\int\limits_{\Omega}\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}+$
(42)
$\displaystyle\left[\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})\right]_{\partial\Omega}$
(43)
$\displaystyle-\int\limits_{\Omega}\mathrm{d}x_{1}\mathrm{d}x_{2}\dfrac{\mathrm{d}}{\mathrm{d}x_{1}}\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})+$
(44)
$\displaystyle\left[\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})\right]_{\partial\Omega}$
(45)
$\displaystyle-\int\limits_{\Omega}\dfrac{\mathrm{d}}{\mathrm{d}x_{2}}\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\eta(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}.$
(46)

The terms (43) and (45) vanish because $\eta(x_{1},x_{2})=0$ in the border of
the domain, i. e., for $(x_{1},x_{2})\in\partial\Omega$. Therefore, the
condition for the stationary point is

$\displaystyle\int\limits_{\Omega}\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})-\dfrac{\mathrm{d}}{\mathrm{d}x_{1}}\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\right.$
$\displaystyle\left.-\dfrac{\mathrm{d}}{\mathrm{d}x_{2}}\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})\right]\eta(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}=0.$
(47)

Since this condition has to be fulfilled for all $\eta(x_{1},x_{2})$

$\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})-\dfrac{\mathrm{d}}{\mathrm{d}x_{1}}\dfrac{\partial}{\partial\widehat{y}_{x_{1}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})-\dfrac{\mathrm{d}}{\mathrm{d}x_{2}}\dfrac{\partial}{\partial\widehat{y}_{x_{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x_{1}},\widehat{y}_{x_{2}},x_{1},x_{2})=0.$
(48)

This result can be generalized for more than two variables following the same
procedure, obtaining

$\dfrac{\partial}{\partial
y}\mathcal{F}(y,y_{x_{1}},\ldots,y_{x_{n}},x_{1},\ldots,x_{n})-\sum\limits_{j=1}^{n}\dfrac{\mathrm{d}}{\mathrm{d}x_{j}}\dfrac{\partial}{\partial
y_{x_{j}}}\mathcal{F}(y,y_{x_{1}},\ldots,y_{x_{n}},x_{1},\ldots,x_{n})=0,$
(49)

where we have used that $\widehat{y}(\Omega,t)=y(\Omega,t)$ for $\xi=0$. ∎

It is important to remark that it is not necessary to know $y$ in the entire
border if the ELE guarantee the existence and uniqueness of the solution.

Finally, let us now consider a functional depending on a function of high
order derivatives as
$\mathcal{F}(\widehat{y},\\{\widehat{y}_{x_{j}}\\}_{j=1}^{n},\\{\widehat{y}_{x_{j}^{2}}\\}_{j=1}^{n},\ldots,\\{\widehat{y}_{x_{j}^{m}}\\}_{j=1}^{n},x_{1},\ldots,x_{n})$,
where $y_{x_{j}^{k}}\equiv\cfrac{\partial^{k}y}{\partial x^{k}_{j}}$. Thus,
the ELE for functional with higher order derivatives read

$\dfrac{\partial}{\partial
y}\mathcal{F}+\sum_{k=1}^{m}(-1)^{k}\sum\limits_{j=1}^{n}\dfrac{\mathrm{d^{k}}}{\mathrm{d}x^{k}_{j}}\dfrac{\partial}{\partial
y_{x^{k}_{j}}}\mathcal{F}=0,$ (50)

where we have not explicitly indicated the arguments of the functional for the
sake of clarity. A detailed proof can be found in B.

###  A.2 Adding constraints

Adding a constraint may be mandatory to solve a given problem, and it is done
by using the Lagrange multipliers [3]. In the context of calculus, to find the
stationary point of a function $y(x_{1},x_{2})$ fulfilling the condition
$g(x_{1},x_{2})=0$ we have to find the stationary point of

$f(x_{1},x_{2},\lambda)=y(x_{1},x_{2})+\lambda g(x_{1},x_{2}).$ (51)

The stationary points fulfill

$\dfrac{\partial}{\partial
x_{1}}f(x_{1},x_{2},\lambda)=\dfrac{\partial}{\partial
x_{2}}f(x_{1},x_{2},\lambda)=\dfrac{\partial}{\partial\lambda}f(x_{1},x_{2},\lambda)=0,$
(52)

where $\lambda$ is the  _Lagrange_ multiplier and the last equality ensures
that the constraint $g(x_{1},x_{2})=0$.

The extension to the theory of variations consists on taking a function as
Lagrange multiplier. In this case, we restrict to functionals with first
derivatives only. The generalization can be done following the steps described
in B. The multiplier function, $\lambda(x_{1},x_{2},\ldots,x_{n})$, may depend
on an arbitrary number of variables, depending on the constraint. Let us show
how to construct the functionals based on these constrains with some
illustrative examples in a two dimensional space.

####  A.2.1 Constraint at each $(x_{1},x_{2})\in\Omega$

First, we consider a constraint at each $(x_{1},x_{2})\in\Omega$, i. e.,
$g(y,y_{x_{1}},y_{x_{2}},\ldots,x_{1},x_{2})=0$. Then, we sum the term

$\int\limits_{\Omega}\lambda
g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}\Omega.$ (53)

to the functional to be minimized,
$\int\limits_{\Omega}\mathcal{F}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}\Omega$.
Now, we apply the ELE (48) to $\lambda(x_{1},x_{2})$ obtaining
$g(y,y_{x_{1}},y_{x_{2}},\ldots,x_{1},x_{2})=0$. Note that we are not taking
into account
$\int\limits_{\Omega}\mathcal{F}(y,y_{x_{1}},y_{x_{2}},\ldots,x_{1},x_{2})\mathrm{d}\Omega$
here, since it does not depend on $\lambda$. Finally, in order to obtain the
ELE corresponding to $y$’s we use the functional with constrains

$\displaystyle\mathcal{G}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2},\lambda,\epsilon,\eta)=$
$\displaystyle\int\limits_{\Omega}\left[\mathcal{F}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})+\lambda
g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\right]\mathrm{d}\Omega$ (54)

####  A.2.2 Constraint given by an integral in one variable of a function of
two variables

Now we analyze another remarkable kind of constrains, given by _an integral in
one variable of a function of two variables_ , i. e.,
$\int_{\Omega_{1}}g(y,\\{y_{x_{i}}\\},\\{y_{x_{i},x_{j}}\\},\ldots,x_{1},x_{2})\mathrm{d}x_{1}=0$.
For example, the conservation of the norm of the wavefunction (integration in
the position space of the norm squared) for all time $t$ is commonly applied
to many-body physics such as Multiconfigurational Time-Dependent Hartree-Fock
and related self-consistent methods [11, 12, 10].

In this case we include the functional

$\int\limits_{\Omega_{1}}\int\limits_{\Omega_{2}}\lambda(x_{2})g(y,y_{x_{1}},y_{x_{1}},x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2},$
(55)

where we assume that $\Omega$ can be split in $\Omega_{1}$ and $\Omega_{2}$,
with $x_{j}\in\Omega_{j}$. Here we can not apply directly the ELE as it is
expressed above, since the multiplier depends on less variables than the
dimensions of $\Omega$. However, we can still use Lagrange’s methodology by
adding a variation $\xi\eta(x_{2})$ to the multiplier and deriving with
respect to $\xi$.

$\displaystyle\mathcal{G}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2},\xi,\eta)$
$\displaystyle=$
$\displaystyle\int\limits_{\Omega_{2}}\int\limits_{\Omega_{1}}(\lambda(x_{2})+\xi\eta(x_{2}))g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2},$

and then, deriving with respect to $\xi$

$\displaystyle\dfrac{\partial}{\partial\xi}\mathcal{G}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2},\xi,\eta)$
$\displaystyle=$
$\displaystyle\int\limits_{\Omega_{2}}\int\limits_{\Omega_{1}}\eta(x_{2})g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}$
(56) $\displaystyle=$
$\displaystyle\int\limits_{\Omega_{2}}\eta(x_{2})\int\limits_{\Omega_{1}}g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}.$
(57)

Thus,
$\dfrac{\partial}{\partial\xi}\mathcal{G}(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2},\xi,\eta)=0$
for all $\eta(x_{2})$ if and only if the constrain
$\int\limits_{\Omega_{1}}g(y,y_{x_{1}},y_{x_{2}},x_{1},x_{2})\mathrm{d}x_{1}=0$
is fulfilled.

####  A.2.3 General case: Constraint given by an integral in $m$ variables of
a function of $n$ variables

The strategy consists on including a Lagrange multiplier in the $n-m$
variables which are not integrated, as can be generalized from the previous
cases.

##  Appendix B Derivation of Euler-Lagrange type equation for higher order
derivatives

In this section we derive the ELE for a functional with higher order
derivatives. First, we will derive the case of second and first order
derivatives and one variable. Then, we will generalize the result.

###### Proof.

Let us take the functional
$\int\limits_{\Omega}\mathcal{F}(y,y_{x},y_{x^{2}},x)\mathrm{d}x$. We follow a
similar procedure as in A.1 to compute the ELE. First, we evaluate the
functional with the function $\widehat{y}(x)=y(x)+\xi\eta(x)$, where $y(x)$ is
the extremal, $\eta(x)$ is a function such that
$\eta(x)|_{\partial\Omega}=\eta_{x}(x)|_{\partial\Omega}=0$ and
$\xi\in\mathds{R}$

$\displaystyle\mathcal{I}[(y+\xi\eta,y_{x}+\xi\eta_{x},y_{x^{2}}+\xi\eta_{x^{2}})]=$
$\displaystyle\int\limits_{\Omega}\mathcal{F}(y+\xi\eta,y_{x}+\xi\eta_{x},y_{x^{2}}+\xi\eta_{x^{2}},x)\mathrm{d}x.$
(58)

Since $y(x)$ is an extremal of the functional $\mathcal{I}$, then
$\frac{\mathrm{d}}{\mathrm{d}\xi}\mathcal{I}[(y+\xi\eta,y_{x}+\xi\eta_{x},y_{x^{2}}+\xi\eta_{x^{2}})]|_{\xi=0}=0$.
The derivative reads as

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}\xi}\mathcal{I}[(y+\xi\eta,y_{x}+\xi\eta_{x},y_{x^{2}}+\xi\eta_{x^{2}})]=$
$\displaystyle\int\limits_{\Omega}\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x}^{2},x)\eta(x)+\dfrac{\partial}{\partial\widehat{y}_{x}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x}^{2},x)\eta_{x}(x)\right.$
$\displaystyle\dfrac{\partial}{\partial\widehat{y}_{x}^{2}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)\eta_{x^{2}}(x)\left.\right]\mathrm{d}x,$
(59)

where we have used that $\widehat{y}(x)=y(x)+\xi\eta(x)$ and therefore
$\frac{\mathrm{d}}{\mathrm{d}\xi}\widehat{y}(x)=\eta(x)$. Now, we integrate by
parts to remove any derivative of $\eta(x)$

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}\xi}\mathcal{I}[(y+\xi\eta,y_{x}+\xi\eta_{x},y_{x^{2}}+\xi\eta_{x^{2}})]=$
$\displaystyle\int\limits_{\Omega}\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)-\dfrac{\mathrm{d}}{\mathrm{d}x}\dfrac{\partial}{\partial\widehat{y}_{x}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)+\right.$
$\displaystyle\left.\dfrac{\mathrm{d^{2}}}{\mathrm{d}x^{2}}\dfrac{\partial}{\partial\widehat{y}_{x^{2}}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)\right]\eta(x)\mathrm{d}x+$
$\displaystyle\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)\eta(x)\right]_{\partial\Omega}+\left[\dfrac{\partial}{\partial\widehat{y}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)\eta_{x}(x)\right]_{\partial\Omega}$
$\displaystyle-\left[\dfrac{\mathrm{d}}{\mathrm{d}x}\dfrac{\partial^{2}}{\partial\widehat{y}_{x}^{2}}\mathcal{F}(\widehat{y},\widehat{y}_{x},\widehat{y}_{x^{2}},x)\eta(x)\right]_{\partial\Omega}.$
(60)

Note that we integrated the last term of Eq. (59) twice by parts to obtain
$\eta(x)$ under the integral. We obtain a plus sign in front of the third term
of the integral. Since
$\eta(x)|_{\partial\Omega}=\eta_{x}(x)|_{\partial\Omega}=0$, the three last
terms in Eq. (60) vanish. Finally, since Eq. (60) must be zero for all
$\eta(x)$, we obtain, for $\xi=0$

$\displaystyle\dfrac{\partial}{\partial
y}\mathcal{F}(y,y_{x},y_{x^{2}},x)-\dfrac{\mathrm{d}}{\mathrm{d}x}\dfrac{\partial}{\partial
y_{x}}\mathcal{F}(y,y_{x},y_{x^{2}},x)+\dfrac{\mathrm{d^{2}}}{\mathrm{d}x^{2}}\dfrac{\partial}{\partial
y_{x^{2}}}\mathcal{F}(y,y_{x},y_{x^{2}},x)=0.$ (61)

This proof can be generalized straightforwardly by imposing that $\eta(x)$ and
its $n-1$ derivatives vanish at $\partial\Omega$ for a functional involving
$n$ order derivatives. Then we get the Euler-Lagrange for a $n$-order
functional

$\displaystyle\dfrac{\partial}{\partial
y}\mathcal{F}(y,y_{x},\ldots,y_{x^{n}},x)+\sum_{k=1}^{n}(-1)^{k}\dfrac{\mathrm{d^{k}}}{\mathrm{d}x^{k}}\dfrac{\partial}{\partial
y_{x^{k}}}\mathcal{F}(y,y_{x},\ldots,y_{x^{n}},x)=0$

∎

##  Appendix C Independence of the complex conjugate in the functional

Deriving one function with respect to its complex conjugate in a functional is
a _confusing_ situation that we physicists have to face from time to time.
Usually we try to bypass this uncomfortable dilemma by seeking for an
alternative way, as in the case of variations which involve the wave function
and its complex conjugate. This cases are specially common in path integral,
QED [59] or many-body theory [60, 10], among others. Below, we show that in
the derivation of the Euler-Lagrange Equations from the functional (3) we can
consider $\Psi(\Omega,t)$ and ${\Psi}^{*}(\Omega,t)$ as independent.

###### Proof.

Let us define $\mathcal{G}(\Phi,\Psi,\chi,\zeta,\varepsilon)$ as

$\displaystyle\mathcal{G}(\Phi,\Psi,\chi,\zeta,\varepsilon)$ $\displaystyle=$
$\displaystyle-i\int\limits_{0}^{\widehat{T}}\left[\int_{\Omega}\Phi(\Omega,t)\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\chi(\Omega,t)\mathrm{d}\Omega\right]\mathrm{d}t$
(62)
$\displaystyle+i\int\limits_{0}^{\widehat{T}}\left[\int_{\Omega}\zeta(\Omega,t)\left(i\frac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\Psi(\Omega,t)\mathrm{d}\Omega\right]\mathrm{d}t$
$\displaystyle-\alpha\int\limits_{0}^{\widehat{T}}\left[\varepsilon(t)-\varepsilon_{\text{ref}}(t)\right]^{2}\mathrm{d}t+\int_{\Omega}\Phi
O\Psi\mathrm{d}\Omega$

Applying the Euler-Lagrange equations for $\chi,\zeta,\Phi,\Psi$ and
$\varepsilon$ we obtain

$\displaystyle i\frac{\partial}{\partial t}\Phi(\Omega,t)$ $\displaystyle+$
$\displaystyle H(\Omega,\varepsilon)\Phi(\Omega,t)=0$ (63) $\displaystyle
i\frac{\partial}{\partial t}\Psi(\Omega,t)$ $\displaystyle-$ $\displaystyle
H(\Omega,\varepsilon)\Psi(\Omega,t)=0$ (64) $\displaystyle
i\frac{\partial}{\partial t}\chi(\Omega,t)$ $\displaystyle-$ $\displaystyle
H(\Omega,\varepsilon)\chi(\Omega,t)=-iO\Psi(\Omega,t)$ (65) $\displaystyle
i\frac{\partial}{\partial t}\zeta(\Omega,t)$ $\displaystyle+$ $\displaystyle
H(\Omega,\varepsilon)\zeta(\Omega,t)=-iO\Phi(\Omega,t)$ (66)
$\displaystyle\varepsilon(t)$ $\displaystyle=$
$\displaystyle\varepsilon_{\text{ref}}(t)-\frac{1}{2\alpha}\left[-i\int_{\Omega}\Phi(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\chi(\Omega,t)\mathrm{d}\Omega\right.$
(67)
$\displaystyle\left.+i\int_{\Omega}\zeta(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\Psi(\Omega,t)\mathrm{d}\Omega\right]$

By inspecting Eqs. (63) and (64) we check that ${\Phi}^{*}(\Omega,t)$ and
$\Psi(\Omega,t)$ fulfil the same equation. Furthermore, if we set
$\Phi(\Omega,0)=\beta{\Psi}^{*}(\Omega,0)$, we find that
$\Phi(\Omega,t)=\beta{\Psi}^{*}(\Omega,t)$, being $\beta\in\mathbb{C}$.
Analogously, we get that $\zeta(\Omega,t)=\xi{\chi}^{*}(\Omega,t)$, where
$\xi\in\mathbb{C}$. Next, we plug this in Eq. (67)

$\displaystyle\varepsilon(t)$ $\displaystyle=$
$\displaystyle\varepsilon_{\text{ref}}(t)-\frac{1}{2\alpha}\left[\beta\int_{\Omega}{\Psi}^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\chi(\Omega,t)\mathrm{d}\Omega\right.$
(68)
$\displaystyle\left.+\xi\int_{\Omega}{\chi}^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\Psi(\Omega,t)\mathrm{d}\Omega\right].$

Note that for $\beta=\xi=1$, we obtain

$\displaystyle\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\psi(\Omega,t)=0,$ (69)
$\displaystyle\left(i\dfrac{\partial}{\partial
t}-H(\Omega,\varepsilon)\right)\chi(\Omega,t)=-iO\psi(\Omega,t)\delta(t-T),$
(70)
$\displaystyle\varepsilon(t)=\varepsilon_{\textup{ref}}(t)+\frac{1}{\alpha}\operatorname{Im}\int_{\Omega}\chi^{*}(\Omega,t)\frac{\partial}{\partial\varepsilon}H(\Omega,\varepsilon)\psi(\Omega,t)\mathrm{d}\Omega,$
(71)

and their complex conjugates.

Summing up, by setting $\Phi(\Omega,0)=\Psi^{*}(\Omega,0)$ and
$\zeta(\Omega,0)=\chi^{*}(\Omega,0)$ we may derive the QOCT equations
(23)-(25) as done in Sec. 2 by assuming that $\Psi(\Omega,t)$ and
$\chi(\Omega,t)$ and their complex conjugates are independent. ∎

## References

  * [1] Anthony P Peirce, Mohammed A Dahleh, and Herschel Rabitz.  Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications.  Phys. Rev. A, 37(12):4950–4964, jun 1988.
  * [2] H H Goldstine.  A history of the calculus of variations from the 17th through the 19th century.  Studies in the history of mathematics and physical sciences. Springer-Verlag, 1980.
  * [3] J Stewart.  Calculus.  Cengage Learning, 2011.
  * [4] H Goldstein, C P Poole, and J L Safko.  Classical Mechanics.  Addison Wesley, 2002.
  * [5] M Mesterton-Gibbons.  A Primer on the Calculus of Variations and Optimal Control Theory.  Student mathematical library. American Mathematical Society, 2009.
  * [6] John David Jackson.  Classical electrodynamics.  Wiley, New York, NY, 3rd ed. edition, 1999.
  * [7] F P Miller, A F Vandome, and J McBrewster.  Einstein-Hilbert Action: General Relativity, Action (physics), Einstein Field Equations, Equations of Motion, David Hilbert, Maxwell’s Equations, Noether’s Theorem.  Alphascript Publishing, 2010.
  * [8] A. Raab.  On the Dirac-Frenkel/McLachlan variational principle.  Chem. Phys. Lett., 319(March):674–678, 2000.
  * [9] H.-D. Meyer, U. Manthe, and L. S. Cederbaum.  The multi-configurational time-dependent Hartree approach.  Chemical Physics Letters, 165:73, 1989.
  * [10] H. Miyagi and L. B. Madsen.  Time-dependent restricted-active-space self-consistent-field theory for laser-driven many-electron dynamics.  Phys. Rev. A, 87(6):062511, jun 2013.
  * [11] M. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer.  The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets.  Phys. Rep., 324(1):1–105, jan 2000.
  * [12] T. Sato and K. L. Ishikawa.  Time-dependent complete-active-space self-consistent-field method for multielectron dynamics in intense laser fields.  Phys. Rev. A, 88(2):023402, aug 2013.
  * [13] D. Hochstuhl, C. Hinz, and M. Bonitz.  Time-dependent multiconfiguration methods for the numerical simulation of photoionization processes of many-electron atoms.  The European Physical Journal Special Topics, 223(2):177–336, feb 2014.
  * [14] R Bulirsch, A. Miele, J. Stoer, and K. H. Well.  Optimal Control: Calculus of Variations, Optimal Control Theory, and Numerical Methods.  International series of numerical mathematics. Springer, 1993.
  * [15] Constantin Brif, Raj Chakrabarti, and Herschel Rabitz.  Control of quantum phenomena: past, present and future.  New Journal of Physics, 12(7):075008, jul 2010.
  * [16] Ignacio R. Sola, Bo Y. Chang, Svetlana A. Malinovskaya, and Vladimir S. Malinovsky.  Quantum Control in Multilevel Systems.  In Advances in Atomic, Molecular and Optical Physics, volume 67, pages 151–256. Academic Press Inc., jan 2018.
  * [17] Alexander S Poznyak.  Chapter 22 - Variational Calculus and Optimal Control.  In Alexander S Poznyak, editor, Advanced Mathematical Tools for Automatic Control Engineers: Deterministic Techniques, pages 647–711. Elsevier, Oxford, 2008.
  * [18] A Prados.  Optimizing the relaxation route with optimal control.  Physical Review Research, 3:023128, 2021.
  * [19] M. Artamonov and T. Seideman.  Optimal-control approach to field-free three-dimensional alignment of polyatomic molecules.  Phys. Rev. A, 82(2):23413, 2010.
  * [20] J Salomon, C M Dion, and G Turinici.  Optimal molecular alignment and orientation through rotational ladder climbing.  J. Chem. Phys., 123:144310, 2005.
  * [21] L. H. Coudert.  Optimal orientation of an asymmetric top molecule with terahertz pulses.  The Journal of Chemical Physics, 146(2):024303, 2017.
  * [22] Ignacio R Sola, Jesus Santamaria, and David J Tannor.  Optimal Control of Multiphoton Excitation: A Black Box or a Flexible Toolkit?  The Journal of Physical Chemistry A, 102(23):4301–4309, jun 1998\.
  * [23] Karsten Sundermann and Regina de Vivie-Riedle.  Extensions to quantum optimal control algorithms and applications to special problems in state selective molecular dynamics.  The Journal of Chemical Physics, 110(4):1896–1904, 1999.
  * [24] Rong Wu, Ignacio R Sola, and Herschel Rabitz.  Optimal quantum control with multi-polarization fields.  Chemical Physics Letters, 400(4):469–475, 2004.
  * [25] Daniel Basilewitsch, Christiane P Koch, and Daniel M Reich.  Quantum Optimal Control for Mixed State Squeezing in Cavity Optomechanics.  Advanced Quantum Technologies, 2:1800110, 2019.
  * [26] Diego S. Acosta Coden, Sergio S. Gómez, Alejandro Ferrón, Omar Osenda, D S Acosta Coden, Sergio S. Gómez, Alejandro Ferrón, and Omar Osenda.  Controlled quantum state transfer in XX spin chains at the Quantum Speed Limit.  Physics Letters A, 387:127009, may 2021.
  * [27] Jesper Hasseriis Mohr Jensen, Frederik Skovbo Møller, Jens Jakob Sørensen, and Jacob Friis Sherson.  Exact Gradients and Hessians for Quantum Optimal Control and Applications in Many-Body Matrix Product States.  may 2020.
  * [28] Jesper Hasseriis Mohr Jensen, Frederik Skovbo Møller, Jens Jakob Sørensen, and Jacob Friis Sherson.  Driving High-Fidelity Many-Body Matrix Product State Dynamics with Exact Gradient Quantum Optimal Control.  pages 1–15, 2020.
  * [29] Henrik R. Larsson and David J. Tannor.  Control of concerted back-to-back double ionization dynamics in helium.  may 2020.
  * [30] Hongling Yu, Tak-San Ho, and Herschel Rabitz.  Optimal control of orientation and entanglement for two dipole–dipole coupled quantum planar rotors.  Physical Chemistry Chemical Physics, 20(18):13008–13029, 2018.
  * [31] Alicia Magann, Linhan Chen, Tak-San Ho, and Herschel Rabitz.  Quantum optimal control of multiple weakly interacting molecular rotors in the time-dependent Hartree approximation.  The Journal of Chemical Physics, 150(16):164303, 2019.
  * [32] N. Dupont, G. Chatelain, L. Gabardos, M. Arnal, J. Billy, B. Peaudecerf, D. Sugny, and D. Guéry-Odelin.  Quantum state control of a Bose-Einstein condensate in an optical lattice.  may 2021.
  * [33] S Amri, R Corgier, D Sugny, E M Rasel, N Gaaloul, and E Charron.  Optimal control of the transport of Bose-Einstein condensates with atom chips.  Scientific Reports, 9(1):5346, 2019.
  * [34] J-F Mennemann, D Matthes, R-M Weishäupl, and T Langen.  Optimal control of Bose–Einstein condensates in three dimensions.  New Journal of Physics, 17(11):113027, nov 2015.
  * [35] Allon Bartana, Ronnie Kosloff, and David J Tannor.  Laser cooling of internal degrees of freedom. II.  The Journal of Chemical Physics, 106(4):1435–1448, 1997.
  * [36] Yukiyoshi Ohtsuki, Wusheng Zhu, and Herschel Rabitz.  Monotonically convergent algorithm for quantum optimal control with dissipation.  The Journal of Chemical Physics, 110(20):9825–9832, 1999.
  * [37] Bernard Bonnard, Monique Chyba, and Dominique Sugny.  Time-minimal control of dissipative two-level quantum systems: The generic case.  IEEE Transactions on Automatic Control, 54(11):2598–2610, nov 2009\.
  * [38] Bernard Bonnard and Dominique Sugny.  Time-Minimal Control of Dissipative Two-Level Quantum Systems: The Integrable Case.  SIAM Journal on Control and Optimization, 48(3):1289–1308, jan 2009\.
  * [39] J. Werschnik and E. K. U. Gross.  Quantum optimal control theory.  Journal of Physics B: Atomic, Molecular and Optical Physics, 40(18):R175, 2007.
  * [40] Gabriel G. Balint-Kurti, Shiyang Zou, and Alex Brown.  Optimal Control Theory for Manipulating Molecular Processes, pages 43–94.  John Wiley & Sons, Inc., 2008.
  * [41] Michael H. Goerz, Sebastián C. Carrasco, and Vladimir S. Malinovsky.  Quantum Optimal Control via Semi-Automatic Differentiation.  Quantum, 6:871, December 2022.
  * [42] V. F. Krotov and I. N. Feldmann.  An iterative method for solving optimal-control problems.  ENGINEERING CYBERNETICS, 21(2):123–130, 1983.
  * [43] David J Tannor, Vladimir Kazakov, and Vladimir Orlov.  Control of Photochemical Branching: Novel Procedures for Finding Optimal Pulses and Global Upper Bounds.  In J Broeckhove and L Lathouwers, editors, Time-Dependent Quantum Molecular Dynamics, pages 347–360. Springer US, Boston, MA, 1992.
  * [44] József Somlói, Vladimir A Kazakov, and David J Tannor.  Controlled dissociation of I2 via optical transitions between the X and B electronic states.  Chemical Physics, 172(1):85–98, 1993.
  * [45] Reuven Eitan, Michael Mundt, and David J. Tannor.  Optimal control with accelerated convergence: Combining the Krotov and quasi-Newton methods.  Phys. Rev. A, 83(5):053426, may 2011.
  * [46] Ido Schaefer and Ronnie Kosloff.  Optimization of high-order harmonic generation by optimal control theory: Ascending a functional landscape in extreme conditions.  Phys. Rev. A, 101(2):023407, feb 2020.
  * [47] Wusheng Zhu and Herschel Rabitz.  A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator.  The Journal of Chemical Physics, 109(2):385–391, 1998.
  * [48] Wusheng Zhu, Jair Botina, and Herschel Rabitz.  Rapidly convergent iteration methods for quantum optimal control of population.  The Journal of Chemical Physics, 108(5):1953–1963, 1998.
  * [49] Yvon Maday and Gabriel Turinici.  New formulations of monotonically convergent quantum control algorithms.  The Journal of Chemical Physics, 118(18):8191–8196, 2003.
  * [50] Yukiyoshi Ohtsuki, Gabriel Turinici, and Herschel Rabitz.  Generalized monotonically convergent algorithms for solving quantum optimal control problems.  Journal of Chemical Physics, 120(12):5509–5517, 2004.
  * [51] Yukiyoshi Ohtsuki, Yoshiaki Teranishi, Peter Saalfrank, Gabriel Turinici, and Herschel Rabitz.  Monotonically convergent algorithms for solving quantum optimal control problems described by an integrodifferential equation of motion.  Physical Review A, 75(3):033407, 2007.
  * [52] Tak-San Ho and Herschel Rabitz.  Accelerated monotonic convergence of optimal control over quantum dynamics.  Physical Review E, 82(2):026703, 2010.
  * [53] Sheng-Lun Liao, Tak-San Ho, Shih-I. Chu, and Herschel Rabitz.  Fast-kick-off monotonically convergent algorithm for searching optimal control fields.  Physical Review A, 84(3):031401, 2011.
  * [54] Daniel M Reich, Mamadou Ndong, and Christiane P Koch.  Monotonically convergent optimization in quantum control using Krotov’s method.  The Journal of Chemical Physics, 136(10):104103, 2012.
  * [55] O. V. Morzhin and A. N. Pechen.  Krotov method for optimal control of closed quantum systems.  Russian Mathematical Surveys, 74(5):851–908, 2019.
  * [56] Michael H. Goerz, Daniel Basilewitsch, Fernando Gago-Encinas, Matthias G. Krauss, Karl P. Horn, Daniel M. Reich, and Christiane P. Koch.  Krotov: A Python implementation of Krotov’s method for quantum optimal control.  SciPost Phys., 7(6):80, feb 2019.
  * [57] Herschel A. Rabitz, Michael M. Hsieh, and Carey M. Rosenthal.  Quantum optimally controlled transition landscapes.  Science, 303(5666):1998–2001, 2004.
  * [58] Katharine W. Moore and Herschel Rabitz.  Exploring quantum control landscapes: Topology, features, and optimization scaling.  Physical Review A, 84(1):012109, 2011.
  * [59] Michael E Peskin and Daniel V Schroeder.  An Introduction to quantum field theory.  Addison-Wesley, Reading, USA, 1995.
  * [60] H.-D. Meyer and G. A. Worth.  Quantum molecular dynamics: propagating wavepackets and density operators using the multiconfiguration time-dependent Hartree method.  Theor. Chem. Acc., 109(5):251–267, 2003.