# Connection between inverse engineering and optimal control in shortcuts to
adiabaticity

Qi Zhang  International Center of Quantum Artificial Intelligence for Science
and Technology (QuArtist) and Department of Physics, Shanghai University,
200444 Shanghai, China  Laboratoire Collisions, Agrégats, Réactivité, IRSAMC,
Université de Toulouse, CNRS, UPS, France  Xi Chen<EMAIL_ADDRESS>International Center of Quantum Artificial Intelligence for Science and
Technology (QuArtist) and Department of Physics, Shanghai University, 200444
Shanghai, China  Department of Physical Chemistry, University of the Basque
Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  David Guéry-Odelin
<EMAIL_ADDRESS>Laboratoire Collisions, Agrégats, Réactivité, IRSAMC,
Université de Toulouse, CNRS, UPS, France

###### Abstract

We consider fast high-fidelity quantum control by using a shortcut to
adiabaticity (STA) technique and optimal control theory (OCT). Three specific
examples, including expansion of cold atoms from the harmonic trap, atomic
transport by moving harmonic trap, and spin dynamics in the presence of
dissipation, are explicitly detailed. Using OCT as a qualitative guide, we
demonstrate how STA protocols designed from inverse engineering method, can
approach with very high precision optimal solutions built about physical
constraints, by a proper choice of the interpolation function and with a very
reduced number of adjustable parameters.

††preprint: APS/123-QED

##  I Introduction

The last ten years witnessed the huge development of “shortcuts to
adiabaticity” (STA) with wide applications ranging from atomic, molecular, and
optical physics (AMO) to quantum information transfer or processing review1 ;
review2 . The concept of STA was originally proposed to speed up the adiabatic
processes in quantum control. Nowadays, STA become versatile toolboxes for
controlling the dynamics and transformation in quantum physics review1 ;
review2 , statistical physics DGO14 ; Mart1 , integrated optics ShuoyenEPL ,
and classical physics Faure1 ; jarzynski1 ; jarzynski2 ; lianao . In this
context, the most popular STA techniques are the fast-forward scaling masuda1
; masuda2 , the counterdiabatic driving rice1 ; rice2  (or transitionless
quantum algorithm berry ; chenprl2 ; adofolprl ; deffner ), and the invariant-
based inverse engineering chenprl , and their variants. These three techniques
can be shown to be mathematically equivalent Xipra11 ; Erikpra12 . However,
the diversity of the designs of shortcut protocols or their combination may be
required for a realistic experimental implementation multiple . Furthermore,
some counterdiabatic hamiltonians turn out to be unfeasible 3STIRAP , or some
systems cannot be treated by means of invariant-based engineering.

STA method provides a useful toolbox for fast and robust quantum controls with
applications in a wide variety of quantum platforms such as cold atoms Schaff1
; Schaff2 , NV center spin NV1 ; NV2  including for their use as a quantum
sensor jorge , trapped ion Kim , and superconducting qubit Yingyi1 ; Yingyi2 ;
Antti ; Dapeng  to name a few.

Such controls have also a clear added value to quantum optimal control in
quantum information processing and quantum computing QuantumOC , in terms of
analytical tools, numerical tools, and a combination of these two. Numerical
optimal control such as the gradient ascent pulse engineering (GRAPE)
algorithm works to some extent as a black box. The dynamics and the structure
of the control field are not easily predictable Assemat . STA techniques based
on a clear physical picture deliver a more easily understandable framework but
are mostly addressing problems of low complexity. However, these techniques
have recently been combined with deep machine learning for more involved
physical problems HesonPNAS ; Nature ; Sels ; Ding .

Interestingly, shortcut protocols can be readily engineered to accommodate for
various physical constraints Larocca2020 , or to mitigate an environmental
noise. In this respect, the combination of inverse engineering methods and
optimal control theory (OCT) has been particularly fruitful Stefanatostime-
optimal ; StepanatosTF ; XiaojingPRA2014 ; Xipra11optimal ; Qipra ; QiJPB ;
PRL2013D ; njp2012 ; inverse13 . Most STA techniques provides solutions that
are robust against a small variation of the duration of the parameter
engineering. In Ref. Martikyan , it is shown how OCT solutions can be adapted
to accommodate for extra boundary conditions to ensure a similar robustness.
Alternatively, the STA technique of Ref. DavidMuga  provides an explicit
solution for linear control problems fulfilling the Kalman criterium
Martikyan2 .

In this article, we compare systematically the inverse engineering method with
the result of optimal control theory on three specific examples that can be
addressed analytically in both formalisms: expansion of cold atoms from the
harmonic trap, atomic transport by moving harmonic trap, and spin dynamics in
the presence of dissipation. Our aim is to provide a pedagogical introduction
and comparison between a simple if not the simplest Shortcut To Adiabaticity
technique, the direct inverse engineering of the equation of motion of the
dynamical variables, and the optimal control theory. STA techniques are built
about the boundary conditions while OCT involves the minimization of a cost
function. To facilitate the comparison we therefore discuss how inverse
engineered (IE) solutions can be modified in order to minimize a cost function
and mimic OCT solutions. Similarly to the variational method in quantum
mechanics, and as illustrated in the following, the family of functions over
which the minimization is performed play a crucial role. In the following, we
also show how a simple ansatz having just a few tunable parameters can
approach very precisely the optimal solution obtained for a given physical
constraint.

##  II Fast cooling in time-varying harmonic traps

Fast frictionless coooling for ultracold and Bose-Einstein condensates belongs
to the first experimental demonstrations of STA techniques Schaff1 ; Schaff2 .
Such techniques have been subsequently adapted and applied to cold-atom
mixtures Choi , Tonks-Girardeau gas CampoPRA11 ; deffner , Fermi gases Deng ;
DGO14 , and many-body systems Jorge15 .

In this section, we address the problem of fast atomic cooling in a time-
dependent harmonic trap chenprl . We derive the time-dependence of the trap
frequency by an inverse engineer procedure on an Ermakov equation and using
OCT. We subsequently compare the two types of solutions. Interestingly, the
tunability inherent to the inverse engineering method provides the required
flexibility to shape the inverse-engineered trajectories to minimize a cost
function. We show how such solutions can be simply adapted to get results very
close to optimal solutions for a time-averaged energy cost function
XiMugaPRA2010 ; chenprl ; Stefanatostime-optimal .

###  II.1 Model, Hamiltonian, and the Inverse Engineering Approach

More specifically, we consider in the following the fast decompression of a
one-dimensional (1D) harmonic potential from an initial angular frequency
$\omega(0)=\omega_{0}$ to the final target one $\omega(t_{f})=\omega_{f}$,
($\omega_{f}<\omega_{0}$). The problem amounts to finding the time-dependent
solution of the Schrödinger equation that ensures the transformation from the
ground state of the initial trap to the ground state of the final trap in a
finite amount of time $t_{f}$:

$i\hbar\frac{\partial\psi}{\partial
t}=\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{2}m\omega^{2}(t)x^{2}\right]\psi.$ (1)

For this purpose, we look for a scaling solution of the form
$\psi(x,t)=\exp[-\beta(t)]\exp[-\alpha(t)x^{2}]f(\rho(t)=x/b(t),t)$. The first
factor accounts for the normalization, the second factor for the evolution of
the phase (we show below that it is purely imaginary) and the last one for the
desired scaling dynamics. By plugging such an ansatz into the Schrödinger
equation, we find how the different parameters are related:

$\displaystyle i\hbar\partial_{t}f$ $\displaystyle=$
$\displaystyle\left(i\hbar\dot{\beta}+\frac{\hbar^{2}}{m}\alpha\right)f+\left(\frac{2\hbar^{2}}{m}\alpha+i\hbar\frac{\dot{b}}{b}\right)\rho\partial_{\rho}f$
(2) $\displaystyle+$
$\displaystyle\\!\\!\left(\\!i\hbar\dot{\alpha}\\!-\\!\frac{2\hbar^{2}}{m}\\!+\\!\frac{1}{2}m\omega^{2}\\!\right)\\!b^{2}\rho^{2}f-\frac{\hbar^{2}}{2mb^{2}}\partial_{\rho\rho}f.$

By introducing the renormalized time
$\tilde{t}(t)=\int^{t}_{0}dt^{\prime}/b(t^{\prime})^{2}$ and for the choice
$\alpha=(-im/2\hbar)\dot{b}/b$ and $\beta=\ln b/2$, the effective wave
function $\Psi(\rho,\tilde{t})=f(\rho,t)$ obeys a _time-independent_
Schrödinger equation:

$i\hbar\frac{\partial\Psi}{\partial\tilde{t}}=\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{2}m\omega^{2}_{0}\rho^{2}\right]\Psi,$
(3)

provided that the scaling parameter $b(t)$ satisfies the following Ermakov
equation

$\ddot{b}+\omega^{2}(t)b=\frac{\omega^{2}_{0}}{b^{3}}.$ (4)

Interestingly, this latter equation is amenable to a set of linear equations.
Indeed, it is the equation of an effective 2D oscillator in polar coordinates,
the $1/b^{-3}$ is nothing but the centrifugal barrier which acts as a
repulsive force that prohibits the access to a zero value of $b$.
Alternatively, the very same result can be obtained by using Lewis-Riesenfeld
dynamical invariant chenprl . The ground state wave function in such a time-
dependent harmonic trap reads

