# Effects of geomagnetic field perturbations on the power supply of
transoceanic fiber optic cables

Antonio Mecozzi A. Mecozzi is with the Department of Physical and Chemical
Sciences, University of L’Aquila, L’Aquila 67100, Italy

###### Abstract

There is a growing concern that a big coronal mass ejection event will induce
perturbations on the power supply of fiber optic transoceanic cables that may
produce a global internet blackout. In this paper we give the expression of
the voltage variations that a transient change of the geomagnetic field
induces on the voltage of the power supply of a transoceanic fiber optic
cable. We show that the transient voltage change is proportional to the
magnitude of the magnetic field deviations and not to its time derivative as a
direct application of Faraday’s law would imply, and this suggests design
criteria to protect transoceanic fiber optic systems against big geomagnetic
storm events. The presented analysis also enables the classification of
existing systems into some that are less sensitive to the weakening of the
geomagnetic field occurring during strong geomagnetic storms and others that
are more prone to experience an outage when a weakening of the geomagnetic
field occurs.

###### Index Terms:

Optical communications, Transoceanic optical cables

##  I Introduction

The interaction of the geomagnetic field with charged particles during coronal
mass ejection (CME) may induce significant electromagnetic perturbations,
known as geomagnetic storms. There is a growing concern that a big CME event
will induce perturbations on fiber optic transoceanic systems that may end up
in a global internet blackout [1].

In this paper, we present a theory of the interaction between the power line
feeding a transoceanic fiber optic system and a time-varying magnetic field.
We show that the combined effect of the frequency dependent attenuation of the
electromagnetic perturbation in the conducting seawater and the reflection
from the bottom, has the effect that the induced electric field in not
proportional to the derivative of the magnetic field as a direct application
of Faraday’s law would suggest, but to the amplitude of the magnetic field
perturbation. The theory predicts that the effective area that was used to
quantify the effect of the magnetic perturbation on the induced voltage
variations becomes linearly proportional to the duration of the perturbation,
and this result is used to check the validity of the analysis with
experimental data extracted from the literature [2]. The predicted
proportionality of the voltage variations to the amplitude of the magnetic
field perturbations has also very practical implications. In the presence of a
sudden drop of the magnetic field, like the ones that are likely to occur
during large magnetic storms [3, 4], the induced voltage variations will
follow the magnetic field perturbation, and undergo an increase or a decrease
depending on the direction of the electric current. This result, beside
providing a simple criterion to estimate the vulnerability of an existing
fiber optic plan, provides design criteria to reduce the sensitivity of newly
deployed systems to possible future strong geomagnetic storm events.

##  II The model

In a transoceanic optical cable, the in-line equipment is fed by a positive
and a negative (or zero) DC voltage from the cable terminals using a wire that
is grounded at its ends. Effectively, the return line is replaced by the
earth. The voltage balance is

$V_{1}+V_{+}-V_{\mathrm{sys}}-V_{-}=V_{2},$ (1)

where $V_{1}$ is the potential of the earth at the transmitter (or the
receiver) side, $V_{+}$ is the positive voltage of the feed at the transmitter
(or the receiver) side, $V_{\mathrm{sys}}$ is the voltage drop due to the line
resistance and the voltage drop of the in-line devices, $V_{-}$ is the
absolute value of the negative voltage of the feed at the receiver (or the
transmitter) side, and $V_{2}$ is the potential of the earth at the receiver
(or the transmitter) side. Writing

$V_{+}-V_{-}=V_{\mathrm{sys}}+\Delta V,$ (2)

where $\Delta V=V_{2}-V_{1}$, one sees right away that the required feeding
voltage is equal to the voltage drop of the system plus the difference of the
earth potential at the receiver and the transmitter side (or vice versa). In a
static situation, the earth voltage is uniform. This is not the case when the
geomagnetic field that surrounds the earth changes with time. Under the action
of a time varying magnetic fields, electric fields are generated by Faraday’s
law. For slowly varying fields like the ones involved in geomagnetic storm
events, the quasi static approximation is always valid, so that an electric
potential can still be defined. In this case, however, the presence of a time-
varying magnetic field produces a non-uniform earth potential. Therefore, if
one assumes that the voltage drop of the system is constant, the voltage
variations of the power supply of transoceanic transmission systems monitor
the variations of the difference of the earth potentials at the transmitter
and receiver sides.

In order to analyze this problem quantitatively, let us write Faraday, Ampère-
Maxwell and Ohm laws in the ocean, in the quasi static approximation, as
follows

$\displaystyle\vec{\nabla}\times\vec{E}$ $\displaystyle=$
$\displaystyle-\frac{\partial\vec{B}}{\partial t},$ (3)
$\displaystyle\vec{\nabla}\times\vec{B}$ $\displaystyle=$
$\displaystyle\mu_{0}\vec{J}+\mu_{0}\epsilon\frac{\partial\vec{E}}{\partial
t},$ (4) $\displaystyle\vec{J}$ $\displaystyle=$
$\displaystyle\sigma\left(\vec{E}+\vec{v}\times\vec{B}\right).$ (5)

Here, $\sigma$ is the electrical conductivity of the water, and we neglect the
very small perturbation introduced by the cable wire. Let us consider the case
$v=0$ first, and a time varying magnetic field, and assume the quasi static
approximation neglecting the Maxwell term proportional to the time derivative
of the electric field in the Ampère-Maxwell law (4). Using then Eq. (5) into
Eq. (4) we obtain

$\vec{\nabla}\times\vec{B}=\mu_{0}\sigma\vec{E}.$ (6)

Inserting the traveling wave solution with $\omega>0$ and
$\vec{B}(-\omega)=\vec{B}^{*}(\omega)$ and
$\vec{E}(-\omega)=\vec{E}^{*}(\omega)$

$\displaystyle\vec{B}$ $\displaystyle=$
$\displaystyle\vec{B}_{0}\exp\left(i\vec{k}\cdot\vec{r}-i\omega t\right),$ (7)
$\displaystyle\vec{E}$ $\displaystyle=$
$\displaystyle\vec{E}_{0}\exp\left(i\vec{k}\cdot\vec{r}-i\omega t\right),$ (8)

we obtain

$\displaystyle\vec{k}\times\vec{E}_{0}$ $\displaystyle=$
$\displaystyle\omega\vec{B}_{0},$ (9) $\displaystyle\vec{k}\times\vec{B}_{0}$
$\displaystyle=$ $\displaystyle-i\mu_{0}\sigma\vec{E}_{0}.$ (10)

The vectors $\vec{E}_{0}$, $\vec{B}_{0}$ and $\vec{k}$ are a right-handed
orthogonal set of vectors. We have

$\vec{k}\times\left(\vec{k}\times\vec{E}_{0}\right)=-i\mu_{0}\sigma\omega\vec{E}_{0},$
(11)

that is

$-k^{2}\vec{E}_{0}=-i\mu_{0}\sigma\omega\vec{E}_{0}\Longrightarrow
k^{2}=i\mu_{0}\sigma\omega.$ (12)

Figure 1: Schematic of the model.

Let us now approximate the ocean as a flat water layer of width $H$ and
resistivity $\sigma$, with air on top and lying on a flat semi-infinite
substrate (the mantle) of resistivity $\sigma_{1}$ [6], and use a reference
frame defined by the base unit vectors $\hat{e}_{z}$ directed downwards,
$\hat{e}_{x}$ parallel to $\vec{E}_{0}$, and $\hat{e}_{y}$ parallel to
$\vec{B}_{0}$ to complete a right-handed reference frame, see Fig. 1. Assume
that the cable is lying at the interface of water and the flat semi-infinite
substrate, at depth $z=H$. Assume that the electromagnetic field is
propagating in the water along the $z$ direction so that $\vec{k}=\pm
k\hat{e}_{z}$. The above equation yields two attenuated traveling wave
solutions proportional to

$\exp\left(ikz-i\omega t\right)=\exp(\pm\kappa\,z-i\omega t).$ (13)

where

$\kappa=\sqrt{-i\mu_{0}\sigma\omega},$ (14)

and we use the convention that the real part of $\kappa$ is positive. With
this choice if we take the minus sign we have the attenuated progressive wave
propagating in the direction $\hat{e}_{z}$

$\vec{E}_{+}=\vec{E}_{0}\exp(-\kappa\,z-i\omega
t)=\vec{E}_{0}\exp\left[(i-1)\frac{z}{\delta(\omega)}-i\omega t\right]$ (15)

and if we take the positive sign the attenuated regressive wave propagating in
the direction opposite to $\hat{e}_{z}$

$\vec{E}_{-}=\vec{E}_{0}\exp(\kappa\,z-i\omega
t)=\vec{E}_{0}\exp\left[-(i-1)\frac{z}{\delta(\omega)}-i\omega t\right],$ (16)

where

$\delta(\omega)=\sqrt{\frac{2}{\mu_{0}\sigma\omega}}$ (17)

is the skin attenuation depth. Let us define the Fourier amplitudes as
$E(z,t)=E(z)\exp(-i\omega t)$ and $B(z,t)=B(z)\exp(-i\omega t)$, assume

$\displaystyle\vec{E}(z)$ $\displaystyle=$ $\displaystyle E(z)\hat{e}_{x},$
(18) $\displaystyle\vec{B}(z)$ $\displaystyle=$ $\displaystyle
B(z)\hat{e}_{y},$ (19)

and decompose that the fields as a combination of forward and backward
propagating components, we have

$\displaystyle E(z)$ $\displaystyle=$ $\displaystyle
E_{+}\exp(-\kappa\,z)+E_{-}\exp(\kappa\,z),$ (20) $\displaystyle B(z)$
$\displaystyle=$ $\displaystyle B_{+}\exp(-\kappa\,z)-B_{-}\exp(\kappa\,z).$
(21)

Using Eqs. (3) and (4) in the quasi static approximation we obtain

$\frac{E_{+}}{B_{+}}=\frac{E_{-}}{B_{-}}=\frac{ik}{\mu_{0}\sigma}=\sqrt{\frac{-k^{2}}{\mu_{0}^{2}\sigma^{2}}}=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}=\frac{Z}{\mu_{0}},$
(22)

