\section{Quantum optimal control}

In recent years, we have added to Octopus some of the key advancements of quantum
optimal-control theory (QOCT)~\cite{Brif2010,Werschnik2007}. In this section,
we will briefly summarize what this theory is about, overview the current
status of its implementation, and describe some of the results that have been
obtained with it until now.

\textit{Quantum control} can be loosely defined as the manipulation of physical
processes at the quantum level. We are concerned here with the theoretical
branch of this discipline, whose most general formulation is precisely
QOCT. This is, in fact, a particular case of the general mathematical field of
``optimal control'', which studies the optimization of dynamical
processes in general. The first applications of optimal control in the quantum
realm appeared in the 80s~\cite{Shi1988,Peirce1988,Kosloff1989a}, and the
field has rapidly evolved since then. Broadly speaking, QOCT attempts to
answer the following question: given a quantum process governed by a
Hamiltonian that depends on a set of parameters, what are the values of those
parameters that maximize a given observable that depends on the behavior of
the system? In mathematical terms: let a set of parameters $u_1,\dots,u_M \equiv u$ determine the Hamiltonian of a system $\hat{H}[u,t]$, so that the
evolution of the system also depends on the value taken by those parameters:
\begin{align}
{\rm i}\frac{{\rm d}}{{\rm d}t}\vert\psi(t)\rangle & = 
\hat{H}[u,t]\vert\psi(t)\rangle\,,
\\
\vert\psi(0)\rangle & = \vert\psi_0\rangle\,,
\end{align}
\textit{i.e.}\ the solution of the Schr{\"{o}}dinger equation determines a map $u
\longrightarrow \psi[u]$. Suppose we wish to optimize a
functional of the system $F=F[\psi]$. QOCT is about finding the extrema of
$G(u)=F[\psi[u]]$. Beyond this search, QOCT also studies topics such
as the robustness of the optimal solutions for those parameters, the number of
solutions, or the construction of suitable algorithms to compute them.

Perhaps the most relevant result of QOCT is the equation for the gradient of $G$, which allows use of the various maximization algorithms available. For the simple formulation given above, this gradient is given by
\begin{equation}
\frac{\partial G}{\partial u_m}(u) = 2\ {\rm Im} \int_0^T\!\!{\rm d}t\;
\left< \chi(t)\left|\frac{\partial \hat{H}}{\partial u_m}[u,t]\right|
\psi(t) \right>\,,
\label{eq:qoctgradient}
\end{equation}
where $\chi$ is the \emph{costate}, an auxiliary wave function that is defined through the following equation of motion:
\begin{align}
{\rm i}\frac{\rm d}{{\rm d}t}\vert\chi(t)\rangle & =  \hat{H}^\dagger[u,t]
\vert\chi(t)\rangle\ ,
\\
\vert \chi(T) \rangle & =  \frac{\delta F}{\delta \psi^*(T)}\ .
\end{align}
This equation assumes, in order to keep this description short, that the target functional $F$ depends on the state of the system only at the final time of the propagation $T$, \textit{i.e.}\ it is a functional of $\psi(T)$. Note the presence of a boundary value equation at the final time of the propagation, as opposed to the equation of motion for the ``real'' system $\psi$, which naturally depends on an initial value condition at time zero. With these simple equations, we may already summarize what is needed from an implementation point of view in order to perform basic QOCT calculations:

The first step is the selection of the parameters $u$, that constitute the \emph{search space}. Frequently, these parameters are simply the values that the \emph{control function} (typically, the electric-field amplitude) takes at the time intervals that are used to discretize the propagation interval, \textit{i.e.} it is a ``real-time parametrization''. However, more sophisticated parametrizations allow fine-tuning of the search space, introducing constraints and penalties into the formulation.

Then, one must choose an algorithm for maximizing multi-dimensional
functions such as $G$. One possibility is the family of \emph{gradient-less} algorithms, which only require a procedure to compute the value of the function, and do not need the gradient. In this case, the previous equations are obviously not needed. One only has to propagate the system forwards in time, which is what Octopus can do best. The value of the function $G$ can then be computed from the evolution of $\psi$ obtained with this propagation, and fed into the optimization procedure. A few gradient-less algorithms are implemented in Octopus.

The most efficient optimizations can be obtained if information about the gradient is employed. In that case, we can use standard schemes, such as the family of conjugate-gradient algorithms, or the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton scheme -- we use the implementation of these algorithms included in the GSL mathematical library~\cite{Galassi2009}. Some \textit{ad hoc} algorithms, developed explicitly for QOCT, exist. These may in some circumstances be faster than the general purpose ones. Some of those are implemented in Octopus as well~\cite{Zhu1998,Zhu1998a,Ohtsuki1999}.

In order to compute the gradient, one must implement a \emph{backwards-propagation} scheme for the costate, which does not differ from the ones used for the normal forwards propagation~\cite{Castro_2004}. Note, however, that in some cases the backwards propagation does not have the exact same simple linear form than the forwards propagation, and may include inhomogeneous or non-linear terms. The final step is the computation of the gradient from the integral given in eq.~(\ref{eq:qoctgradient}).

The formulation of QOCT we have just sketched out is quite generic; in our case the quantum systems are
those that can be modeled with Octopus (periodic systems are not supported at the moment), and the
\textit{handle} that is used to control the system is a time-dependent electric field, such as the
ones that can be used to model a laser pulse. The set of parameters \(
\left\{u\right\}_i\) define the shape of this electric field; for example, they can be the Fourier coefficients of the field amplitude.

The usual formulation of QOCT assumes the linearity of quantum mechanics. However, the time-dependent KS equations are not linear, making both the theory and the numerics more complicated. We have extended the basic theory previously described to handle the TDDFT equations, and implemented the resulting equations in Octopus~\cite{Castro2012a}.

We conclude this section by briefly describing some of the applications of the QOCT machinery included in Octopus, which can give an idea of the range of possibilities that can be attempted.  The study presented in Ref.~\cite{Rasanen2007} demonstrates the control of single-electron states in a two-dimensional semiconductor quantum-ring model. The states whose transitions are manipulated are the current-carrying states, which can be populated or de-populated with the help of circularly polarized light. 

Reference~\cite{Rasanen2008} studies double quantum dots, and shows how the electron state of these systems can be manipulated with the help of electric fields tailored by QOCT.

Another interesting application is how to tailor the shape of femtosecond laser pulses in order to obtain maximal ionization of atoms and molecules~\cite{Castro2009e}. The system chosen to demonstrate this possibility is the H$_2^+$ molecule, irradiated with short ($\approx 5$~fs) high-intensity laser pulses.

The feasibility of using the electronic current to define the target functional of the QOCT formalism is considered in Ref.~\cite{Kammerlander2011a}. 

Finally, a series of works has studied the use of optimal control for photo-chemical control: the tailoring of laser pulses to create or break selected bonds in molecules. The underlying physical model should be based on TDDFT, and on a mixed quantum/classical scheme (within Octopus, Ehrenfest molecular dynamics). Some first attempts in this area were reported in Refs~\cite{Krieger2011,Castro14}. However, these works did not consider a fully consistent optimal control theory encompassing TDDFT and Ehrenfest dynamics. This theory has been recently presented~\cite{Castro2014}, and the first computations demonstrating its feasibility will be reported soon.