$\psi(x,t)=\frac{\mathscr{N}}{\sqrt{b}}\exp\left(\frac{im\dot{b}}{2\hbar
b}x^{2}\right)\exp\left(-\frac{x^{2}}{2a^{2}_{0}b^{2}}\right),$ (5)

where $\mathscr{N}$ accounts for the normalization and
$a_{0}=\sqrt{\hbar/(m\omega_{0})}$. The self-consistent boundary conditions
for a smooth continuous interpolation function are chenprl :

$\displaystyle b(0)$ $\displaystyle=$ $\displaystyle
1,\;\;\dot{b}(0)=0,\;\;\ddot{b}(0)=0,$ $\displaystyle b(t_{f})$
$\displaystyle=$
$\displaystyle\gamma=\sqrt{\omega_{0}/\omega_{f}},\;\;\dot{b}(t_{f})=0,\;\;\mbox{and}\;\;\ddot{b}(t_{f})=0.$
(6)

As a simple example, one can choose for the scaling factor $b(t)$ a fifth
order polynomial ansatz that fulfills the above six boundary conditions
chenprl :

$b(\tau)=1+(\gamma-1)(10\tau^{3}-15\tau^{4}+6\tau^{5}).$ (7)

In view of the comparison with optimal protocols, we calculate hereafter the
mean energy associated to the ground state wave function (5) XiMugaPRA2010 :

$\displaystyle\overline{E}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{t_{f}}\int_{0}^{t_{f}}E(t)dt=\frac{1}{t_{f}}\int_{0}^{t_{f}}\langle\psi(t)|H(t)|\psi(t)\rangle
dt$ (8) $\displaystyle=$
$\displaystyle\frac{\hbar}{2\omega_{0}}\frac{1}{t_{f}}\int_{0}^{t_{f}}\left(\dot{b}^{2}+\frac{\omega_{0}^{2}}{b^{2}}\right)dt,$

where $E(t)=K(t)+E_{p}(t)$ is the sum of the kinetic energy
$K(t)=\hbar(\dot{b}^{2}+\omega_{0}^{2}/b^{2})/(4\omega_{0})$ and the potential
energy $E_{p}=\hbar\omega^{2}(t)b^{2}/(4\omega_{0})$. The mean energies obey
the Virial theorem: $\overline{E_{p}}=\overline{K}=\overline{E}/2$. For any
$b(t)$ trajectory that fulfills the boundary conditions, one can infer
$\omega(t)$ from Eq. (4) and calculate explicitly the mean energies. Using the
available freedom to shape the scaling factor $b(t)$, the inverse-engineered
solutions can be tuned so to minimize the time-averaged energy as discussed in
section II.3.

###  II.2 Optimal control theory

Shortcut To Adiabaticity protocols such as inverse engineering are built about
the boundary conditions. We have provided an example using a polynomial
interpolation. Optimal control theory (OCT) offers an alternative to find a
path between two states but shall be built about a cost function. We propose
hereafter to use OCT on the Ermakov equation (4).

For this purpose, we recast Eq. (4) into a set of first order nonlinear
coupled equations, $\dot{\textbf{x}}=\textbf{f}(\textbf{x}(t),u)$ by defining
the x components as $x_{1}=b(t)$ and $x_{2}=\dot{b}/\omega_{0}$, and
introducing the (scalar) control function,
$u(t)=\omega^{2}(t)/\omega_{0}^{2}$:

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle x_{2},$ (9)
$\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle-
ux_{1}+\frac{1}{x_{1}^{3}}.$ (10)

In the following, we work out two OCT solutions associated to the minimization
of the final time and then of the mean energy. As a result of the nonlinear
character of the set of Hamiltonian equations, the Pontryagin maximum
principle only gives a necessary condition to get an extremum.

####  II.2.1 Time-optimal solution

The so-called time-optimal solution amounts to minimizing the cost function

$J=\int^{t_{f}}_{0}1dt.$ (11)

with the boundary conditions (6) which translates on the x vector components
as $x_{1}(0)=1$, $x_{2}(0)=0$ and $x_{1}(t_{f})=\gamma$ and $x_{2}(t_{f})=0$.
We furthermore choose the constraint $|u|\leq 1$ Hoffman ; Stefanatostime-
optimal . We stress that we let the possibility for the control parameter to
be either positive or negative. When it is negative, the curvature of the
harmonic confinement is reversed. Atoms are therefore transiently expelled
which provides a method to accelerate the desired transformation.

To minimize the cost function (11), we apply the Pontryagin maximum principle
which states that there exists non-zero, continuous vector p with components
($p_{0},p_{1},p_{2}$), fulfilling Hamilton’s equations Stefanatostime-optimal
; StepanatosTF : $\dot{\textbf{x}}=\partial H_{c}/\partial\textbf{p}$ and
$\dot{\textbf{p}}=-\partial H_{c}/\partial\textbf{x}$. With the cost function
$J$, the control Hamiltonian $H_{c}$ reads

$H_{c}=p_{0}+p_{1}x_{2}+p_{2}\left(-x_{1}u+\frac{1}{x_{1}^{3}}\right),$ (12)

where $p_{0}$ is a non-zero normalization constant, and $p_{1}$ and $p_{2}$
are generalized Lagrange multipliers. The Pontryagin’s maximum principle
states that at any instant ($0\leq t\leq t_{f}$), the values of the control
function $u$ maximize $H_{c}$. As $H_{c}$ is linear in the control function
$u$ and since $x_{1}>0$, the sign of the factor in front of $u$,
($-p_{2}x_{1}$) is fully determined by the sign of $-p_{2}$. This latter
parameter plays the role of a switching function for “bang-bang” type control
as discussed in the literature XiaojingPRA2014 ; Hoffman ; Stefanatostime-
optimal ; StepanatosTF . The fact that the Hamilton equations are nonlinear
enables the possibility to have multiple bang-bang solutions Stefanatostime-
optimal . We consider in the following the simplest solution with analytical
expression. This “bang-bang” solution has a single intermediate time (see Fig.
1):

Figure 1: Fast cooling in time-varying harmonic traps: The 3-jump “bang-bang”
control function, $u(t)=\omega^{2}(t)/\omega_{0}^{2}$.