where $Z=E_{+}/H_{+}=\sqrt{-i\omega\mu_{0}/\sigma}$ is the impedance of the
ocean water, so that Eq. (21) becomes

$B(z)=\sqrt{\frac{\mu_{0}\sigma}{-i\omega}}\left[E_{+}\exp(-\kappa\,z)-E_{-}\exp(\kappa\,z)\right],$
(23)

Defining the reflection coefficient $\Gamma$ of an incident wave on the
substrate as

$\Gamma=\frac{E_{-}\exp(\kappa\,H)}{E_{+}\exp(-\kappa\,H)},$ (24)

and using the preservation of the parallel components of $\vec{E}$ and
$\vec{B}$ at the interface and the absence of a backward propagating field
returning from the semi-infinite medium, it is possible to show that $\Gamma$
has the expression

$\Gamma=\frac{Z_{1}-Z}{Z_{1}+Z}=\frac{(\sigma_{1})^{-1/2}-(\sigma)^{-1/2}}{(\sigma_{1})^{-1/2}+(\sigma)^{-1/2}},$
(25)

where we used that the impedance of the mantle is
$Z_{1}=\sqrt{-i\omega\mu_{0}/\sigma_{1}}$. The use in Eq. (20) and in Eq. (23)
calculated for $z=0$ of the expression for $E_{-}$ obtained from Eq. (24),
that is $E_{-}=\Gamma\exp(-2\kappa\,H)E_{+}$, yields

$\displaystyle\frac{E(z)}{E_{+}}$ $\displaystyle=$
$\displaystyle\exp(-\kappa\,z)+\Gamma\exp[\kappa\,(z-2H)],$ (26)
$\displaystyle\frac{B(0)}{E_{+}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\mu_{0}\sigma}{-i\omega}}\left[1-\Gamma\exp(-2\kappa\,H)\right].$
(27)

The ratio of the first to the second of the equations above gives

$\frac{E(z)}{B(0)}=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\,\frac{\exp[\kappa\,(H-z)]+\Gamma\exp[-\kappa\,(H-z)]}{\exp(\kappa\,H)-\Gamma\exp(-\kappa\,H)},$
(28)

or, using Eq. (25) and after straightforward algebra111This equation
generalizes the combination of Eqs. (27) and (28) of ref. [6] that gives
$Z=\mu_{0}E(0)/B(0)$, and becomes identical to that for $z=0$, after the
replacements $\sigma\to\sigma_{1}$, $\sigma_{1}\to\sigma_{3}$ and $H\to
h_{1}$, and with $-i$ replaced by $i$ because of the different Fourier
transform definition.

$\frac{E(z)}{B(0)}\\!=\\!\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\frac{\cos[\kappa\,(H-z)]+\left(\sigma_{1}/\sigma\right)^{1/2}\sinh[\kappa\,(H-z)]}{\sinh(\kappa\,H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)}.$
(29)

For $z=H$, Eq. (29) becomes

$E(H)=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\,\frac{B(0)}{\sinh(\kappa\,H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)},$
(30)

where the magnetic field $B(0)$ is the total magnetic field evaluated at the
interface between water and air. Using now Eq. (18) for $z=H$ and Eq. (19) for
$z=0$, and that $\hat{e}_{x}=\hat{e}_{y}\times\hat{e}_{z}$, we obtain the
final result [7]

$\vec{E}(H)=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\,\frac{\vec{B}(0)\times\hat{e}_{z}}{\sinh(\kappa\,H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)}.$
(31)

Outside the water, the propagation constant is $k=\omega/c$ so that the
wavelength is $\lambda=cT$ where $T=2\pi/\omega$ is the period of the
perturbation. For $T$ larger than one second, the wavelength is larger than
$300,000$ km, larger than the size of the magnetosphere. This means that the
plane wave approximation cannot be applied in air and the quantity of interest
is the total magnetic field at the air-water interface, not a forward-
propagating magnetic field in air, which cannot be defined outside the
framework of a plane wave approach. For the generation of a planar wave in the
ocean layer, instead, we need to assume that the amplitude and phase of the
geomagnetic field is uniform on the water surface over the length scale set by
the wavelength in the water $2\pi/\mathrm{imag}(k)=2\pi\delta(\omega)$. For
time scales $\Delta\tau=2\pi/\omega$ ranging from $10^{3}$ to $10^{5}$ s,
which are typical time scales of geomagnetic field fluctuations,
$\delta(\omega)$ lies between $9$ to $90$ kilometers for
$\sigma=3\,\Omega^{-1}$m-1.

Let us now consider a cable parametrized by the curvilinear coordinate $s$ and
be ${\operatorname{d}}\vec{s}$ its increment. We will allow that the direction
of the magnetic field varies along the cable over a length scale long enough
that the plane wave approximation still holds in the water layer, and we also
allow a variability of $H$ over the same length scale. Let us assume that the
return path of the circuit goes deep into the earth mantle in a region where
$\vec{E}\simeq 0$. The induced electromotive force in the circuit acts only in
the cable section and its component at frequency $\omega$ is

${\cal
E}=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\int_{0}^{L}\frac{\hat{e}_{z}\times{\operatorname{d}}\vec{s}\cdot\vec{B}(0)}{\sinh(\kappa\,H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)},$
(32)

where we used that
$\vec{B}(0)\times\hat{e}_{z}\cdot{\operatorname{d}}\vec{s}=\hat{e}_{z}\times{\operatorname{d}}\vec{s}\cdot\vec{B}(0)$.
Let us now use a local reference frame defined by the base unit vectors
$\hat{e}_{z}$, $\hat{e}_{x}$ and $\hat{e}_{y}$ where $\hat{e}_{x}$ is defined
by ${\operatorname{d}}\vec{s}=\hat{e}_{x}{\operatorname{d}}s$ and
$\hat{e}_{y}$ by the property that it completes a right-handed orthogonal
reference frame together with the vertical $\hat{e}_{z}$ pointing downwards,
and assume that the cable is lying at depth $z=H(s)$. In this reference frame
Eq. (32) becomes

${\cal
E}=\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\int_{0}^{L}\frac{B_{y}(\omega){\operatorname{d}}s}{\sinh(\kappa\,H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)},$
(33)

where $B_{y}(\omega)=\hat{e}_{y}\cdot\vec{B}(0)$, and we made explicit the
dependence of the magnetic field on the angular frequency $\omega$ and
implicit the fact that the magnetic field is referred to the interface ocean-
air, that is at $z=0$. A constant intensity in the cable requires that the
total voltage $V$ needed for the nominal current flow in the system be
constant in the presence of the perturbation ${\cal E}$, that is $V+\Delta
V+{\cal E}=V$, so that the voltage change $\Delta V$ in response to the
generation of the electromotive force $\cal E$ is $\Delta V=-{\cal E}$, that
is

$\Delta
V=-\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\int_{0}^{L}\frac{B_{y}(\omega){\operatorname{d}}s}{\sinh(\kappa
H)+\left(\sigma_{1}/\sigma\right)^{1/2}\cosh(\kappa\,H)}.$ (34)

When the current flows in the circuit in the positive direction of $s$ and of
$\hat{e}_{x}$, $V$ is positive and hence a positive $\Delta V$ corresponds to
an increase of the voltage $V$ of the power supply. On the other hand, when
the current flows in the opposite direction, $V$ is negative and hence a
positive $\Delta V$ corresponds to a decrease of the absolute value of $V$.

##  III Approximate analytical expressions

Equation (34) shows that $\Delta V(\omega)=0$ for $\omega=0$ and hence $\Delta
V(t)=0$ for a constant $B_{y}(t)$. We are interested into events which have a
finite time duration over a constant background. We may therefore separate the
constant background using

$B_{y}(t)=B_{y,0}+\Delta B_{y}(t),$ (35)

where $\Delta B_{y}(t)\to 0$ for $|t|\to\infty$ sufficiently fast that its
spectrum is finite. Because $\Delta V(\omega)$ does not change if $B_{y}(t)$
is replaced by $\Delta B_{y}(t)$ by subtracting a constant bias, it is
legitimate to replace in Eq. (34) $B_{y}$ by $\Delta B_{y}$ in the following
analysis.

In order to obtain some useful analytical relation between the voltage and the
magnetic field variations, let us expand the denominator of Eq. (34) the
arguments of the hyperbolic functions as

$\Delta
V(\omega)=-\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\int_{0}^{L}\frac{\Delta
B_{y}(\omega){\operatorname{d}}s}{\alpha+\alpha^{3}/6+\left(\sigma_{1}/\sigma\right)^{1/2}(1+\alpha^{2}/2)},$
(36)

where $\alpha=\kappa\,H$. If we assume that $\Delta B_{y}(\omega)$ is
bandwidth limited to the frequency $\Omega=2\pi/\Delta\tau$, then we have
$|\alpha|<1$ for $\Delta\tau>2\pi\mu_{0}\sigma H^{2}$. Assuming
$\sigma=3\,\Omega^{-1}$ m-1 for the value of the conductivity of the ocean
water, and $H=4$ km, we obtain $\Delta\tau>380$ s. For events that fulfill
this condition we may then write

$\Delta
V(\omega)=-\sqrt{\frac{-i\omega}{\mu_{0}\sigma}}\int_{0}^{L}\frac{\Delta
B_{y}(\omega){\operatorname{d}}s}{\alpha+\left(\sigma_{1}/\sigma\right)^{1/2}}.$
(37)

The value of the conductivity of the substrate is much smaller than that of
the water, of the order of $0.01\,\Omega^{-1}$ m-1 or less. Assuming the value
$\sigma_{1}=0.01\,\Omega^{-1}$ m-1 the second term at the denominator is
dominant for $\mu_{0}\sigma(2\pi/\Delta\tau)H^{2}>(\sigma_{1}/\sigma)$, that
is for $\Delta\tau<(\sigma/\sigma_{1})\mu_{0}\sigma 2\pi H^{2}$ which, with
our values of the parameters is $\Delta\tau<10^{5}$ s. For $\Delta\tau>380$ s
and $\Delta\tau<10^{5}$ both the first order expansion and the neglecting of 1
are legitimate and we have