$u(t)=\left\\{\begin{aligned} 1&,&~{}t\leq 0\\\
-(\omega_{1}/\omega_{0})^{2}&,&~{}0<t<t_{1}\\\
(\omega_{2}/\omega_{0})^{2}&,&~{}t_{1}<t<t_{f}\\\
(\omega_{f}/\omega_{0})^{2}&,&~{}t\geq t_{f}\end{aligned}\right.$ (13)

With such a control function, we infer the value of the scaling factor $b(t)$
from the Ermakov equation and find the following solution for “bang-bang”
control that fulfills the boundary conditions (6):

$\displaystyle
b(t)\\!=\\!\left\\{\begin{array}[]{ll}\sqrt{1+\frac{\omega_{1}^{2}+\omega_{0}^{2}}{\omega_{1}^{2}}\sinh^{2}(\omega_{1}t)},&0\\!\leq\\!t\\!\leq\\!t_{1}\\\
\sqrt{\gamma^{2}+\frac{\omega_{0}^{2}-\gamma^{4}\omega_{2}^{2}}{\gamma\omega_{2}^{2}}\sin^{2}[\omega_{2}(t_{f}-t)]},&t_{1}\\!\leq\\!t\\!\leq\\!t_{f}.\end{array}\right.$
(16)

It is worth noticing that the Ermakov equation implies that the quantity
$x_{2}^{2}+ux_{1}^{2}+x_{1}^{-2}=c$ is constant. The value of the constant $c$
is fixed by the initial conditions for $0<t<t_{1}$ and by the final conditions
for $t_{1}<t<t_{f}$. Using the continuity of $b(t)$ at $t_{1}$ and $t_{f}$ due
to the second derivative in the Ermakov equation, we find the explicit
expression for both times Hoffman ; Stefanatostime-optimal :

$\displaystyle t_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\omega_{1}}{\rm{arcsinh}}\sqrt{\frac{\omega_{1}^{2}(\gamma^{2}-1)(\gamma^{2}\omega_{2}^{2}-\omega_{0}^{2})}{\gamma^{2}(\omega_{1}^{2}+\omega_{0}^{2})(\omega_{2}^{2}+\omega_{1}^{2})}},$
(17) $\displaystyle t_{f}$ $\displaystyle=$ $\displaystyle
t_{1}+\frac{1}{\omega_{2}}\arcsin{\sqrt{\frac{\omega_{2}^{2}(\gamma^{2}-1)(\gamma^{2}\omega_{1}^{2}+\omega_{0}^{2})}{(\omega_{1}^{2}+\omega_{2}^{2})(\gamma^{4}\omega_{2}^{2}-\omega_{0}^{2})}}}.$
(18)

As the time $t_{1}$ shall remain real, we deduce from Eq. (17) that
$\omega_{2}\geq\omega_{0}/\gamma>\omega_{f}$. The last inequality is naturally
satisfied because of the cooling constraint $\omega_{0}>\omega_{f}$. The first
inequality requires $\omega_{0}/\gamma\leq\omega_{2}\leq\omega_{0}$. In Fig.
2, we plot the normalized final time $s_{f}=t_{f}\omega_{0}$ as a function of
$\omega_{1}/\omega_{0}$ and $\omega_{2}/\omega_{0}$ in their accessible
domains. We conclude that the shortest normalized final time $s_{f}$ is
obtained for the largest $\omega_{1}$ and $\omega_{2}$. With the choice
$\omega_{2}=\omega_{1}=\omega_{0}$, we obtain the shortest time

$s_{f}^{min}=\frac{\pi}{4}+\frac{1}{2}\ln\left(\frac{\omega_{0}}{\omega_{f}}\right).$
(19)

The lowest bound for $\omega_{2}$, namely, $\omega_{2}=\omega_{0}/\gamma$
provides the upper bound for final time

$s_{f}=\frac{\pi}{2}\gamma,$ (20)

where the first period of time is reduced to $t_{1}=0$, so that only two jumps
are needed.

Figure 2: Fast cooling in time-varying harmonic traps: 2D color plot of the
final normalized time $t_{f}\omega_{0}$ for a 3-jump “bang-bang” control as a
function of the first pulse amplitude $\omega_{1}/\omega_{0}$, and the second
pulse amplitude $\omega_{2}/\omega_{0}$.

Figure 3:  Fast cooling in time-varying harmonic traps: (a)Example of time-
optimal trajectory of $b(t)$ from Eq. (16). Parameters:
$\omega_{f}^{2}=\omega_{0}^{2}/5$. (b) The time-averaged energy as a function
of the normalized final time $s_{f}=\pi\gamma/2$, obtained from the time-
optimal control solution.

In this latter range of parameter, the scaling factor reads

$b(\tau)=\sqrt{\gamma^{2}+(1-\gamma^{2})\sin^{2}\left[\frac{\pi(1-\tau)}{2}\right]},$
(21)

with $\tau=t/t_{f}$. In Fig. 3 (a), we plot such an example of the time
evolution of $b(t)$. The solution that corresponds to the upper bound for the
final time also provides the minimum time-averaged energy. Using Eq. (8), we
calculate this latter quantity:

$\overline{E_{p}}=\frac{\varepsilon}{2}\left(1+\frac{1}{\gamma^{2}}\right)=\frac{\varepsilon}{2}\left(1+\frac{\pi^{2}}{4s_{f}^{2}}\right),$
(22)

where $\varepsilon=\hbar\omega_{0}/4$. In Fig. 3 (b), we plot this time-
averaged energy $\overline{E_{p}}$ as a function of the final time $s_{f}$. It
is worth noticing that $\omega_{f}$ and $t_{f}$ are not independent since
$s_{f}=\pi\gamma/2=\pi\sqrt{\omega_{0}/\omega_{f}}/2$.

####  II.2.2 Time-averaged energy minimization

In this section, we consider optimal control solution associated to the
minimization of time-averaged energy with unbounded constraint XiMugaPRA2010 .
The lower bound for the time-averaged potential (total) energy in Eq. (8)
reads XiMugaPRA2010 ; Cui :

$\displaystyle\overline{E_{p}}^{op}$ $\displaystyle=$
$\displaystyle\varepsilon\left[\left(\frac{B}{s_{f}}\right)^{2}-1-\frac{2}{s_{f}}\operatorname{arctanh}\left(\frac{B^{2}+B-s_{f}^{2}}{s_{f}}\right)\right.$
(23) $\displaystyle+$
$\displaystyle\left.\frac{2}{s_{f}}\operatorname{arctanh}\left(\frac{B}{s_{f}}\right)\right],$

with the following solution of
$b(\tau)=\sqrt{(B^{2}-s_{f}^{2})\tau^{2}+2B\tau+1}$ and
$B=-1+\sqrt{s_{f}^{2}+\gamma^{2}}$. In Fig. 4, we plot this lower bound for
optimized time-averaged energy as a blue dashed line.

Figure 4: Fast cooling in time-varying harmonic traps: Comparison of time-
averaged potential energy for different optimal protocols: (1) energy-
minimization (blue dashed line), (2) time-optimal protocol (with $\omega_{f}$
fixed, $\overline{E_{p}}$ constant) (red point); and inversed-engineered
protocols: (1) 0-freedom polynomial in Eq. (7) (black dotted line), (2)
polynomial IE solution with two free parameters optimized for a given
normalized final time (stars) (see Table 1): $s_{f}=1.1$ (green line),
$s_{f}=\pi\gamma/2$ (red solid line), and $s_{f}=4$ (orange solid line).
Parameters: $\omega_{f}^{2}=\omega_{0}^{2}/5$.

###  II.3 Comparison between IE and OCT

In the previous subsections, we have reviewed the streamline of IE and OCT
protocols to ensure a fast frictionless decompression in a harmonic trap whose
strength can be time-engineered. As already discussed, there is a lot of
freedom to design inverse-engineered protocols since the only requirements
concern the boundary conditions. However, the question of the mean energy cost
of such protocols may be relevant since a real potential always exhibits some
anharmonicity when the potential energy becomes too large. In what follows, we
propose to design IE protocols having a minimal mean potential energy. We will
show how we can readily approach the optimal results.

The IE solution exhibited in Eq. (7) relies on a fifth-order polynomial that
fulfills the six boundary conditions. In Fig. 4, we plot the corresponding
mean potential energy $\overline{E_{p}}(s_{f})$ using a black dotted line
which turns out to be quite far from the optimal solution (dashed blue line).

To reduce $\overline{E_{p}}(s_{f})$, we remove the constraints on $\dot{b}$
and $\ddot{b}$ at initial and final time since they are not strictly speaking
necessary neither fulfilled by the optimal solution. We also enlarge the
parameter space for $b(\tau)$ using a third-order polynomial ansatz
$b(\tau)=\sum\nolimits_{n=0}^{3}a_{n}\tau^{n}$ to keep some free parameters.
The two boundary conditions yields $a_{0}=1$ and
$a_{1}=-1-a_{2}-a_{3}+\gamma$. For different normalized final time $s_{f}$, we
can therefore minimize the time-averaged energy with respect to the two
parameters $a_{2}$ and $a_{3}$. In Table 1, we provide the optimal values
$a_{2}$ and $a_{3}$ that minimizes the mean potential energy for the three
cases with $s_{f}=1.1$, $s_{f}=\pi\gamma/2$, and $s_{f}=4$.

Figure 5: Fast cooling in time-varying harmonic traps: Comparison of time-dependent normalized variable $b(\tau)$ obtained from optimal control theory (averaged energy optimization) (blue dashed line) and from inverse engineered solutions optimized to minimize the time-averaged energy (red solid line) for different final times (a) $s_{f}=1.1$, (b) $s_{f}=\pi\gamma/2$, and (c) $s_{f}=4$. The corresponding optimal values of polynomial functions are detailed in Table 1. Parameters: $\omega_{f}^{2}=\omega_{0}^{2}/5$. $s_{f}$ |  $a_{2}$ |  $a_{3}$
---|---|---
1.1 |  -0.44893 |  0.10996
$\pi\gamma/2$ |  -1.47741 |  0.34535
4 |  -2.86194 |  0.62841
Table 1:  Optimal values of the free parameters $a_{2}$ and $a_{3}$ in the
three-order polynomial ansatz for the IE protocol that minimize the time-
averaged energy. Parameter $\omega_{f}^{2}=\omega_{0}^{2}/5$.

The results are represented as stars in Fig. 4. They nearly coincide with the
result of the optimal control theory. This is confirmed by plotting the
scaling functions for both protocols (see Fig. 5). We conclude that the IE
trajectories inspired by the OCT solutions can be readily designed to approach
with an impressive accuracy the exact OCT solutions.

##  III Fast transport of atoms in moving harmonic traps

STA techniques have also been applied to high-fidelity fast quantum transport
of neutral atoms DGO08  or charged ions Walther12 ; Bowler12  using a moving
trap. Such developments have a wide range of applications from quantum
information processing Erikion ; ErikBEC  to atom fountain clock, atom chip
manipulation CorgierNJP ; Becker ; Amri  or atomic interferometry Dupont . In
recent closely related works, optimal trajectories that minimize the
excitation in ion shuttling in the presence of stochastic noise have been
designed by combining invariant-based inverse engineering, perturbation
theory, and optimal control Xiaojingpra ; Xiaopra18 .

In this section, we address the problem of the fast transport of a single atom
based on a moving 1D harmonic potential. The particle is supposed to be
initially in the ground state and shall remain in the ground state at the
final time. We follow the same kind of presentation as previously: We first
design inverse-engineered protocols, we then derive the OCT protocols for time
Xipra11optimal  and mean-energy optimization QiJPB , and eventually compare
both approaches.

###  III.1 Classical and quantum inverse-engineered solutions

The time-dependent Hamiltonian of atomic transport using a moving harmonic
trap reads

$H(t)=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}_{0}[\hat{x}-x_{0}(t)]^{2},$
(24)