$\Delta V(\omega)=-\int_{0}^{L}\frac{\Delta
B_{y}(\omega){\operatorname{d}}s}{\mu_{0}\sigma H(s)},$ (38)

and in time domain

$\Delta V(t)=-\frac{L}{\mu_{0}\sigma}\left[\frac{1}{L}\int_{0}^{L}\frac{\Delta
B_{y}(t)\,{\operatorname{d}}s}{H(s)}\right].$ (39)

For constant $H(s)=H$ and constant $\Delta B_{y}(t)$ we then obtain the simple
relation

$\Delta V(t)=-\frac{L}{\mu_{0}\sigma H}\,\Delta B_{y}(t),$ (40)

Assuming $H=4$ km and $\sigma=3\,\Omega^{-1}$m-1, we obtain $\Delta
V/(L\,\Delta B_{y})\simeq 66.3$ mV/(Mm nT) (Mm stands for megameter, a handy
unit for transoceanic system lengths). Notice that, contrary to what could be
anticipated from a direct application of the Faraday’s law, the voltage
variations are proportional to the amplitude of the magnetic field
perturbations, not to its time derivative. This can be explained by the fact
that the derivative becomes multiplication by $-i\omega$ in frequency domain,
and that the coupling between the induced electric field and the derivative of
the magnetic field is inversely proportional to $\sqrt{-i\omega\mu_{0}\sigma}$
and the enhancement originated by the reflection of the substrate is in turn,
for high enough frequencies, inversely proportional to
$\sqrt{-i\omega\mu_{0}\sigma}\,H$, so that the effective coupling between the
voltage and the magnetic field perturbations turns out to be independent of
$\omega$.

If we use in Eq. (40) the dimensionless time $\tau=t/\Delta\tau$, we obtain

$\Delta V(\tau)=-\frac{L\,\Delta B_{y}(\tau)}{\mu_{0}\sigma H},$ (41)

which shows that a consequence of the dependence of the voltage variations on
the amplitude of the magnetic field perturbations is that the proportionality
coefficient is independent on the time scale $\Delta\tau$ of the magnetic
field perturbation.

For very low frequencies, corresponding to $\Delta\tau>10^{5}$ s, we are in
the opposite limit in which we may instead neglect the term proportional to
$H$ in Eq. (37). If we do so, we obtain

$\Delta V(\omega)=-\sqrt{\frac{-i\omega}{\mu_{0}\sigma_{1}}}\int_{0}^{L}\Delta
B_{y}(\omega){\operatorname{d}}s.$ (42)

An elegant and very useful formal expression in time domain of the above
equation can be obtained in terms of the derivative of order one-half of the
magnetic field (the derivative of order one-half is a linear operator that
applied twice produces a first derivative [8])

$\Delta
V(t)=-\sqrt{\frac{1}{\mu_{0}\sigma_{1}}}\int_{0}^{L}\left(\frac{{\operatorname{d}}}{{\operatorname{d}}t}\right)^{1/2}\,\Delta
B_{y}(t){\operatorname{d}}s,$ (43)

In this case, the use of the dimensionless time $\tau=t/\Delta\tau$ gives

$\Delta
V(\tau)=-\sqrt{\frac{1}{\mu_{0}\sigma_{1}\Delta\tau}}\int_{0}^{L}\left(\frac{{\operatorname{d}}}{{\operatorname{d}}\tau}\right)^{1/2}\,\Delta
B_{y}(\tau){\operatorname{d}}s.$ (44)

For a gaussian $\Delta B_{y}(\tau)$ the ratio between the maxima of the half-
derivative and the derivative is about 3, and it is in all cases of the order
of unit. Equation (44) shows an inverse proportionality of $\Delta V(\tau)$
with the square root of the time scale of the perturbation $\Delta\tau$.

Direct application of Faraday’s law suggests the definition of the effective
area $A$ [5, 2]

$\Delta V(t)=-A\frac{{\operatorname{d}}\Delta B_{y}(t)}{{\operatorname{d}}t},$
(45)

that is, using the dimensionless time $\tau=t/\Delta\tau$,

$A=-\Delta\tau\Delta V(\tau)\left[\frac{{\operatorname{d}}\Delta
B_{y}(t)}{{\operatorname{d}}t}\right]^{-1}.$ (46)

Using for $V(\tau)$ the expression given by Eq. (41), valid for
$\Delta\tau>380$ s and $\Delta\tau<10^{5}$ s, and expanding $A$ to lowest
order in $\Delta\tau$, we obtain that $A$ is linearly dependent on the time
scale of the magnetic field perturbation $\Delta\tau$

$A=A_{0}\Delta\tau,$ (47)

with proportionality coefficient

$\displaystyle A_{0}$ $\displaystyle=$
$\displaystyle\int_{0}^{L}\frac{{\operatorname{d}}s}{\mu_{0}\sigma
H(s)}\left[\left(\frac{{\operatorname{d}}}{{\operatorname{d}}\tau}\right)\,\ln\Delta
B_{y}(\tau)\right]^{-1}$ (48) $\displaystyle=$
$\displaystyle\frac{L}{\mu_{0}\sigma
H}\left[\left(\frac{{\operatorname{d}}}{{\operatorname{d}}\tau}\right)\,\ln\Delta
B_{y}(\tau)\right]^{-1},$

where the second equality holds for constant $H(s)=H$ and $\Delta
B_{y}(\tau)$. Using in Eq. (46) the expression for $V(\tau)$ given by Eq.
(44), valid for $\Delta\tau>10^{5}$, we obtain that $A$ is proportional to the
square root $\sqrt{\Delta\tau}$ of the time scale of the magnetic field
perturbation

$A=A_{1}\sqrt{\Delta\tau},$ (49)

with proportionality coefficient

$A_{1}=L\sqrt{\frac{1}{\mu_{0}\sigma_{1}}}\left[\frac{\left({\operatorname{d}}/{\operatorname{d}}\tau\right)^{1/2}\,\Delta
B_{y}(\tau)}{\left({\operatorname{d}}/{\operatorname{d}}\tau\right)\,\Delta
B_{y}(\tau)}\right].$ (50)

##  IV Comparison with literature data

In Fig. 2 we plot the values taken from Table 1 of ref. [2] of $A$ in square
meters vs. the time change of the magnetic field $\Delta\tau$ in seconds for
geomagnetic storms events affecting the TAT-8 system (blue squares), the TAT-7
system (red dots), and the TAT-6 system (magenta triangle). A rigorous
comparison between the $A$ obtained with our theory and that defined in Refs.
[5] and [2] is of course impossible because of the qualitative definition of
$\Delta\tau$ and the presence of the derivative of $\ln(B_{y})$ with respect
to $\tau=t/\Delta\tau$ in Eq. (48). However, it is reasonable to assume that
the absolute value of the derivative of $\ln(B_{y})$ is of the order of one
(for instance, the derivative of the $\ln\mathrm{sech}(\tau)$ is
$-\tanh(\tau)$ which spans the interval from $-1$ to 1), so that an estimate
of the scaling of $A$ with the event duration $\Delta\tau$ may be obtained by
replacing in Eq. (48) the derivative of the logarithm of $B_{y}$ with $1$.
Using then $\sigma=6$ $\Omega^{-1}$m-1 for the electrical conductivity of the
seawater, $L=6300$ km and $H=4$ km we obtain $A_{0}=2.1\,10^{9}$ m2s-1. The
plot of $A=A_{0}\Delta\tau$ is reported as a dashed blue line in the figure,
showing good agreement with the experimental points, with the exception of two
outliers for small and large values of $\Delta\tau$.

Figure 2: Values taken from Table 1 of ref. [2] of $A$ in square meters vs.
the time change of the magnetic field $\Delta\tau$ is seconds for geomagnetic
storms events affecting the TAT-8 system (blue squares), the TAT-7 system (red
dots), and the TAT-6 system (magenta triangle). Dashed blue line, plot of the
scaling $A=A_{0}\,\Delta\tau$ (see text).

Let us consider a cable whose path is a straight line east-west, and use a
three dimensional right-handed frame with origin on the ocean surface and base
unit vectors $\hat{e}_{y}$ horizontal pointing north, $\hat{e}_{z}$ vertical
oriented downwards, and $\hat{e}_{x}$ horizontal parallel to the cable and
pointing west. In this reference frame, a positive $\Delta B_{y}(t)$ produces
a negative $\Delta V$. In both TAT-6 and TAT-8 systems the current flows west-
east, and hence the voltage $V$ is negative and a negative $\Delta V$
corresponds to an increase of the absolute value of the voltage. Equation (40)
predicts for a positive $\Delta B_{y}(t)$ a negative $\Delta V(t)$, which
corresponds in our case to an increase of the voltage $V$, and this is
consistent with the experimental observations that show that an increase of
the south-north magnetic field component always corresponds to an increase of
the absolute value of the voltage [5, 2].

Let us now perform a more detailed comparison between the prediction of the
theory and the experimental observations reported in Refs. [5] and [2].

The left panel of Fig. 3 shows the south-north magnetic field horizontal
component averaged over five days, digitized from Fig. 2 of ref. [5]. The
right panel of the same figure shows by a solid blue line line the voltage
variation in the TAT-6 system, of length $L=6300$ km, obtained from Fig. 3 of
the same reference, and by a dashed red line the prediction of Eq. (34) with
$\sigma=6\,\Omega^{-1}$m-1 and $\sigma_{1}=0.01\,\Omega^{-1}$m-1, and by a
dotted green line the prediction of its approximated version Eq. (40) with
$\sigma=7\,\Omega^{-1}$m-1.

Figure 3: Left panel: South-north magnetic field horizontal component averaged
over five days, digitized from Fig. 2 of ref. [5]. Right panel: the voltage
variation in the TAT-6 system obtained from Fig. 3 of the same reference
(solid blue line), and the prediction of Eq. (34) with
$\sigma=6\,\Omega^{-1}$m-1 and $\sigma_{1}=0.01\,\Omega^{-1}$m-1 (dashed red
line), and of Eq. (40) with $\sigma=7\,\Omega^{-1}$m-1 (dotted green line).

Figure 4 compares the results of the theory with the observations recorded
during the event reported in Fig. 2 of ref. [2]. The right panel shows good
agreement of the voltage variation observed experimentally (solid blue line)
with the predictions of Eq. (34) using $\sigma=1.5\,\Omega^{-1}$m-1 and
$\sigma_{1}=0.001\,\Omega^{-1}$m-1 (dashed red line) and with the predictions
of Eq. (40) using $\sigma=1.7\,\Omega^{-1}$m-1 (dotted green line).

Figure 4: Left panel: South-north magnetic field horizontal component
digitized from Fig. 2 of ref. [2]. Right panel: the voltage variation in the
TAT-8 system also obtained from Fig. 2 of ref. [2] (solid blue line), and the
prediction of Eq. (34) with $\sigma=1.5\,\Omega^{-1}$m-1 and
$\sigma_{1}=0.001\,\Omega^{-1}$m-1 (dashed red line) and of Eq. (40) with
$\sigma=1.7\,\Omega^{-1}$m-1 (dotted green line)

Figure 5 compares the results of the theory with the event reported in Fig. 4
of ref. [2]. In this case, the right panel shows good agreement of the voltage
variation observed experimentally (solid blue line) with the predictions of
Eq. (34) using $\sigma=6\,\Omega^{-1}$m-1 and
$\sigma_{1}=0.01\,\Omega^{-1}$m-1 (dashed red line), and with the predictions
of Eq. (40) using $\sigma=6\,\Omega^{-1}$m-1 (dotted green line), using the
east-west component of the magnetic field, which shows a magnetic field
variation about three times larger than the south-north component. Notice that
the geometry of the TAT-8 cable, which runs southwest to northeast, is such
that the south-north and east-west magnetic field components have components
of the same sign on the direction orthogonal to the cable.

As a general comment, the spread of the value of the conductivity $\sigma$ and
$\sigma_{1}$ in Figs. 3–5 is caused by the fact that the magnetic field to be
used in Eq. (34) is the spatial average of the component orthogonal to the
local cable direction, while we have used a single component, south-north in
Figs. 3 and 4 and east-west in Fig. 5, at a fixed location. The three cases
analyzed however show very clearly that the shape of the voltage variation is
not proportional to the time derivative of the magnetic field variation, which
would produce positive voltage variations in the leading edge and negative in
the trailing, but to the magnitude of its deviations. This is evident for the
good agreement between the observations and the shape predicted by Eq. (40),
especially for the shortest time scale events of Figs. 4 and 5, where the
approximations over which Eq. (40) is based are more accurate.

Figure 5: Left panel: East-west magnetic field horizontal component digitized
from Fig. 4 of ref. [2]. Right panel: the voltage variation in the TAT-8
system also obtained from Fig. 4 of ref. [2] (solid blue line), and the
prediction of Eq. (34) with $\sigma=6\,\Omega^{-1}$m-1 and
$\sigma_{1}=0.01\,\Omega^{-1}$m-1 (dashed red line) and of Eq. (40) with
$\sigma=6\,\Omega^{-1}$m-1 (dotted green line).

We have shown above that in both TAT-6 and TAT-8 systems the sign of the
electromotive force is consistent with the prediction of the theory. It is
interesting to notice that the choice of the direction of the current produces
qualitative differences in the sensitivity of a cable to geomagnetic storms.
In particular, the same increase of the south-north component would produce an
increase of the voltage or a decrease depending upon the direction of the
current. A large increase of the voltage is the worst scenario, because the
power supply may not be able to deliver the excess voltage required to feed
the in-line components with the prescribed current, whereas it can easily
reduce its voltage when necessary. Our finding that the voltage is
proportional to the deviations of the magnetic field rather than to their time
derivative implies that during a positive or negative magnetic field peak the
voltage does not show both positive and negative voltage variations. Our
analysis gives indications on the choice of the positive and negative side of
the line to reduce the sensitivity of the systems to geomagnetic storms.
Indeed, in the strongest geomagnetic storm recorded to date, the Carrington
event in 1859 [4], the magnetic field underwent a sharp decrease, and this is
the case also of the event that occurred on the 13th of March 1989 reported in
Fig. 5. In addition, there is experimental evidence that the strongest
geomagnetic storms will be associated to a weakening of the geomagnetic field
[3]. Therefore, it would be worth designing cables where a decrease of the
magnetic field is associated to a decrease of the required voltage, which does
not harm system operations. In east-west cables, this would correspond to
having the positive end westbound and the negative (or zero) eastbound, with
the current in the cable going west-east. Note that if the direction of the
current in the TAT-8 cable had been east-west instead of west-east, powering
the cable during the magnetic storm of March 13 1989 would have required a
voltage of 200 V in excess of the nominal value.

##  V Conclusions

We have analyzed the response of the voltage of the power supply of a
transoceanic transmission system to geomagnetic field transients. We have
shown that the voltage perturbation are proportional to the magnitude of the
geomagnetic field perturbation and not to its time derivative. This result
lead us to suggest design criteria to reduce the sensitivity of fiber optics
cables to strong geomagnetic storm events.

## References

  * [1] S. A. Jyothi, “Solar superstorms: Planning for an internet apocalypse,” in Proceedings of the 2021 ACM SIGCOMM 2021 Conference, 2021 692–704 (2021).
  * [2] L. V. Medford, L. J. Lanzerotti, J. S. Kraus, and C. G. Maclennan, “Transatlantic earth potential variations during the March 1989 magnetic storms,” Geophysical Research Letters 16, 1145–1148 (1989).
  * [3] P. K. Mohanty, K. P. Arunbabu, T. Aziz, S. R. Dugad, S. K. Gupta, B. Hariharan, P. Jagadeesan, A. Jain, S. D. Morris, B. S. Rao, Y. Hayashi, S. Kawakami, A. Oshima, S. Shibata, S. Raha, P. Subramanian, and H. Kojima, “Transient Weakening of Earth’s Magnetic Shield Probed by a Cosmic Ray Burst,” Phys. Rev. Lett. 117, 171101 (2016).
  * [4] E. W. Cliver and W. F. Dietrich, “The 1859 space weather event revisited: limits of extreme activity,” J. Space Weather Space Clim. 3 A31 (2013).
  * [5] L. V. Medford, A. Meloni, L. J. Lanzerotti, and G. P. Gregori, “Geomagnetic induction on a transatlantic communications cable,” Nature 290, 392–393 (1981).
  * [6] J. R. Wait, “Theory of Magneto-Telluric Fields,” Journal of Research of the National Bureau of Standards-D. Radio Propagation 66D, 509–541 (1962).
  * [7] D. H. Boteler, and R. J. Pirjola, The magnetic and electric fields produced in the sea during geomagnetic disturbances, Pure Appl. Geophys. 160, 1695–1716 (2003).
  * [8] I. Podlubny, R. L. Magin, and I. Trymorush, “Niels Henrik Abel and the birth of fractional calculus,” Fractional Calculus and Applied Analysis 20, 1068–1075 (2017).


#  High Fidelity Quantum State Transfer by Pontryagin Maximum Principle

Nahid Binandeh Dehaghani  Fernando Lobo Pereira  SYSTEC - Research Center for
Systems and Technologies, FEUP - Faculty of Engineering, Porto University, Rua
Dr. Roberto Frias sn, i219, 4200-465 Porto, Portugal (e-mail: nahid@ fe.up.pt,
flp@fe.up.pt).

###### Abstract

High fidelity quantum state transfer is an essential part of quantum
information processing. In this regard, we address the problem of maximizing
the fidelity in a quantum state transformation process satisfying the
Liouville-von Neumann equation. By introducing fidelity as the performance
index, we aim at maximizing the similarity of the final state density operator
with the one of the desired target state. Optimality conditions in the form of
a Maximum Principle of Pontryagin are given for the matrix-valued dynamic
control systems propagating the probability density function. These provide a
complete set of relations enabling the computation of the optimal control
strategy.

###### keywords:

quantum optimal control, quantum fidelity measures, Pontryagin maximum
principle

††thanks: The authors acknowledge the support of FCT for the grant
2021.07608.BD, the ARISE Associated Laboratory LA/P/0112/2020, and the R$\&$D
Unit SYSTEC - Base - UIDB/00147/2020 and Programmatic - UIDP/00147/2020 funds,
and also the support of projects SNAP - NORTE-01-0145-FEDER-000085, HARMONY -
NORTE-01-0145-FEDER-031411, funded by ERDF — COMPETE2020 — FCT/MEC — PT2020 —
NORTE2020, POCI and PIDDAC.

##  1 Introduction

Quantum control theory, (Dehaghani and Lobo Pereira, 2021), studies how the
dynamics of an atomic or subatomic-level system can be manipulated by means of
appropriate external electromagnetic fields or forces, designated by control,
in order to maximize a given performance criterion to the system. For many
quantum control protocols, the control law is required to be open-loop, that
is, only a function of time and does not require state feedback data provided,
for example, by measurements. In this context, optimal control theory provides
powerful tools, (Glaser et al., 2015). By means of Quantum Optimal Control
(QOC) theory, it is possible to formulate quantum control problems in order to
seek a set of admissible controls satisfying the system dynamics while
minimizing a cost functional in order to obtain a control law. In other words,
QOC aims to compute the shape and sequence of control pulses to achieve a
given task in an optimum way. Moreover, the high versatility of the optimal
control problem formulation provides advantages such as the incorporation of
diverse experimental constraints or limitations, while the optimality leads to
the exploitation of physical limits of the driven dynamics, (Boscain et al.,
2021). The investigation of optimal controls to quantum systems has been done
through adaptation of traditional optimal control tools such as the
variational method, (Kogut and Leugering, 2011), Pontryagin maximum principle,
(Lin et al., 2019; Yang et al., 2017; D’alessandro and Dahleh, 2001),
convergent iterative algorithms, (Wilhelm et al., 2020), among others.