where $\omega_{0}$ is the constant trap angular frequency, and $x_{0}(t)$ the
time-dependent position of the trap center. This problem amounts to finding
the appropriate driving of this harmonic oscillator. The exact mapping between
the classical and quantum solutions enables one to solve the classical problem
to get a solution valid quantum mechanically DavidMuga . The time-evolution of
the coordinate, $x(t)$, of a classical particle under the time-dependent
Hamiltonian (24) is given by

$\ddot{x}+\omega_{0}^{2}(x-x_{0}(t))=0.$ (25)

A smooth perfect transport i.e. a transport without any residual oscillations
at final can be obtained using inverse engineering by imposing the six
boundary conditions:

$\displaystyle x(0)=x_{0}(0)=0,\;\;\dot{x}(0)=0,\;\;\ddot{x}(0)=0,\;\;$
$\displaystyle
x(t_{f})\\!=\\!x_{0}(t_{f})\\!=\\!d,\;\;\dot{x}_{0}(t_{f})\\!=\\!0,\;\;\mbox{and}\;\;\ddot{x}_{0}(t_{f})\\!=\\!0.$
(26)

Any interpolation function $x(t)$ that fulfills these boundary conditions
provides a possible solution of our problem. For instance, one can use the
following fifth order polynomial interpolation function:

$x(t)=d\left[10\left(\frac{t}{t_{f}}\right)^{3}-15\left(\frac{t}{t_{f}}\right)^{4}+6\left(\frac{t}{t_{f}}\right)^{5}\right].$
(27)

Once $x(t)$ is known, the trajectory of the trap center $x_{0}(t)$ can be
directly inferred from Eq. (25). A similar result can be derived quantum
mechanically using the properties of dynamical invariants Erikion ;
Xipra11optimal . In view of the optimization that we will perform later on, it
is worth working out the instantaneous average potential energy

$\langle V(t)\rangle=\frac{\hbar\omega_{0}}{4}+E_{p}(t),$ (28)

where the first term accounts for the zero-point energy contribution and
$E_{p}(t)=\frac{1}{2}m\omega_{0}^{2}(x(t)-x_{0}(t))^{2}$ i.e. the
instantaneous potential energy for the effective classical particle. The time-
averaged potential energy is defined by

$\overline{E_{p}}=\frac{1}{t_{f}}\int^{t_{f}}_{0}E_{p}(t)dt.$ (29)

###  III.2 Optimal control theory

To recast this problem as an optimal problem, we define the variables
$x_{1}(t)=x(t)$ and $x_{2}(t)=\dot{x}$, and the control function
$u(t)=x(t)-x_{0}(t)$. The control function corresponds to the relative
position of the effective particle with respect to the trap center. The
equation of motion (25) for the effective particle can be encapsulated in the
following set of linearly coupled first order differential equations
$\dot{\textbf{x}}=\textbf{f}[\textbf{x}(t),u]$:

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle x_{2},$ (30)
$\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle-\omega_{0}^{2}u.$
(31)

Interestingly, for this linear system, the solution deduced from the
Pontryagin formalism provides the unique control solution $u$ that minimizes
the cost function.

####  III.2.1 Time minimization

In this section, we solve the time-optimal problem with an upper bound on the
relative displacement $|u|\leq\delta$. The cost function to minimize $t_{f}$
is

$J=\int^{t_{f}}_{0}1dt.$ (32)

The corresponding Pontryagin Hamiltonian reads
$H_{c}=p_{0}+p_{1}x_{2}-\omega^{2}p_{2}u$, where the Lagrange multipliers
$p_{1}$ and $p_{2}$ fulfill $\dot{p}_{1}=0$ and $\dot{p}_{2}=-p_{1}$. We
deduce $p_{1}=c_{1}$ and $p_{2}=-c_{1}t+c_{2}$ where $c_{1}$ and $c_{2}$ are
constants to be determined. The Hamiltonian $H_{c}$ is a linear function of
the bounded control function $u(t)$. As a result, the sign of $p_{2}$ sets the
sign of $u(t)$ to maximize $H_{c}$. The parameter $p_{2}$ being a linear
function of time, the sign of $p_{2}$ can only change once. By considering the
initial and final boundary conditions, the appropriate control sequence taking
into account the upper bound for $|u(t)|$ is a (three-jump) “bang-bang”
control