The power of QOC has been instrumental in the design of several experiments,
such as preparation of motional states of Bose-Einstein condensate with
optimized control sequences, (van Frank et al., 2014), and loading improvement
of ultracold atoms in an optical lattice, (Rosi et al., 2013), to name just a
few. Pulses computed by using QOC allowed for error-resistant single qubit
gates, (Timoney et al., 2008). QOC of a single qubit also led to the design of
high-dynamic-range imaging of nanoscale magnetic fields, (Häberle et al.,
2013). Virtual transitions and leakage outside the qubit manifold was also
overcome in quantum processors based on superconducting circuits by means of
QOC results, (Lucero et al., 2010; Motzoi et al., 2009). QOC has proved to be
a key tool for quantum engineering in complex Hilbert spaces, (Larrouy et al.,
2020). These experiments and several others were preceded by a large number of
theoretical considerations on how QOC may impact on the control design of
quantum operations including state preparation, (Rojan et al., 2014; Günther
et al., 2021), ranging from squeezed states, (Grond et al., 2009), and
cluster-states, (Fisher et al., 2009), to many-body entangled state, (Caneva
et al., 2012), state transfer problem, (Ying-Hua et al., 2016; Zhang et al.,
2016; Guo et al., 2018), and quantum gate synthesis, (Huang and Goan, 2014;
Chou et al., 2015; Berrios et al., 2012).

Amongst the above-mentioned problems, quantum state transfer has received
considerable attention due to the high importance of trustworthy information
transfer in quantum networks. In this regard, high-quality quantum information
transport is crucial for practical models used in quantum computation,
(DiVincenzo, 1995). A protocol for performing fast and high fidelity quantum
state transfer in quantum spin $\displaystyle-\frac{1}{2}$ chains has been
proposed in (Zhang et al., 2016), by using QOC based on the Krotov algorithm.
Another research has shown that QOC allows a success rate of more than $98\%$
in an arbitrary state transfer process for a non-Markovian system implemented
with an optimal control, (Ying-Hua et al., 2016). However, despite a number of
researches based on QOC algorithms, the Pontryagin Maximum Principle is still
far from being fully exploited for quantum systems, and, even more for the
case when the dynamics described the evolution of the matrix-valued
probability density function are considered. In this context, we consider a
closed system described by Liouville-von Neumann equation, and formulate the
control problem of quantum state transfer aiming at the maximization of the
fidelity criterion. Moreover, we derive a shooting algorithm to solve the two
point boundary value problem entailed by the Maximum Principle of Pontryagin
by using a multiple shooting method.

The paper is organized as follows. In Section 2, after a brief description of
pure and mixed quantum states, we describe the models used for closed quantum
systems. The current state of the art concerning the use of the Pontryagin
Maximum Principle in the context of quantum systems are described in Section
3. In the next section, we formulate a quantum state transfer in the context
of optimal control subjected to Liouville-von Neumann equation with the aim of
maximizing fidelity. In section 5, we present an iterative algorithm to solve
the optimal control problem. The article ends with brief conclusions on the
reported research as well as a brief overview on prospective research
challenges.

##  2 Quantum Control Modelling

In the mathematical formulation of quantum mechanics, the state of a closed
quantum system is described by a unit vector $\left|\psi\right\rangle$ in a
complex Hilbert space $\mathbb{H}$ by using the Dirac representation, (Dirac,
1981). In quantum information theory, information is coded by means of qubits,
which can be is represented by

$\left|\psi\right\rangle=\cos\frac{\theta}{2}\left|0\right\rangle+{{e}^{i\varphi}}\sin\frac{\theta}{2}\left|1\right\rangle$
(1)

where $\theta\in\left[0,\pi\right]$, and $\varphi\in\left[0,2\pi\right]$.
Here, $\left|0\right\rangle$, and $\left|1\right\rangle$ correspond to the
states $0$ and $1$ for a classical bit, respectively, (Paul, 2007). Such
quantum state represented by the wave function $\left|\psi\right\rangle$ is
called a pure state.

In the density operator formalism, quantum states are described by density
operators $\rho:\mathbb{H}\to\mathbb{H}$ on the system’s Hilbert space
$\mathbb{H}$. For an ensemble
$\left\\{{{p}_{j}},\left|{{\psi}_{j}}\right\rangle\right\\}$ of pure states,
the density operator is defined as

$\rho\equiv\sum\limits_{j}{{{p}_{j}}\left|{{\psi}_{j}}\right\rangle\left\langle{{\psi}_{j}}\right|}$
(2)

in which
$\left\langle{{\psi}_{j}}\right|={{\left(\left|{{\psi}_{j}}\right\rangle\right)}^{\dagger}}$
and $\sum\limits_{j}{{{p}_{j}}=1}$, (Cong, 2014). For a pure quantum state,
$\rho=\left|\psi\right\rangle\left\langle\psi\right|$ and
$tr\left({{\rho}^{2}}\right)=1$, while for a mixed state
$tr\left({{\rho}^{2}}\right)<1$. Here $tr\,A$ denotes the trace of the matrix
$A$. An arbitrary state $\rho$ of a qubit can be written as a linear
combination of the so-called Pauli matrices
$\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$, which provide a basis for
$2\times 2$ self-adjoint matrices as

$\rho=\frac{1}{2}\left(I+\vec{r}.\vec{\sigma}\right)$ (3)

where the real vector $\textbf{r}=\left({{r}_{x}},{{r}_{y}},{{r}_{z}}\right)$
forms the coordinates of a point within the Bloch Sphere, and $I$ indicates
the $2\times 2$ identity matrix. Hence, for pure states $\left\|r\right\|=1$,
while $\left\|r\right\|<1$ represents mixed states, (Wilde, 2011).

There are several approaches for modelling a quantum system to be controlled.
In the Schrödinger model of quantum mechanics, bilinear models, including
Schrödinger and quantum Liouville equations, are used to describe closed
quantum systems. Schrödinger equation implies the evolution of the state
vector $\left|\psi(t)\right\rangle$ as

$i\hbar\frac{\partial}{\partial
t}\left|\psi(t)\right\rangle={H(t)}\left|\psi(t)\right\rangle,\quad\left|\psi(t)\right\rangle_{|t=0}=\left|{{\psi}_{0}}\right\rangle$
(4)

where ${H(t)}$, the Hamiltonian of the system, is a Hermitian operator on
$\mathbb{H}$, and $\hbar$ is the reduced Planck’s constant, considered as a
unit for convenience. The system control can be realized by a set of control
functions $u_{k}(t)\in\mathbb{R}$, which is coupled to the quantum system via
time independent interaction Hamiltonians ${{H}_{k}}(k=1,2,\dots)$. Therefore,
the total Hamiltonian
$H\left(t\right)={{H}_{d}}+\sum\limits_{k}{{{u}_{k}}(t)}{{H}_{k}}$ determines
the controlled evolution. Here, ${{H}_{d}}$ indicates the drift Hamiltonian,
and the ${{H}_{k}}$’s are the interaction Hamiltonians.

For systems in mixed states, quantum statistical ensemble can be characterized
via the statistical operator
$\rho\left(t\right)=U\left(t,t_{0}\right)\rho\left({{t}_{0}}\right){{U}^{\dagger}}\left(t,t_{0}\right)$,
in which
$\rho\left({{t}_{0}}\right)=\sum\limits_{\alpha}{{{w}_{\alpha}}\left|{{\psi}_{\alpha}}\left({{t}_{0}}\right)\right\rangle}\left\langle{{\psi}_{\alpha}}\left({{t}_{0}}\right)\right|$
with ${{w}_{\alpha}}$ indicating positive weights and
$\left|{{\psi}_{\alpha}}\left({{t}}\right)\right\rangle$ being the normalized
state vector evolving in time according to (4). By differentiating
$\rho\left(t\right)$ with respect to time, the equation of motion for density
operator, referred to as Liouville-von Neumann equation,

$\dot{\rho}\left(t\right)=\frac{-i}{\hbar}\left[{H}\left(t\right),\rho\left(t\right)\right]$
(5)

is obtained, where $\rho\left(t\right)$ is the variable to be controlled.
Equation (5) can be written in a form analogous to the classical Liouville
equation as
$\dot{\rho}\left(t\right)=\mathcal{L}\left(t\right)\rho\left(t\right)$, in
which $\mathcal{L}$ is Liouville super-operator, (Cong, 2014; Breuer et al.,
2002).

##  3 Short Overview on Quantum Optimal Control

The question that we address is: how we can drive (4) from an initial state
$\left|{{\psi}_{0}}\right\rangle$ to a desired target state
$\left|{{\psi}_{f}}\right\rangle$ while minimizing the cost functional? There
are several methods to design the optimal controller in quantum systems, which
vary according to the choice of the cost function, the construction of the
Pontryagin-Hamilton function, and the computation scheme using the Maximum
Principle conditions, (Cong, 2014).

A general optimal control problem can be formulated as follows: Given a set
$\mathcal{X}$ of state functions $x:\mathbb{R}\to{{\mathbb{R}}^{n}}$, and a
set $\mathcal{U}$ of control functions $u:\mathbb{R}\to{{\mathbb{R}}^{m}}$,
find the functions $x\in\mathcal{X}$ and $u\in\mathcal{U}$, which minimize a
cost functional $J:\mathcal{X}\times\mathcal{U}\to\mathbb{R}$ and satisfy the
dynamical constraint $\dot{x}=f(x,u)$. We can formulate any optimal control
problem as a specific case of the mentioned general problem, (d’Alessandro,
2021). Hence, for a system with state vector $x$ driven by controls $u$ over a
certain (fixed or variable) time interval $\left[0,T\right]$, the scalar real-
valued objective functional $J$ in the form of a problem of Bolza, is
expressed by

$J:=\varPhi\left(x\left(T\right),T\right)+\int\limits_{0}^{T}{L\left(x\left(t\right),u(t),t\right)dt}$
(6)

with $\varPhi$ and $L$ smooth functions
${{\mathbb{R}}^{n}}\times\mathbb{R}\to\mathbb{R},{{\mathbb{R}}^{n}}\times{{\mathbb{R}}^{m}}\times\mathbb{R}\to\mathbb{R}$,
respectively. The task is to maximize $J$ subjected to the condition that the
system dynamics are satisfied, with $x\left(0\right)={{x}_{0}}$, and $u(t)$
restricted to the set of admissible controls $\mathcal{U}$, (Dehaghani and
Lobo Pereira, 2021).

For many practical control problems in quantum setting, optimal control theory
is a powerful tool to achieve quantum control objectives, (Peirce et al.,
1988; Werschnik and Gross, 2007). The key concept of optimal control is that
the control law will be obtained via minimizing a cost function that drives
system (4) from
$\left|{{\psi}(0)}\right\rangle=\left|{{\psi}_{0}}\right\rangle$ to the
desired $\left|{{\psi}(T)}\right\rangle=\left|{{\psi}_{f}}\right\rangle$. Let
consider the problem of determining the control fields
$u\in{{L}^{2}}\left(\mathbb{C},\left[0,T\right]\right)$ while satisfying
equation (4). Suppose that we want to fulfill the following optimal criteria:

  * •

The control sequence brings the system at time $T$ to the desired state
${{\psi}_{d}}\in{{\mathbb{C}}^{n}}$.

  * •

The limited laser resources are taken into account through a minimization of
the control field effort.

  * •

The population of intermediate states that suffer strong environment losses
may need to be suppressed.

The above mentioned constraints can be summarized in the cost functional,
(Peirce et al., 1988; James, 2021; Borzi et al., 2002),

$\displaystyle{{J}(\left|\psi\right\rangle,u)}=\frac{1}{2}\left\|\left|\psi\left(T\right)\right\rangle-\left|{{\psi}_{f}}\right\rangle\right\|_{{{\mathbb{C}}^{n}}}^{2}$
(7) $\displaystyle\hskip
36.98866pt+\frac{\gamma}{2}\left\|u\right\|_{{{L}^{2}}\left(\mathbb{C},\left[0,T\right]\right)}^{2}+\frac{1}{2}\sum\limits_{j=1}^{n}{{{\alpha}_{j}}}\left\|{{\psi}_{j}}\right\|_{{{L}^{2}}\left(\mathbb{C},\left[0,T\right]\right)}^{2}$

in which $\gamma>0$ and ${{\alpha}_{j}}\geq 0$ are the weighting factors. We
also may add further constraints in the same way. The problem of optimal
control for minimizing cost functional (6) subject to equation(4) under normal
conditions admits a solution
$\left(\bar{\psi},\bar{u}\right)\in{{H}^{1}}\left({{\mathbb{C}}^{n}},\left[0,T\right]\right)\times{{L}^{2}}\left(\mathbb{C},\left[0,T\right]\right)$,
(Peirce et al., 1988; James, 2021; Borzi et al., 2002; Fattorini et al.,
1999). In order to calculate the necessary optimality conditions of first
order, the method of Lagrange multipliers is used. Hence, the Lagrangian
function is defined as

$L\left(\left|\psi\right\rangle,\left|p\right\rangle,u\right)=J\left(\left|\psi\right\rangle,u\right)+\operatorname{Re}\big{\langle}p,i\dot{\psi}-H\psi\big{\rangle}$
(8)

where
$\left\langle\phi,\psi\right\rangle=\int\limits_{0}^{T}{\phi\cdot}{{\psi}^{*}}dt$.
Here, “*” means the complex conjugate and “$\cdot$” is the usual vector-scalar
product in ${{\mathbb{C}}^{n}}$. Consider the minimization problem

$\displaystyle
L(|{\tilde{\psi}}\rangle,\left|\tilde{p}\right\rangle,\tilde{u})=\underset{\left|\psi\right\rangle\in
X^{0},\left|p\right\rangle,u\in{U}}{\mathop{\inf}}\,L\left(\left|\psi\right\rangle,\left|p\right\rangle,u\right)$
(9) $\displaystyle
U={{L}^{2}}\left({{\mathbb{C}}^{n}},\left[0,T\right]\right)$ $\displaystyle
X^{0}={{X}}\bigcap\left\\{\psi:\psi(0)={{\psi}_{0}}\right\\}$

The necessary conditions for problem (9) are obtained by equating to zero the
Fréchet derivative of $L$ with respect to $\left|\psi\right\rangle$,
$\left|p\right\rangle$, and $u$, Dehaghani and Lobo Pereira (2021).

Therefore, the optimality system entails:

$\displaystyle|{\dot{\psi}}\rangle$ $\displaystyle=$
$\displaystyle-i\left({{H}_{d}}+{{H}_{k}}u\right)\left|\psi\right\rangle,\;\;\left|\psi\left(0\right)\right\rangle=\left|{{\psi}_{0}}\right\rangle\vspace{.2cm}$
$\displaystyle\left|{\dot{p}}\right\rangle$ $\displaystyle=$
$\displaystyle-i\left(H_{d}+{{H}_{k}}u\right)\left|p\right\rangle-q,\vspace{.2cm}$
$\displaystyle\left|p\left(T\right)\right\rangle$ $\displaystyle=$
$\displaystyle-i\left(\left|\psi(T)\right\rangle-\left|{{\psi}_{d}}\right\rangle\right)\vspace{.2cm}$
(10) $\displaystyle u$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma}\operatorname{Re}\left[\left|p\right\rangle{{\left(\frac{\partial\left({{H}_{d}}+{{H}_{k}}u\right)}{\partial{{u}_{Re}}}\left|\psi\right\rangle\right)}^{*}}\right]\vspace{.2cm}$
$\displaystyle\hskip
28.45274pt+i\frac{1}{\gamma}\operatorname{Re}\left[\left|p\right\rangle{{\left(\frac{\partial\left({{H}_{d}}+{{H}_{k}}u\right)}{\partial{{u}_{Im}}}\left|\psi\right\rangle\right)}^{*}}\right]$

where ${{q}_{j}}={{\alpha}_{j}}{{\psi}_{j}}$, and $u={{u}_{Re}}+i{{u}_{Im}}$,
(Cong, 2014; Borzi et al., 2002; Zhu and Rabitz, 1998).

##  4 Formulation of the Optimal Control Problem

Consider the following optimal control problem $(P_{1})$

Minimize $\displaystyle-\mathcal{F}(\rho,\sigma)$ subject to
$\displaystyle\dot{\rho}(t)=F\left(\rho(t),u(t)\right)$
$\displaystyle\rho\left(0\right)={{\rho}_{0}}\in{{\mathbb{R}}^{n\times n}}$
$\displaystyle
u\left(t\right)\in\mathcal{U}:=\left\\{u\in{{L}_{\infty}}:u\left(t\right)\in\Omega\subset{{\mathbb{R}}^{m}}\right\\}$

where $\dot{\rho}:=\displaystyle\frac{d\rho}{dt}$, $t\in\left[0,1\right]$
determines the time variable, $\rho$ is the state variable in
${{\mathbb{R}}^{n}}\times{{\mathbb{R}}^{n}}$ supposed to satisfy the
differential constraints $\dot{\rho}=F\left(\rho,u\right)$ according to
Liouville-von Neumann equation (5), and $\rho_{0}$ is the so-called initial
quantum state. $u(\cdot)$ is the measurable bounded function termed as
control. Here, we aim to maximize fidelity $\mathcal{F}$ so that the density
operator $\rho$ has the maximum overlap with the target $\sigma$. The most
widely-used generalization of fidelity that has been indicated in the
literature is the Uhlmann-Jozsa fidelity, (Jozsa, 1994), which represents the
maximal transition probability between the purification of a pair of density
matrices, $\rho$ and $\sigma$, (Liang et al., 2019), and is defined as

$\mathcal{F}\left(\rho,\sigma\right):=\underset{\left|\psi\right\rangle,\left|\varphi\right\rangle}{\mathop{\max}}\,{{\left|\left\langle\psi\left|\varphi\right.\right\rangle\right|}^{2}}={{\left(tr\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)}^{2}}$
(11)

satisfying the following properties, (Liang et al., 2019),

  * i

$\mathcal{F}(\rho,\sigma)=\mathcal{F}(\sigma,\rho)$

  * ii

$0\leq\mathcal{F}\left(\rho,\sigma\right)\leq
1\quad\forall\;\rho,\,\sigma,\quad\mathcal{F}\left(\rho,\rho\right)=1$

  * iii

$\mathcal{F}\left(\rho,\sigma\right)=tr\left(\rho\sigma\right)$ if either
$\rho$ or $\sigma$ is a pure state.

  * iv

$F\left(U\rho{{U}^{\dagger}},U\sigma{{U}^{\dagger}}\right)=F\left(\rho,\sigma\right)$
for all unitary operations $U$.

The fidelity criteria signifies a security level for quantum state
transformation or the effectiveness of quantum gate synthesis, so is of high
importance in quantum systems.

##  5 Necessary Conditions of Optimality in the form of a Maximum Principle

The Pontryagin-Hamilton function ${\mathcal{H}}$ is defined for almost all
$t\in[0,T]$ by introducing the adjoint variable $\pi$, which is a time-varying
multiplier vector designated by costate or adjoint variable of the system.
Thus,

${\mathcal{H}}\left(\rho,u,\pi\right)=tr\left(\pi^{\dagger}F\left(\rho,u\right)\right)$
(12)

According to the Pontryagin’s Maximum Principle, for the optimal state
trajectory $\rho^{*}$ and the corresponding adjoint variable, the matrix
$\pi$, the optimal control $u^{*}(t)$, maximizes the Pontryagin-Hamiltonian
function ${\mathcal{H}}$ , i.e., for almost all $t\in[0,T]$,

${\mathcal{H}}\left(\rho^{*}(t),u,\pi(t)\right)\leq{\mathcal{H}}\left(\rho^{*}(t),u^{*}(t),\pi(t)\right)$
(13)