$\displaystyle u(t)=\left\\{\begin{array}[]{lll}0,&t\leq 0\\\
-\delta,&0<t<t_{1}\\\ \delta,&t_{1}<t<t_{f}\\\ 0,&t\geq
t_{f}.\end{array}\right.$ (37)

With such a control function, the time-optimal solution of Eq. (25) compatible
with the boundary conditions (26) reads

$\displaystyle x(t)=\left\\{\begin{array}[]{lll}0,&t\leq 0\\\ \omega^{2}\delta
t^{2}/2,&0<t<t_{1}\\\ d-\omega^{2}\delta(t-t_{f})^{2}/2,&t_{1}<t<t_{f}\\\
d.&t\geq t_{f}.\end{array}\right.$ (42)

The driving of the trap bottom is then given by
$x_{0}(t)=\ddot{x}(t)/\omega^{2}_{0}+x(t)$. By imposing, the continuity on
$x(t_{1})$ and $\dot{x}(t_{1})$, one gets the explicit expression for the
switching and final times:

$t_{1}=\frac{t_{f}}{2},~{}~{}~{}t_{f}=\frac{2}{\omega_{0}}\sqrt{\frac{d}{\delta}}.$
(43)

According to Eq. (29), the time-averaged potential energy $\overline{E_{p}}$
for this constrained protocol is

$\overline{E_{p}}=\frac{8md^{2}}{\omega_{0}^{2}{t_{f}}^{4}}=\frac{1}{2}m\omega_{0}^{2}\delta^{2}.$
(44)

####  III.2.2 Mean potential energy minimization

In this section, we work out the energy-optimal protocol. We here provide a
solution that minimizes the time-averaged potential energy for a given
transport time $t_{f}$ and distance $d$, with unbounded constraint. According
to the definition of potential energy,
$E_{p}=\frac{1}{2}m\omega^{2}_{0}(x-x_{0})^{2}$, the cost function for this
problem is

$J=\frac{1}{2}m\omega^{2}_{0}\int^{t_{f}}_{0}u^{2}dt,$ (45)

and the Pontryagin Hamiltonian

$H_{c}=\frac{1}{2}m\omega^{2}_{0}p_{0}u^{2}+p_{1}x_{2}-p_{2}\omega^{2}_{0}u.$
(46)

The Hamilton equations give two costate equations similar to those derived in
the previous section. For the normalization, we can choose the constant
parameter $p_{0}=-1/m$, so that the optimal problem amounts to maximizing the
quantity $-u^{2}/2-p_{2}u$.

For convenience, we consider the unbounded case ($u(t)$ is unbounded) which
sets the lowest bound for time-averaged potential energy $\overline{E_{p}}$.
The quantity $-u^{2}/2-p_{2}u$ is maximal for $u=-p_{2}$. This expression for
the control fonction combined to Eq. (25) and the boundary conditions (26)
enables one to determine the optimal trajectory of the center of mass:

$x(t)=\frac{dt^{2}}{{t_{f}}^{2}}\left(3-2\frac{t}{t_{f}}\right),$ (47)

from which we infer the trap center trajectory $x_{0}(t)$ using Eq. (25) with
initial and final boundary conditions $x_{0}(0)=0$ and $x_{0}(t_{f})=d$:

Figure 6:  Fast transport of atoms in a moving harmonic trap: Comparison of
the trajectories of (a) the center of mass $x(t/t_{f})/d$ and (b) the trap
center $x_{0}(t/t_{f})/d$, obtained from the OCT formalism by minimizing the
time-averaged potential energy (blue dashed line) and using the IE approach
(red solid line) based on a fifth-order polynomial ansatz. Parameters:
$\omega_{0}=2\pi\times 50$ Hz and $t_{f}=22$ ms.

$\displaystyle x_{0}(t)\\!=\\!\left\\{\begin{array}[]{lll}0,&t\leq 0\\\
\displaystyle\\!\\!\\!\left(\\!\\!1-\frac{2t}{t_{f}}\\!\right)\\!\frac{6d}{\omega_{0}^{2}{t_{f}}^{2}}\\!+\\!\\!\left(\\!\\!3\\!-\\!\frac{2t}{t_{f}}\\!\right)\\!\\!\frac{t^{2}d}{{t_{f}}^{2}},&0\\!<\\!t\\!<\\!t_{f}\\\
d,&t\geq t_{f}.\end{array}\right.$ (51)

In Fig. 6, we plot the OCT center of mass along with the bottom trap
trajectories for some specific values using blue dashed lines. It is worth
noticing that according to our optimal solution the trap center has to include
two sudden jumps at initial and final time. With such an optimization
performed for an unbounded control function, we get the following lowest time-
averaged potential energy

$\overline{E_{p}}^{(OCT)}=\frac{6md^{2}}{\omega_{0}^{2}{t_{f}}^{4}}.$ (52)

In Fig. 7, we also plot this minimal time-averaged potential energy as a
function of the final time $t_{f}$ as a blue dashed line.

Figure 7:  Fast transport of atoms in a moving harmonic trap: Time-averaged
potential energy $\overline{E_{p}}/\varepsilon$ (normalized to
$\varepsilon=m\omega_{0}^{2}d^{2}/2$) as a function of final time $t_{f}$ by
using different protocols: time-optimal (orange dash-dotted line), energy-
minimization with unbounded constraint (blue dashed line), and IE approaches
with a fifth-order polynomial (black upper solid line), a seventh-order
polynomial (purple solid line), and nineteenth-order polynomial (red lower
solid line). Same parameters as Fig. 6.

###  III.3 Comparison between IE and OCT

In this section, we use the freedom in the interpolation function that enters
IE solutions to approach the solution of the optimal control theory associated
to a minimization of the time-averaged potential energy.

####  III.3.1 IE with polynomial ansatzs

For this purpose, we first enlarge the parameter space of the polynomial
ansatz that fulfills the boundary conditions (26) and search for the optimal
values of the coefficients that minimize the time-averaged potential energy.

To satisfy the six boundary conditions (26) the minimal order of the
polynomial interpolation function is five (see Eq. (27)). In Fig. 6, we plot
the center of mass, $x(t/t_{f})/d$, and trap center, $x_{0}(t)$, trajectories
as a function of time using red solid lines. The corresponding time-averaged
potential energy is $\overline{E_{p}}^{(P5)}=1.42\overline{E_{p}}^{(OCT)}$
which is significantly larger than the minimal potential energy given by Eq.
(52). It is represented as a black solid line in Fig. 7.

In order to further reduce the time-averaged potential energy, we enlarge the
parameter space, while keeping the six boundary conditions satisfied. We
search for a solution of the forme
$x(t)=d[\sum\nolimits_{n=0}^{7}a_{n}(t/t_{f})^{n}$]. By applying the boundary
conditions (26), we have $a_{0}=a_{1}=a_{2}=0$, $a_{5}=21-6a_{3}-3a_{4}$,
$a_{6}=-35+8a_{3}+3a_{4}$, and $a_{7}=15-3a_{3}-a_{4}$. The time-averaged
potential energy can be explicitly worked out:

Figure 8:  Fast transport of atoms in a moving harmonic trap: Comparison of
trajectories of mass of center (a) and trap center (b), calculated from the
OCT formalism (blue dashed line) and the IE approach (red solid line) with a
19th order polynomial ansatz. Same parameters as Fig. 6.

$\displaystyle\overline{E_{p}}(a_{3},a_{4})$ $\displaystyle=$
$\displaystyle\left[7+\frac{17}{77}(a_{3}-21)^{2}+\frac{4}{385}(a_{4}+70)^{2}\right.$
(53) $\displaystyle+$
$\displaystyle\left.\frac{(a_{3}-21)(a_{4}+70)}{11}\right]\frac{md^{2}}{\omega_{0}^{2}{t_{f}}^{4}}.$

The minimization of this energy yields $a_{3}=21$ and $a_{4}=-70$ and
$\overline{E_{p}}^{(P7)}\simeq 1.16\overline{E_{p}}^{(OCT)}$. This curve as a
function of the final time is represented in Fig. 7 as a purple solid line. It
provides a clear improvement with respect to the fifth-order polynomial
solution. A priori, it is possible to further improve the optimization using a
higher order polynomial ansatz. For instance, using a 19th order well-
optimized polynomial, we have found $\overline{E_{p}}^{(P19)}\simeq
1.018\overline{E_{p}}^{(OCT)}$. In Fig. 8, we have plotted the corresponding
time-dependent trajectories $x_{c}(t)$ and $x_{0}(t)$. However, the increase
of the polynomial order requires a minimization with an increasing number of
parameters. This is somehow cumbersome. In the following section, we propose
another type of interpolating function inspired by the OCT solution and
yielding astonishing results.

####  III.3.2 IE with hyperbolic ansatz

Figure 9:  Fast transport of atoms in a moving harmonic trap: Comparison of
trajectories of mass of center (a) and trap center (b), calculated from the
OCT formalism (blue dashed line) and the IE approach with the optimized
hyperbolic-function protocol in Eq. (54) (red solid line). The “magic” values
are $a_{1}=1.2$, $a_{2}=1.25$, and the other parameters are the same as those
in Fig. 6.

In this subsection, we apply the IE approach using the following hyperbolic-
function

$x(t)=\frac{d}{2}\tanh\left\\{a_{1}\tan\left[\frac{\pi}{a_{2}}\left(\frac{t}{t_{f}}-\frac{1}{2}\right)\right]\right\\}+\frac{d}{2},$
(54)

where $a_{2}>1$ to avoid any singularity. Interestingly, the choice of the
parameter $a_{2}$ enables one to mimic a jump at initial and final time. This
class of solution with the possibility of an initial and final offsest and
with similar symmetry as the optimal function provides a very performant class
of functions for the optimization. The freedom provided by the two parameters
$a_{1}$ and $a_{2}$ enables one to reduce the time-averaged potential energy
while satisfying the two boundary conditions $x(0)=0$ and $x(t_{f})=d$. Such
an optimization gives $a_{1}=1.2$ and $a_{2}=1.25$. The corresponding
trajectories $x(t)$ and $x_{0}(t)$ are plotted in Fig. 9, and the mean
potential energy is represented in Fig. 10 with marked red points. This ansatz
provides a solution that nearly coincides with the exact solution,
$\overline{E_{p}}^{(hyp)}\simeq 1.0001\overline{E_{p}}^{(OCT)}$.

Figure 10:  Fast transport of atoms in a moving harmonic trap: Time-averaged
potential energy $\overline{E_{p}}(t_{f})/\varepsilon$ (normalized to
$\varepsilon=m\omega_{0}^{2}d^{2}/2$) calculated from different protocols:
time-optimal (orange dash-dotted line), energy-minimization with unbounded
constraint (blue dashed line), and IE approach based on a hyperbolic-function-
ansatz by choosing the “magic” values $a_{1}=1.2$ and $a_{2}=1.25$ (marked red
point). Same parameters as Fig. 7.

For this transport problem, we have shown how the freedom on the interpolation
ansatz enables one to optimize extra constraints such as the mean energy
whilst fulfilling the boundary conditions. The choice of the ansatz has a
strong impact. One could naively think that a very high order polynomial could
always provide a succesful strategy. However, we have shown on this example
that the convergence may be quite slow with the degree of the polynomial, and
that the investigation of other shapes with a few adjustable parameters can
easily outperform the polynomial interpolation for a give constraint.

##  IV Spin dynamics in the presence of dissipation

In contrast with the previous sections, we address in the following an example
dealing with the control of internal degrees of freedom. Optimal control
provides a powerful tool to solve time-optimal and energy-optimal problems in
quantum two-level and three-level systems JMP ; Damme ; Damme18 ;
SpinStefanatos2009 . Such result can be directly extended to two uncoupled
Damme18  and coupled StefanatosPRA19  spins with similar approach. Using
numerical optimal algorithm, robust optimal control can also be designed that
accounts for inhomogeneous boarding and/or dissipation Owrutsky ; Damme ; Lin
; ShersonNJP . Inverse engineering techniques have also been used for the fast
and robust control of single spin QiScirep  and two-interacting spins Xiaotong
; QiScirep  in the presence of dissipation impens . Systematic error or
perturbation induced from the parameter fluctuatiosn, dephasing noise, bit
flip can be further suppressed using IE and OCT in atomic population transfer
PRL2013D ; njp2012 ; inverse13  and spin flip Xiaotong .

Figure 11: Spin dynamics in the presence of dissipation: Equivalent magnetic
field ($B,B_{c}$) of transverse magnetic field $B_{\perp}$.

Strictly speaking, the presence of dissipation rules out the possibility of an
adiabatic evolution. However, the inverse engineering can still be applied. In
the following, we consider the control of a spin 1/2 (${\bf
S}=(S_{x},S_{y},S_{z})$through the appropriate design of the time-varying
magnetic field components (${\bf B}=(B_{x},B_{y},B_{z})$) for the desired
boundary conditions. More precisely, we address the dissipative evolution of
this spin in the presence of a strong transverse relaxation rate, $R>0$. As is
commonly the case in NMR, the longitudinal relaxation rate is supposed to be
negligible compared to the transverse one, and is here neglected
SpinStefanatos2009 . Under those assumptions, the spin components obey the
Bloch equations:

$\displaystyle\dot{S_{x}}$ $\displaystyle=$ $\displaystyle-RS_{x}-B_{y}S_{z},$
(55) $\displaystyle\dot{S_{y}}$ $\displaystyle=$ $\displaystyle-
RS_{y}+B_{x}S_{z},$ (56) $\displaystyle\dot{S_{z}}$ $\displaystyle=$
$\displaystyle B_{y}S_{x}-B_{x}S_{y}.$ (57)

Following Ref. SpinStefanatos2009 , we recast the Bloch equations using
spherical coordinates. For this purpose, we introduce the angles $\theta(t)$
and $\phi(t)$ such that ${\bf
S}=(r\sin\theta\cos\varphi,r\sin\theta\sin\varphi,r\cos\theta)$ where $r$
denotes de length of the spin $r=\sqrt{S_{x}^{2}+S_{y}^{2}+S_{z}^{2}}$. It is
convenient to decompose the transverse magnetic field
$B_{\perp}=(B_{x},B_{y})$ into $B_{\perp}=(B,B_{c})$, satisfying $B\parallel
S_{\perp}$ and $B_{c}\perp S_{\perp}$ (see Fig. 11):
$B=(B_{x}/R)\cos\phi-(B_{y}/R)\sin\phi$ and
$B_{c}=(B_{x}/R)\sin\phi+(B_{y}/R)\cos\phi$. The Bloch equations can be
readily rewritten with the variables $a=\ln r$,
$\tan\theta=\sqrt{S_{x}^{2}+S_{y}^{2}}/S_{z}$, $\tan\phi=S_{y}/S_{x}$, and the
normalized time $t=Rt^{\prime}$:

$\displaystyle\dot{a}$ $\displaystyle=$ $\displaystyle-[\sin\theta(t)]^{2},$
(58) $\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle
B-\sin\theta(t)\cos\theta(t),$ (59) $\displaystyle\dot{\phi}$ $\displaystyle=$
$\displaystyle B_{c}\cot\theta(t).$ (60)

To ensure a spin rotation from an initial spin-up state to a given final
target state, we shall use the boundary conditions

$\displaystyle\theta(0)$ $\displaystyle=$ $\displaystyle
0,\;\;a(0)=0,\;\;\theta(t_{f})=\theta_{f},\;\;\mbox{and}$ $\displaystyle
a(t_{f})$ $\displaystyle=$ $\displaystyle
a_{f}=-\int^{t_{f}}_{0}[\sin\theta(t)]^{2}dt.$ (61)

It is worth emphasizing the fact that choosing the final spin length $a_{f}$
and orientation $\theta_{f}$ for a given final time $t_{f}$ may have no
solution for finite resources. Indeed, if the driving by the magnetic field is
not sufficiently strong, the dissipation will set an upper limit on the final
spin length.

The field component $B_{c}$ is always perpendicular to $r_{\perp}$ and
therefore only affects the spin rotation about the $z$-axis. The angle
$\theta$ is responsible for the partial or total spin flip. To minimize the
energy cost, the trajectory length shall be minimal. This latter condition
sets the value of $B_{c}$ to zero which means $\phi=constant$. Basically, the
IE technique amounts here to fixing the $\theta(t)$ function in accordance
with the boundary conditions (61), and inferring the external magnetic field
$B(t)$ from Eq. (59).

###  IV.1 Energy minimization by OCT

We consider here a given spin manipulation from $(a(0)=0,\theta(0)=0)$ to
$(a_{f},\theta_{f})$ with the minimum magnetic field amplitude. For this
purpose, we aim at minimizing the cost function

$E=\int^{t_{f}}_{0}\frac{B(t)^{2}}{2}dt.$ (62)

Let’s first recast this problem as a control problem involving a set of
coupled first order equations. By defining the state variables $x_{1}=a$,
$x_{2}=\theta$, and the control function $u(t)=B(t)$, the system equations
(58) and (59) is of the form $\dot{\textbf{x}}=\textbf{f}(\textbf{x}(t),u)$:

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle-\sin^{2}x_{2},$
(63) $\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle u-\sin
x_{2}\cos x_{2},$ (64)

and the cost function is

$J=\int^{t_{f}}_{0}\frac{u(t)^{2}}{2}dt.$ (65)

The corresponding Pontryagin Hamiltonian reads

$H_{c}=-\frac{1}{2}u^{2}-p_{1}\sin^{2}x_{2}+p_{2}(u-\sin x_{2}\cos x_{2}),$
(66)

where $p_{1}$ and $p_{2}$ are the Lagrange multipliers fulfilling
$\dot{\textbf{p}}=-\partial H_{c}/\partial\textbf{x}$ i.e. $\dot{p}_{1}=0$,
$\dot{p}_{2}=p_{1}\sin(2x_{2})+p_{2}\cos(2x_{2})$. The maximum Pontryagin
principle states for an unbounded control $u$ that $\partial H_{c}/\partial
u=0$, i.e. $u=p_{2}$. In the absence of terminal cost, the optimal solution
for this optimization between fixed initial and final states but without
fixing the final time gives the extra condition
$H_{c}[\textbf{p}(t),\textbf{x}(t),u(t)]=0$:

$p_{2}=(\sqrt{2p_{1}+\cos^{2}x_{2}}+\cos x_{2})\sin x_{2}.$ (67)

From Eq. (64), we deduce

$\dot{x}_{2}=\sin x_{2}\sqrt{\cos^{2}x_{2}+2p_{1}}.$ (68)

By combining Eqs. (68) and Eq. (63), we find $dx_{1}=-\sin
x_{2}/\sqrt{2p_{1}+\cos^{2}x_{2}}dx_{2}$. After integration, this relation
gives

$r(\theta)=\frac{\cos\theta+\sqrt{2p_{1}+\cos^{2}\theta}}{1+\sqrt{2p_{1}+1}}.$
(69)

The (constant) value of $p_{1}$ is deduced self-consistently with the boundary
conditions. The final time provided by OCT for an arbitrary target $r_{f}$ is
determined by

$t_{f}=\int^{t_{f}}_{0}dt=\int^{\theta_{f}}_{0}\frac{1}{\dot{\theta}(\theta)}d\theta.$
(70)

We note that the dissipation has an influence on the final time.

###  IV.2 Case I: reaching the horizontal plane of the Bloch sphere

In this subsection, we consider the transfer of the spin from the quantization
axis to the horizontal plane. The boundary conditions are thus $\theta(0)=0$,
$r(0)=1$ and $\theta(t_{f})=\theta_{f}=\pi/2$. This choice sets the value of
the constant $p_{1}$: $p_{1}^{\pi/2}=2{r_{f}}^{2}/(1-r_{f}^{2})^{2}$. To
address a specific example, we consider in the following the final value
$r(t_{f})=r_{f}=e^{-2}$. The final time obtained from Eq. (70) suffers from a
logarithmic divergence. To cure this problem, we shift the initial and final
time by a small quantity $\varepsilon\ll 1$: $\theta(0)=\epsilon$ and
$\theta_{f}=\pi-\epsilon$:

$t_{f}^{\pi/2}=\frac{1-r_{f}^{2}}{1+r_{f}^{2}}\left[\ln\left(\frac{1+r_{f}^{2}}{r_{f}}\right)-\ln\epsilon\right]=8.60481849$
(71)

for $\varepsilon=10^{-3}$. For this specific example, the cost function
associated to this optimal solution (see Eq. (65)) is

${E}^{(OCT)}_{\pi/2}=\frac{1}{1-{r_{f}}^{2}}=1.01866.$ (72)

Figure 12:  Spin dynamics in the presence of dissipation: Energy as a function
of $a_{f}=\ln r_{f}$ for the same target state
$(\theta_{f},r_{f},t_{f})=(\pi/2,0.6,3.6357955)$. We compare the results
obtained from the optimal control theory (red star) with the inverse
engineering results involving two different polynomial ansatz fulfilling the
boundary conditions. The energy curve are plotted for different values of the
polynomial coefficient $a_{1}$: (a) for a second order polynomial ansatz and
(b) for a third order polynomial ansatz with $a_{3}=0.1$. The inset in (b)
shows the proximity of the inverse engineering result with that of the optimal
control theory.

For comparison with the inverse engineering method, we propose, for the very
same $t_{f}$, the following second order polynomial ansatz:

$\theta(t)=a_{1}t-\frac{a_{1}t_{f}-\theta_{f}}{{t_{f}}^{2}}t^{2}.$ (73)

This ansatz fulfills the boundary conditions and has a single free parameter.
The corresponding cost function, $E^{(P2)}$, is minimal for $a_{1}=-0.119582$:
$E^{(P2)}=1.055E^{(OCT)}_{\pi/2}$.

Figure 13:  Spin dynamics in the presence of dissipation: For a $\pi/2$
rotation, we plot (a) the magnetic field $B(t)$ and (b) the corresponding
variable $\theta(t)$ obtained from an inverse engineering technique based on
an optimized third-order polynomial (red solid line) and from the optimal
control theory formalism associated to a mean energy minimization (blue dashed
line). Figure 14: Spin dynamics in the presence of dissipation: Time evolution
of the spin components $S_{z}(t)$ (red solid line), $S_{x}(t)$ (blue dashed
line), and $S_{y}(t)$ (black dotted line) under the magnetic field obtained
from the inverse engineering method. Same parameters as Fig. 13. The inset
depicts the corresponding spin trajectory on the Bloch sphere.

However, our simple polynomial ansatz provides an upper bound on the reachable
values of $r_{f}$. This point is illustrated in Fig. 12 (a) where we plot the
energy as a function of the logarithm of the final radius $a_{f}$ for
different values of the free parameter $a_{1}$. For this example, the
reachable range of values for $r_{f}$ is $[0.055;0.476]$. As a result, a
target such as $r_{f}=0.6$ turns out to be out of reach. It is worth noticing
that this limit is intimately related to the choice of the ansatz. For
instance, we can choose a third-order polynomial ansatz:

$\theta(t)=a_{0}+a_{1}t+a_{2}t^{2}+a_{3}t^{3},$ (74)

where $t_{f}$ is determined as previously ($t_{f}=3.6357955$ for $r_{f}=0.6$)
and, the coefficients $a_{0}=0$ and
$a_{2}=-(a_{1}t_{f}+a_{3}{t_{f}}^{3}-\theta_{f})/{t_{f}}^{2}$ are dictated by
the boundary conditions (61). The extra freedom provided by the $a_{3}$
coefficient enables one to (1) reach the target and (2) minimize the cost
function. With the values $a_{3}=0.1$ and $a_{1}=0.15713222$, the cost
function, $E^{(P3)}$, is quite close to the optimal value:
$E^{(P3)}=1.03E^{(OCT)}_{\pi/2}$. In Fig. 12 (b), we plot the energy as a
function of the free parameter $a_{1}$ for $a_{3}=0.1$. This curve defines a
new interval of reachable $r_{f}$: $[0.043;0.608]$. The variable $\theta(t)$
and its corresponding magnetic field $B(t)$ obtained from the latter IE method
are depicted in Fig. 13, and the associated spin trajectory on the Bloch
sphere along with the spin components in Fig. 14. Our results can be a priori
further improved using an optimization on an even larger order polynomial.

###  IV.3 Case II: spin flip

Figure 15:  Spin dynamics in the presence of dissipation: (a) The magnetic
field $B(t)$ and (b) the variable $\theta(t)$ as a function of time for a
minimal-energy spin flip. An optimal ninth-order polynomial has been used for
$\theta(t)$ to apply the inverse engineering method (red solid line). The
optimal solution is plotted as a blue dashed line.

Figure 16:  Spin dynamics in the presence of dissipation: In the case of spin
flip (b) obtained with magnetic field (a), compared with OCT (blue dashed
line), an $\tanh$ ansatz (instead of a polynomial) used in IE approach (red
solid line) is chosen to reduce energy to $E=4.028$, with “magic” values
$a_{5}=1.1$ and $a_{1}=3.104678$.

In this section, we consider a spin flip ($\theta_{f}=\pi$) for which the
constant $p_{1}$ parameter is $p_{1}^{\pi}=2r_{f}/(1-r_{f})^{2}$. With the
same notations as previously, the final time reads (we use $r_{f}=0.6$ in the
following)

$t_{f}^{\pi}=\frac{1-r_{f}}{1+r_{f}}\left[\ln\left(\frac{(1+r_{f})^{2}}{r_{f}}\right)-2\ln\epsilon\right]=3.8165858.$
(75)

The cost function associated to the optimal solution (see Eq. (65)) is

${E}^{(OCT)}_{\pi}=\frac{1+r_{f}}{1-r_{f}}=4.0.$ (76)

This optimal solution is plotted as a blue line in Fig. 15. The optimal
solution exhibits a smooth variations of $\theta(t)$ at initial and final and
a symmetry about $t_{f}/2$. This suggest to add the following extra condition
to the polynomial ansatz for $\theta(t)$ for the inverse engineered solution:

$\displaystyle\theta(0)$ $\displaystyle=$ $\displaystyle
0,\;\;\theta(t_{f}/2)=\pi/2,\;\;\theta(t_{f})=\pi,$
$\displaystyle\dot{\theta}(0)$ $\displaystyle=$
$\displaystyle\dot{\theta}(t_{f})=0,\;\;\mbox{and}\;\;\ddot{\theta}(0)=\ddot{\theta}(t_{f})=0.$
(77)

We have used a ninth-order polynomial to accommodate for the 7 boundary
conditions listed above, an extra parameter is fixed by the final target
radius, $r_{f}$. The remaining two free parameters are used to minimize the
energy. Knowing $\theta(t)$, we infer the magnetic field to be applied to
drive the spin in accordance with our boundary conditions. As explicitly shown
in Fig. 15, we find a bell shape for the magnetic field $B(t)$ associated to
this $\theta(t)$. However, the curves remain relatively far from the optimal
result. We find $E^{(P9)}=1.13E^{(OCT)}$. The ripples in the polynomial ansatz
increase the energy and are difficult to remove by increasing the polynomial
order. The convergence towards the optimal solution is therefore once again
slow with the polynomial order.

Alternatively, the shape obtained from OCT suggests that the following ansatz
could be worth trying:

$\theta(t)=\frac{\pi}{2}\tanh\left\\{a_{1}\tan\left[\frac{\pi}{a_{5}t_{f}}(t-\frac{t_{f}}{2})\right]\right\\}+\frac{\pi}{2}.$
(78)

Minimizing the energy, we find $E=1.007E^{(OCT)}$ with $a_{5}=1.1$ and
$a_{1}=3.104678$. The comparison of this solution with its optimal counterpart
confirms the proximity between the two approaches (see Fig. 16).

##  V Conclusion

In summary, we have investigated different implementations of the inverse
engineering method and compare them with solutions deduced from the OCT for a
given cost function. We have addressed in this manner the fast atomic cooling
in harmonic trap, the atomic transport with a moving harmonic trap, and the
spin control in the presence of dissipation. We have shown how the freedom on
the ansatz inherent to inverse engineering techniques provide enough
tunability to minimize a cost function while fulfilling the boundary
conditions. We have systematically found class of functions with few
adjustable parameters approaching the optimal control result with a relative
excess of energy below one percent. Inverse engineered solutions are usually
search as continuous and analytical functions which is a priori an asset for
their practical use. However, we have also exhibit the possibility to design
inverse engineered trajectories having initial and final jump to mimic the
optimal control solution yielding solutions that are nearly undistinguishable
from their optimal counterpart.

###### Acknowledgements.

This research was funded by NSFC (12075145), STCSM (2019SHZDZX01-ZX04,
18010500400 and 18ZR1415500), Program for Eastern Scholar, Spanish Government
PGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE), Basque Government IT986-16. X. C.
acknowledges Ramón y Cajal program (RYC-2017-22482).

## References

  * (1) E. Torrontegui, S. Ibánez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, in Advances in atomic, molecular, and optical physics (Elsevier, 2013), pp. 117.
  * (2) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Reviews of Modern Physics 91, 045001 (2019).
  * (3) D. Guéry-Odelin, J. Muga, M. J. Ruiz-Montero, and E. Trizac, Physical review letters 112, 180602 (2014).
  * (4) I. A. Martínez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, and S. Ciliberto, Nature physics 12, 843 (2016).
  * (5) H.-C. Chung, S. Martínez-Garaot, X. Chen, J. Muga, and S.-Y. Tseng, EPL (Europhysics Letters) 127, 34001 (2019).
  * (6) S. Faure, S. Ciliberto, E. Trizac, and D. Guéry-Odelin, American Journal of Physics 87, 125 (2019).
  * (7) C. Jarzynski, Physical Review A 88, 040101 (2013).
  * (8) C. Jarzynski, S. Deffner, A. Patra, and Y. Subaşı, Physical Review E 95, 032122 (2017).
  * (9) Y. Li, L.-A. Wu, and Z. Wang, Physical Review A 83, 043804 (2011).
  * (10) S. Masuda and K. Nakamura, Physical Review A 78, 062108 (2008).
  * (11) S. Masuda and K. Nakamura, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 466, 1135 (2010).
  * (12) M. Demirplak and S. A. Rice, The Journal of Physical Chemistry A 107, 9937 (2003).
  * (13) M. Demirplak and S. A. Rice, The Journal of Physical Chemistry B 109, 6838 (2005).
  * (14) M. V. Berry, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009).
  * (15) X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. Muga, Physical review letters 105, 123003 (2010).
  * (16) A. del Campo, Physical review letters 111, 100502 (2013).
  * (17) S. Deffner, C. Jarzynski, and A. del Campo, Physical Review X 4, 021013 (2014).
  * (18) X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Physical review letters 104, 063002 (2010).
  * (19) X. Chen, E. Torrontegui, and J. G. Muga, Physical Review A 83, 062116 (2011).
  * (20) E. Torrontegui, S. Martínez-Garaot, A. Ruschhaupt, and J. G. Muga, Physical Review A 86, 013601 (2012).
  * (21) S. Ibáñez, X. Chen, E. Torrontegui, J. G. Muga, and A. Ruschhaupt, Physical review letters 109, 100403 (2012).
  * (22) Y.-X. Du, Z.-T. Liang, Y.-C. Li, X.-X. Yue, Q.-X. Lv, W. Huang, X. Chen, H. Yan, and S.-L. Zhu, Nature communications 7, 1 (2016).
  * (23) J.-F. Schaff, X.-L. Song, P. Vignolo, and G. Labeyrie, Physical Review A 82, 033430 (2010).
  * (24) J.-F. Schaff, X.-L. Song, P. Capuzzi, P. Vignolo, and G. Labeyrie, EPL (Europhysics Letters) 93, 23001 (2011).
  * (25) J. Zhang et al., Physical review letters 110, 240501 (2013).
  * (26) B. B. Zhou et al., Nature Physics 13, 330 (2017).
  * (27) C. Munuera-Javaloy, Y. Ban, X. Chen, and J. Casanova, Physical Review Applied 14, 054054 (2020).
  * (28) S. An, D. Lv, A. Del Campo, and K. Kim, Nature communications 7, 12999 (2016).
  * (29) T. Wang et al., New Journal of Physics 20, 065003 (2018).
  * (30) T. Wang et al., Physical Review Applied 11, 034030 (2019).
  * (31) A. Vepsäläinen, S. Danilin, and G. S. Paraoanu, Science advances 5, eaau5999 (2019).
  * (32) T. Yan et al., Physical review letters 122, 080501 (2019).
  * (33) S. J. Glaser et al., The European Physical Journal D 69, 1 (2015).
  * (34) E. Assémat, L. Attar, M.-J. Penouilh, M. Picquet, A. Tabard, Y. Zhang, S. Glaser, and D. Sugny, Chemical Physics 405, 71 (2012).
  * (35) B. M. Henson, D. K. Shin, K. F. Thomas, J. A. Ross, M. R. Hush, S. S. Hodgman, and A. G. Truscott, Proceedings of the National Academy of Sciences 115, 13216 (2018).
  * (36) J. J. W. Sørensen et al., Nature 532, 210 (2016).
  * (37) D. Sels, Physical Review A 97, 040302 (2018).
  * (38) Y. Ding, Y. Ban, J. D. Martín-Guerrero, E. Solano, J. Casanova, and X. Chen, arXiv preprint arXiv:2009.04297 (2020).
  * (39) M. Larocca, E. Calzetta, and D. A. Wisniacki, Physical Review A 101, 023410 (2020).
  * (40) D. Stefanatos, J. Ruths, and J.-S. Li, Physical Review A 82, 063422 (2010).
  * (41) D. Stefanatos and J.-S. Li, Physical Review A 86, 063602 (2012).
  * (42) X.-J. Lu, X. Chen, J. Alonso, and J. Muga, Physical Review A 89, 023627 (2014).
  * (43) X. Chen, E. Torrontegui, D. Stefanatos, J.-S. Li, and J. Muga, Physical Review A 84, 043415 (2011).
  * (44) Q. Zhang, X. Chen, and D. Guéry-Odelin, Physical Review A 92, 043410 (2015).
  * (45) Q. Zhang, J. G. Muga, D. Guéry-Odelin, and X. Chen, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 125503 (2016).
  * (46) A. Ruschhaupt, X. Chen, D. Alonso, and J. Muga, New Journal of Physics 14, 093040 (2012).
  * (47) X.-J. Lu, X. Chen, A. Ruschhaupt, D. Alonso, S. Guerin, and J. G. Muga, Physical Review A 88, 033406 (2013).
  * (48) D. Daems, A. Ruschhaupt, D. Sugny, and S. Guerin, Physical Review Letters 111, 050404 (2013).
  * (49) V. Martikyan, D. Guéry-Odelin, and D. Sugny, Physical Review A 101, 013423 (2020).
  * (50) D. Guéry-Odelin and J. Muga, Physical Review A 90, 063425 (2014).
  * (51) V. Martikyan, A. Devra, D. Guéry-Odelin, S. Glaser, and D. Sugny, Physical Review A 102, 053104 (2020).
  * (52) S. Choi, R. Onofrio, and B. Sundaram, Physical Review A 84, 051601 (2011).
  * (53) A. Del Campo, Physical Review A 84, 031606 (2011).
  * (54) S. Deng, P. Diao, Q. Yu, A. del Campo, and H. Wu, Physical Review A 97, 013628 (2018).
  * (55) W. Rohringer, D. Fischer, F. Steiner, I. E. Mazets, J. Schmiedmayer, and M. Trupke, Scientific reports 5, 9820 (2015).
  * (56) X. Chen and J. G. Muga, Physical Review A 82, 053403 (2010).
  * (57) K. Hoffmann, P. Salamon, Y. Rezek, and R. Kosloff, EPL (Europhysics Letters) 96, 60015 (2011).
  * (58) Y.-Y. Cui, X. Chen, and J. Muga, The Journal of Physical Chemistry A 120, 2962 (2016).
  * (59) A. Couvert, T. Kawalec, G. Reinaudi, and D. Guéry-Odelin, EPL (Europhysics Letters) 83, 13001 (2008).
  * (60) A. Walther, F. Ziesel, T. Ruster, S. T. Dawkins, K. Ott, M. Hettrich, K. Singer, F. Schmidt-Kaler, and U. Poschinger, Physical Review Letters 109, 080501 (2012).
  * (61) R. Bowler, J. Gaebler, Y. Lin, T. R. Tan, D. Hanneke, J. D. Jost, J. Home, D. Leibfried, and D. J. Wineland, Physical review letters 109, 080502 (2012).
  * (62) E. Torrontegui, S. Ibáñez, X. Chen, A. Ruschhaupt, D. Guéry-Odelin, and J. Muga, Physical Review A 83, 013415 (2011).
  * (63) E. Torrontegui, X. Chen, M. Modugno, S. Schmidt, A. Ruschhaupt, and J. Muga, New Journal of Physics 14, 013031 (2012).
  * (64) R. Corgier, S. Amri, W. Herr, H. Ahlers, J. Rudolph, D. Guéry-Odelin, E. M. Rasel, E. Charron, and N. Gaaloul, New Journal of Physics 20, 055002 (2018).
  * (65) D. Becker et al., Nature 562, 391 (2018).
  * (66) S. Amri, R. Corgier, D. Sugny, E. M. Rasel, N. Gaaloul, and E. Charron, Scientific Reports 9, 5346 (2019).
  * (67) M. Dupont-Nivet, C. Westbrook, and S. Schwartz, New Journal of Physics 18, 113012 (2016).
  * (68) X.-J. Lu, J. G. Muga, X. Chen, U. Poschinger, F. Schmidt-Kaler, and A. Ruschhaupt, Physical Review A 89, 063414 (2014).
  * (69) X.-J. Lu, A. Ruschhaupt, and J. G. Muga, Physical Review A 97, 053402 (2018).
  * (70) U. Boscain, G. Charlot, J.-P. Gauthier, S. Guérin, and H.-R. Jauslin, Journal of Mathematical Physics 43, 2107 (2002).
  * (71) L. Van Damme, Q. Ansel, S. Glaser, and D. Sugny, Physical Review A 95, 063403 (2017).
  * (72) L. Van Damme, Q. Ansel, S. Glaser, and D. Sugny, Physical Review A 98, 043421 (2018).
  * (73) D. Stefanatos, Physical Review A 80, 045401 (2009).
  * (74) D. Stefanatos and E. Paspalakis, Physical Review A 99, 022327 (2019).
  * (75) P. Owrutsky and N. Khaneja, Physical Review A 86, 022315 (2012).
  * (76) C. Lin, D. Sels, and Y. Wang, Physical Review A 101, 022320 (2020).
  * (77) H. L. Mortensen, J. J. W. Sørensen, K. Mølmer, and J. F. Sherson, New Journal of Physics 20, 025009 (2018).
  * (78) Q. Zhang, X. Chen, and D. Guéry-Odelin, Scientific Reports 7, 15814 (2017).
  * (79) X.-T. Yu, Q. Zhang, Y. Ban, and X. Chen, Physical Review A 97, 062317 (2018).
  * (80) F. Impens and D. Guéry-Odelin, Scientific Reports 9, 4048 (2019).


* [96] Desroches M, Krauskopf B, Osinga HM. 2010 Numerical continuation of canard orbits in slow–fast dynamical systems. Nonlinearity 23, 739–765.
  * [97] Desroches M, Krauskopf B, Osinga HM. 2008 The geometry of slow manifolds near a folded node. SIAM Journal on Applied Dynamical Systems 7, 1131–1162.
  * [98] Perryman CG. 2015 How fast is too fast? Rate-induced bifurcations in multiple time-scale systems. PhD thesis University of Exeter.
  * [99] Marchionne A, Ditlevsen P, Wieczorek S. 2018 Synchronisation vs. resonance: Isolated resonances in damped nonlinear oscillators. Physica D: Nonlinear Phenomena 380, 8–16.
  * [100] Yuan H, Restuccia F, Rein G. 2021 Spontaneous ignition of soils: A multi-step reaction scheme to simulate self-heating ignition of smouldering peat fires. International Journal of Wildland Fire 30, 440.
  * [101] Sierra CA, Trumbore SE, Davidson EA, Vicca S, Janssens I. 2015 Sensitivity of decomposition rates of soil organic matter with respect to simultaneous changes in temperature and moisture. Journal of Advances in Modeling Earth Systems 7, 335–356.
  * [102] 2022 Resources: water@leeds. See https://water.leeds.ac.uk/our-missions/mission-1/digibog/resources/. Accessed: 2022-09-12.
  * [103] Waddington JM, Morris PJ, Kettridge N, Granath G, Thompson DK, Moore PA. 2014 Hydrological feedbacks in northern peatlands. Ecohydrology 8, 113–127.
  * [104] Ritchie P, Sieber J. 2017 Probability of noise- and rate-induced tipping. Phys. Rev. E 95, 052209.
  * [105] Slyman K, Jones CKRT. 2022 Rate and Noise Tipping Working in Concert. In Tipping Phenomena: Nonlinear Theory and Climate Applications. SIAM Conference on Mathematics of Planet Earth (MPE22).
  * [106] 2022 Joint UK Land Environment Simulator (JULES). See https://jules.jchmr.org/. Accessed: 2022-09-12.
  * [107] 2022 Digibog: water@leeds. See https://water.leeds.ac.uk/our-missions/mission-1/digibog/. Accessed: 2022-09-12.