for all admissible control values $u\in\Omega$. Additionally, the adjoint
equation, and its terminal conditions imply that, Lebesgue a.e.,

$(-\dot{\pi}^{\dagger}(t),\dot{\rho}^{*}(t))=\nabla_{(\rho,\pi)}{\mathcal{H}}(\rho^{*}(t),u^{*}(t),\pi(t))$
(14) $\pi^{\dagger}(1)=\nabla_{\rho}\mathcal{F}(\rho^{*}(1),\sigma(1))$ (15)

From (14) we can obtain

$\displaystyle-\dot{\pi}^{\dagger}(t)$ $\displaystyle=$
$\displaystyle\nabla_{\rho}tr(\pi^{\dagger}F(\rho^{*}(t),u^{*}(t)))$ (16)
$\displaystyle=$
$\displaystyle\frac{-i}{\hbar}\frac{\partial}{\partial\rho}tr(\pi^{\dagger}(t)\left[H^{*}(t),\rho(t)\right])$
$\displaystyle=$
$\displaystyle\frac{-i}{\hbar}(\pi^{\dagger}(t)H^{*}(t)-H^{*}(t)\pi^{\dagger}(t))$
$\displaystyle=$
$\displaystyle\frac{i}{\hbar}\left[H^{*}(t),\pi^{\dagger}(t)\right].$

Remark that the dependence of $H$ on $t$ is through the control variable, and,
hence, $H^{*}(t)$ denotes $H$ evaluated at each time along the optimal control
$u^{*}$. It is straighforward to conclude that the differential equation (16)
has the formal solution

$\pi^{\dagger}(t)=e^{i\int_{t}^{1}H^{*}(s)ds}\pi^{\dagger}(1)e^{-\int_{t}^{1}H^{*}(s)ds}$
(17)

and the boundary condition at the final time is obtained by computing

$\pi^{\dagger}(1)=\nabla_{\rho}\left(tr\sqrt{\sqrt{\rho(1)}\sigma(1)\sqrt{\rho(1)}}\right)^{2}$
(18)

for which by means of Taylor expansion at the point $\rho=I$, we have

$\sqrt{\rho}=\sum_{k=0}^{\infty}\frac{1}{k!}\frac{d^{k}}{d\rho^{k}}\sqrt{\rho}_{|\rho=I}(\rho-I)^{k}$
(19)

which, in turn, under some conditions, the application of the Cayley–Hamilton
theorem yields, for a certain coefficients $\alpha_{k}$, $k=0,\ldots,n-1$,
being $n$ the dimension of the square matrix $\rho$,

$\sqrt{\rho}=\sum_{k=0}^{n-1}\alpha_{k}(\rho-I)^{k}.$

Thus, by using matricial calculus, we obtain

$\displaystyle\pi^{\dagger}(1)$ $\displaystyle=$ $\displaystyle
2tr\sqrt{\rho(1)\sigma(1)}\nabla_{\rho}(tr\sqrt{\rho(1)\sigma(1)})$ (20)
$\displaystyle=$ $\displaystyle
2tr\sqrt{\rho(1)\sigma(1)}\sum\limits_{k=0}^{n-1}{\alpha}_{k}\sum_{i=0}^{k-1}\bar{\rho}(1)^{i}\sqrt{\sigma(1)}\bar{\rho}(1)^{k-i-1}$

where $\bar{\rho}(1)=\rho(1)-I$.

##  6 Application of the Pontryagin Maximum Principle

Let us consider a spin $\displaystyle-\frac{1}{2}$ particle in an invariable
magnetic field $B_{0}$ along z-axis as the controlled system, the control
magnetic fields on the x–y plane is given by

$\displaystyle\gamma{{B}_{x}}=u\left(t\right)\cos\left(\omega t+\phi\right)$
(21) $\displaystyle\gamma{{B}_{y}}=-u\left(t\right)\sin\left(\omega
t+\phi\right)$

where $\gamma$ is the magnetic ratio of spin particle, and $u(t)$ is a real
valued number indicating the Rabi frequency of the particle. Hence, The system
Hamiltonian $H(t)$, composing of the free Hamiltonian $H_{d}$ and the control
Hamiltonian $H_{c}$, is expressed as

$\displaystyle H\left(t\right)=$
$\displaystyle\frac{-\gamma\hbar}{2}\left({{B}_{0}}{{\sigma}_{z}}+{{B}_{x}}\left(t\right){{\sigma}_{x}}+{{B}_{y}}\left(t\right){{\sigma}_{y}}\right)$
(22) $\displaystyle=$
$\displaystyle-\frac{\hbar}{2}\left(\begin{matrix}{{\omega}_{0}}&u(t){e}^{i\left(\omega
t+\phi\right)}\\\ u(t){{e}^{-i\left(\omega t+\phi\right)}}&-{{\omega}_{0}}\\\
\end{matrix}\right)$

in which $\gamma{{B}_{0}}=\omega_{0}$, and ${\sigma}_{x}$, ${\sigma}_{y}$, and
${\sigma}_{z}$ are the so-called Pauli matrices given by

${{\sigma}_{x}}=\left(\begin{matrix}0&1\\\ 1&0\\\
\end{matrix}\right),\quad{{\sigma}_{y}}=\left(\begin{matrix}0&-i\\\ i&0\\\
\end{matrix}\right),\quad{{\sigma}_{z}}=\left(\begin{matrix}1&0\\\ 0&-1\\\
\end{matrix}\right)$ (23)

Once the matrix $H(t)$ is known, the matrix $U(t)$ of this system can be
computed by

$U(t)={e}^{-\frac{i}{\hbar}\int_{0}^{t}H(s)ds}$ (24)

The system dynamics are given by the Liouville equation whose solution is
analogous to the one of the adjoint system considered in a previous section,
that is,

$\rho(t)=U\left(t\right)\rho(0){{U}^{\dagger}}\left(t\right)$ (25)

Here, we consider a state transfer problem from $\left|0\right\rangle$ to
$\left|1\right\rangle$, so $\rho\left(0\right)=\left(\begin{matrix}1&0\\\
0&0\\\ \end{matrix}\right)$ and
$\sigma\left(1\right)=\left(\begin{matrix}0&0\\\ 0&1\\\ \end{matrix}\right)$.

In this section, we present a time discretized computational scheme to solve
the optimal control problem associated with $(P_{1})$ by using an indirect
method based on the Maximum Principle. Let $N$ be the number of discrete time
subintervals. Given the smooth properties of the problems data, let us
consider a uniform discretization. Thus, we consider the $N$ points
$\displaystyle t_{k}=\frac{k}{N}$ for $k=0,\ldots,N-1$ and denote the value of
any function $f(t_{k})$ by $f_{k}$. We consider the both the system dynamics,
and the adjoint differential equations approximated by a first order Euler
approximation. Higher order methods yield better approximations but our option
is made to keep the presentation simple. Let $j=0,\ldots$ be the iterations
counter, and we denote the $j^{th}$ iteration of the function $f$ at time
$t_{k}$ by $f_{k}^{j}$. The proposed algorithm is as follows:

  * Step 1 - Initialization.

Let $j=0$, and initialize the values of $u_{k}^{j}$ for $k=1,\ldots,N-1$,
$\omega^{j}$, and $\phi^{j}$ in (22).

  * Step 2 - Computation of the state trajectory.

For $k=0,\dots,N-1$, let
$\displaystyle\rho_{k+1}^{j}=\rho_{k}^{j}+\frac{1}{k}\dot{\rho}_{k}^{j}$.

  * Step 3 - Computation of the adjoint trajectory.

Compute $\pi_{N}^{j}$ by using (20) with $\rho_{N}^{j}$ computed in Step 2.

Compute $\pi_{0}^{j}$ by using the discretized version of (17), that is,
$\displaystyle\pi_{0}^{j}=U^{j}{\pi_{N}^{j}}{U^{j}}^{\dagger}$ where
$\displaystyle U^{j}=e^{i\sum_{k=0}^{N-1}\frac{1}{N}H_{k}^{j}}$, being
$H_{k}^{j}$ the Hamiltonian of the system dynamics with the value control
$u^{j}$ at time $t_{k}$.

For $k=0,\dots,N-1$, let
$\displaystyle\pi_{k+1}^{j}=\pi_{k}^{j}+\frac{1}{k}\dot{\pi}_{k}^{j}$,

  * Step 4 - Computation of the Pontryagin Hamilton function

For $k=0,\dots,N-1$, let

$\displaystyle\mathcal{H}_{k}^{j}(u,\omega,\phi)=tr(-i{\pi^{j}_{k}}^{\dagger}[H(u,\omega,\phi,t_{k}),\rho_{k}^{j}])$

where $\displaystyle H(u,\omega,\phi,t_{k})$, from (22), is given by

$\displaystyle-\frac{1}{2}\left(\begin{matrix}{{\omega}_{0}}&u{e}^{i\left(\omega
t_{k}+\phi\right)}\\\ u{{e}^{-i\left(\omega
t_{k}+\phi\right)}}&-{{\omega}_{0}}\\\ \end{matrix}\right)$

  * Step 5: Update the control function.

For $k=0,\ldots,N-1$, compute the values $u_{k}^{j+1}$, $\omega^{j+1}$, and
$\phi^{j+1}$ that maximize the map

$(u,\omega,\phi)\to\mathcal{H}_{k}^{j}(u,\omega,\phi).$

  * Step 6: Stopping test.

For given small positive numbers $\varepsilon_{\omega}$, $\varepsilon_{\phi}$,
and $\varepsilon_{u}$ (tolerance errors), check whether all the following
inequalities hold:

$\displaystyle\hskip 42.67912pt|\omega^{j}-\omega^{j+1}|<\varepsilon_{\omega}$
$\displaystyle\hskip 42.67912pt|\phi^{j}-\phi^{j+1}|<\varepsilon_{\phi}$
$\displaystyle\hskip
14.22636pt\max_{k=0,\ldots,N-1}\\{|u_{k}^{j}-u_{k}^{j+1}|\\}<\varepsilon_{u}$

If yes, let, $\omega^{*}=\omega^{j+1}$, $\phi^{*}=\phi^{j+1}$, for
$k=0,\ldots,N-1$, $u^{*}(t_{k})=u_{k}^{j+1}$, and exit the algorithm.

Otherwise, let $j=j+1$, go to Step 2.

##  7 Conclusion

In this paper we have shown how the Maximum Principle of Pontryagin can be
applied in order to compute the optimal control of the problem of maximizing
fidelity in the transfer of the state variable, expressed by the matrix-valued
probability density function, between two given state values with the dynamics
of the system the Liouville-von Neumann equation. This context is not very
common neither in Optimal Control nor in Quantum Optimal Control. We obtained
a shooting algorithm to solve the two point boundary value problem resulting
from the application of the Maximum Principle of Pontryagin. Future challenges
consists in exploiting the versatility of the optimal control paradigm further
by including state constraints and other types of constraints. More efficient
algorithms to solve this problem will also be considered.

## References

  * Berrios et al. (2012) Berrios, E., Gruebele, M., Shyshlov, D., Wang, L., and Babikov, D. (2012).  High fidelity quantum gates with vibrational qubits.  _The Journal of Physical Chemistry A_ , 116(46), 11347–11354.
  * Borzi et al. (2002) Borzi, A., Stadler, G., and Hohenester, U. (2002).  Optimal quantum control in nanostructures: Theory and application to a generic three-level system.  _Physical Review A_ , 66(5), 053811.
  * Boscain et al. (2021) Boscain, U., Sigalotti, M., and Sugny, D. (2021).  Introduction to the pontryagin maximum principle for quantum optimal control.  _PRX Quantum_ , 2(3), 030203.
  * Breuer et al. (2002) Breuer, H.P., Petruccione, F., et al. (2002).  _The theory of open quantum systems_.  Oxford University Press on Demand.
  * Caneva et al. (2012) Caneva, T., Calarco, T., and Montangero, S. (2012).  Entanglement-storage units.  _New Journal of Physics_ , 14(9), 093041.
  * Chou et al. (2015) Chou, Y., Huang, S.Y., and Goan, H.S. (2015).  Optimal control of fast and high-fidelity quantum gates with electron and nuclear spins of a nitrogen-vacancy center in diamond.  _Physical Review A_ , 91(5), 052315.
  * Cong (2014) Cong, S. (2014).  _Control of quantum systems: theory and methods_.  John Wiley & Sons.
  * d’Alessandro (2021) d’Alessandro, D. (2021).  _Introduction to quantum control and dynamics_.  Chapman and hall/CRC.
  * D’alessandro and Dahleh (2001) D’alessandro, D. and Dahleh, M. (2001).  Optimal control of two-level quantum systems.  _IEEE Transactions on Automatic Control_ , 46(6), 866–876.
  * Dehaghani and Lobo Pereira (2021) Dehaghani, N.B. and Lobo Pereira, F. (2021).  Quantum control modelling, methods, and applications.  _work submitted for publication_.
  * Dirac (1981) Dirac, P.A.M. (1981).  _The principles of quantum mechanics_.  27\. Oxford university press.
  * DiVincenzo (1995) DiVincenzo, D.P. (1995).  Quantum computation.  _Science_ , 270(5234), 255–261.
  * Fattorini et al. (1999) Fattorini, H.O., Fattorini, H.O., et al. (1999).  _Infinite dimensional optimization and control theory_ , volume 54.  Cambridge University Press.
  * Fisher et al. (2009) Fisher, R., Yuan, H., Spörl, A., and Glaser, S. (2009).  Time-optimal generation of cluster states.  _Physical Review A_ , 79(4), 042304.
  * Glaser et al. (2015) Glaser, S.J., Boscain, U., Calarco, T., Koch, C.P., Köckenberger, W., Kosloff, R., Kuprov, I., Luy, B., Schirmer, S., Schulte-Herbrüggen, T., et al. (2015).  Training schrödinger’s cat: quantum optimal control.  _The European Physical Journal D_ , 69(12), 1–24.
  * Grond et al. (2009) Grond, J., Schmiedmayer, J., and Hohenester, U. (2009).  Optimizing number squeezing when splitting a mesoscopic condensate.  _Physical Review A_ , 79(2), 021603.
  * Günther et al. (2021) Günther, S., Petersson, N.A., and DuBois, J.L. (2021).  Quantum optimal control for pure-state preparation using one initial state.  _AVS Quantum Science_ , 3(4), 043801.
  * Guo et al. (2018) Guo, Y., Dong, D., and Shu, C.C. (2018).  Optimal and robust control of quantum state transfer by shaping the spectral phase of ultrafast laser pulses.  _Physical Chemistry Chemical Physics_ , 20(14), 9498–9506.
  * Häberle et al. (2013) Häberle, T., Schmid-Lorch, D., Karrai, K., Reinhard, F., and Wrachtrup, J. (2013).  High-dynamic-range imaging of nanoscale magnetic fields using optimal control of a single qubit.  _Physical review letters_ , 111(17), 170801.
  * Huang and Goan (2014) Huang, S.Y. and Goan, H.S. (2014).  Optimal control for fast and high-fidelity quantum gates in coupled superconducting flux qubits.  _Physical Review A_ , 90(1), 012318.
  * James (2021) James, M. (2021).  Optimal quantum control theory.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 4, 343–367.
  * Jozsa (1994) Jozsa, R. (1994).  Fidelity for mixed quantum states.  _Journal of modern optics_ , 41(12), 2315–2323.
  * Kogut and Leugering (2011) Kogut, P.I. and Leugering, G.R. (2011).  _Optimal control problems for partial differential equations on reticulated domains: approximation and asymptotic analysis_.  Springer Science & Business Media.
  * Larrouy et al. (2020) Larrouy, A., Patsch, S., Richaud, R., Raimond, J.M., Brune, M., Koch, C.P., and Gleyzes, S. (2020).  Fast navigation in a large hilbert space using quantum optimal control.  _Physical Review. X_ , 10(2).
  * Liang et al. (2019) Liang, Y.C., Yeh, Y.H., Mendonça, P.E., Teh, R.Y., Reid, M.D., and Drummond, P.D. (2019).  Quantum fidelity measures for mixed states.  _Reports on Progress in Physics_ , 82(7), 076001.
  * Lin et al. (2019) Lin, C., Wang, Y., Kolesov, G., and Kalabić, U. (2019).  Application of pontryagin’s minimum principle to grover’s quantum search problem.  _Physical Review A_ , 100(2), 022327.
  * Lucero et al. (2010) Lucero, E., Kelly, J., Bialczak, R.C., Lenander, M., Mariantoni, M., Neeley, M., O’Connell, A., Sank, D., Wang, H., Weides, M., et al. (2010).  Reduced phase error through optimized control of a superconducting qubit.  _Physical Review A_ , 82(4), 042339.
  * Motzoi et al. (2009) Motzoi, F., Gambetta, J.M., Rebentrost, P., and Wilhelm, F.K. (2009).  Simple pulses for elimination of leakage in weakly nonlinear qubits.  _Physical review letters_ , 103(11), 110501.
  * Paul (2007) Paul, T. (2007).  Quantum computation and quantum information.  _Mathematical Structures in Computer Science_ , 17(6), 1115–1115.
  * Peirce et al. (1988) Peirce, A.P., Dahleh, M.A., and Rabitz, H. (1988).  Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications.  _Physical Review A_ , 37(12), 4950.
  * Rojan et al. (2014) Rojan, K., Reich, D.M., Dotsenko, I., Raimond, J.M., Koch, C.P., and Morigi, G. (2014).  Arbitrary-quantum-state preparation of a harmonic oscillator via optimal control.  _Physical Review A_ , 90(2), 023824.
  * Rosi et al. (2013) Rosi, S., Bernard, A., Fabbri, N., Fallani, L., Fort, C., Inguscio, M., Calarco, T., and Montangero, S. (2013).  Fast closed-loop optimal control of ultracold atoms in an optical lattice.  _Physical Review A_ , 88(2), 021601.
  * Timoney et al. (2008) Timoney, N., Elman, V., Glaser, S., Weiss, C., Johanning, M., Neuhauser, W., and Wunderlich, C. (2008).  Error-resistant single-qubit gates with trapped ions.  _Physical Review A_ , 77(5), 052334.
  * van Frank et al. (2014) van Frank, S., Negretti, A., Berrada, T., Bücker, R., Montangero, S., Schaff, J.F., Schumm, T., Calarco, T., and Schmiedmayer, J. (2014).  Interferometry with non-classical motional states of a bose–einstein condensate.  _Nature communications_ , 5(1), 1–6.
  * Werschnik and Gross (2007) Werschnik, J. and Gross, E. (2007).  Quantum optimal control theory.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 40(18), R175.
  * Wilde (2011) Wilde, M.M. (2011).  From classical to quantum shannon theory.  _arXiv preprint arXiv:1106.1445_.
  * Wilhelm et al. (2020) Wilhelm, F.K., Kirchhoff, S., Machnes, S., Wittler, N., and Sugny, D. (2020).  An introduction into optimal control for quantum technologies.  _arXiv preprint arXiv:2003.10132_.
  * Yang et al. (2017) Yang, Z.C., Rahmani, A., Shabani, A., Neven, H., and Chamon, C. (2017).  Optimizing variational quantum algorithms using pontryagin’s minimum principle.  _Physical Review X_ , 7(2), 021027.
  * Ying-Hua et al. (2016) Ying-Hua, J., Ju-Ju, H., Jian-Hua, H., and Qiang, K. (2016).  Optimal control-based states transfer for non-markovian quantum system.  _Physica E: Low-dimensional Systems and Nanostructures_ , 81, 77–82.
  * Zhang et al. (2016) Zhang, X.P., Shao, B., Hu, S., Zou, J., and Wu, L.A. (2016).  Optimal control of fast and high-fidelity quantum state transfer in spin-1/2 chains.  _Annals of Physics_ , 375, 435–443.
  * Zhu and Rabitz (1998) Zhu, W. and Rabitz, H. (1998).  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.