# Robust control of linear systems and shortcut to adiabaticity based on
superoscillations

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

###### Abstract

With the advent of quantum technologies, control issues are becoming
increasingly important. In this article, we address the control in phase space
under a global constraint provided by a minimal energy-like cost function and
a local (in Fourier space) constraint resulting from a robustness criterion.
The inverse engineering Lagrangian formalism developed here generalizes the
one commonly used to describe the superoscillation phenomenon. It is applied
to both non-dissipative and dissipative quantum mechanics, and extended to
stochastic thermodynamics. Interestingly, our approach also allows to improve
the sensing capabilities by an appropriate control of the system.

Quantum control is nowadays a crucial and implicit step in the manipulation of
a wide variety of quantum systems, including, but not limited to, the quantum
computation in the noisy intermediate-scale quantum (NISQ) era EPJD ;
review2022 . Many methods have been developed in this field ranging from
adiabatic control to optimal control theory EPJD ; review2022  and shortcuts
to adiabaticity STARMP . For instance, protocols for high-fidelity fast
transport or shuttling of charged ions Walther12 ; Bowler12 ; Qi2021  and
neutral atoms DGO08 ; davidmugapra  have been worked out, and are essential
for the development of quantum technologies applications including quantum
information processing Erikion ; ErikBEC , atom chip manipulation CorgierNJP ;
Becker ; Amri  and metrology with cold atoms Dupont .

Besides reaching the desired target state, one cannot discard the question of
the robustness of the designed protocols against various sources of noise, an
inaccurate knowledge of the experimental parameters or model approximations
Xipra11optimal ; Xiaojingpra ; an2015 ; QiJPB ; Xiaopra18 ; EntropyDGO ;
Muga2022 . This concerns is particularly relevant to ions or neutral-atoms
transport based on moving traps. Incidentally, such systems share the same
dynamical equations as the load manipulation of mechanical cranes CranePRAppl
; EnergyPRA , and the motion of a charged particle in an external electric
field OTCSTA . These quantum-classical analogies allow to use control
protocols in ubiquitous linear Lagrangian systems and compare the methods
OTCSTA ; controldavidsugny . Furthermore, the above strategies can be applied
to linear dissipative systems. This is illustrated with the motion of
harmonically trapped Brownian particles driven by a time-varying force
Cunuder16 ; TuPRE ; JunPRR ; New18 .

In this article, we address the problem of the optimal control in a one-
dimensional (1D) phase space for a given protocol duration and under a global
constraint provided by an energy cost and a robustness criterion that can be
formulated as a local constraint in Fourier space. Such a control problem
encompasses the transport and shuttling (see STARMP ; Qi2021  and references
therein) problems that have been addressed separately in the literature so
far. The focus of the article is on the resolution of linear systems under
such constraints, in the context of inverse engineering Lagrangian formalism.
However, solving such a classical problem also provides the solution for the
corresponding quantum problem davidmugapra . We exemplify this idea with an
application of our formalism to quantum transport described by coherent state
obeying a Lindblad equation.

In the following, we first provide a short reminder on how optimal control
theory is applied to linear systems. We then identify an intrinsic robustness
against slight final time uncertainties. Later, we translate mathematically
the requirement of robustness against an inaccurate knowledge of the strength
of the moving transport potential, and explain how this extra constraint can
be encapsulated in an inverse engineering Lagrange formulation. Finally, the
detailed applications are exemplified, ranging from classical and quantum
mechanics to stochastic thermodynamics.

Consider a linear systems described by the following set of linear equations
$\dot{\bm{x}}=\bm{A}\bm{x}(t)+\bm{B}\bm{u}(t)$, where
$\bm{x}(t)\in\mathbb{R}^{n}$ denotes the state vector and
$\bm{u}(t)\in\mathbb{R}^{q}$ the control vector. The dynamics is assumed here
to involve two time independent matrices: $\bm{A}\in\mathbb{R}^{n\times n}$
and $\bm{B}\in\mathbb{R}^{n\times q}$. We assume the system to be
controllable, i.e. it fulfills the Kalman’s controllability criterion stating
that the rank of the Kalman’s controllability matrix
$\bm{K}_{c}=[\bm{B}\;\bm{A}\bm{B}\;\bm{A}^{2}\bm{B}\;...\;\bm{A}^{n-1}\bm{B}]$
is $n$, see controldavidsugny . The solution then reads
$\bm{x}(t)=e^{\bm{A}t}\bm{x}_{0}+\bm{I}_{u}(t)$ with
$\bm{I}_{u}(t)=\int_{0}^{t}e^{\bm{A}(t-t^{\prime})}\bm{B}\bm{u}(t^{\prime})dt^{\prime}$.
The control problem consists in designing the time-dependent control vector to
drive the system from an initial state ($t_{0},\bm{x}_{0}$) to the desired
target state ($t_{f},\bm{x}_{f}$) in a time interval that lasts $t_{f}$. The
boundary condition at final time imposes
$\bm{I}_{u}(t_{f})=\bm{x}_{f}-e^{\bm{A}t_{f}}\bm{x}_{0}$. In most physical
situations, extra requirements are needed such as the minimization of a cost
function $E=\int^{t_{f}}_{0}f[\bm{x}(t),\bm{u}(t),t]dt$ that depends on both
the trajectory $\bm{x}(t)$ and the driving $\bm{u}$. such a minimization
problem can be solved thanks to the Lagrangian multiplier formalism:

$\frac{\partial
E}{\partial\bm{u}}-\bm{\mu}_{0}^{T}\frac{\partial\bm{G}_{0}(\bm{u})}{\partial\bm{u}}=0,$
(1)

where
$\bm{G}_{0}(\bm{u})=\bm{I}_{u}(t_{f})-\bm{x}_{f}+e^{\bm{A}t_{f}}\bm{x}_{0}$.
For the cost function, $E_{1}=\int_{0}^{t_{f}}\bm{u}^{T}\bm{u}dt$ the
Lagrangian multiplier is given by applying the controllability Gramian matrix
$\bm{\mu}_{0}=\bm{W}(t_{f})^{-1}\big{(}\bm{x}_{f}-e^{\bm{A}t_{f}}\bm{x}_{0}\big{)}$
Gramian . Alternatively, the very same results can be recovered using the
Pontryagin formalism of optimal control theory.

For a particle in a moving harmonic trap characterized by its time-dependent
center at $x_{0}(t)$ and its constant angular frequency $\omega_{0}$, the
dynamical equation reads
$\ddot{x}+\omega_{0}^{2}x=\omega_{0}^{2}x_{0}(t)\equiv u(t)$, where $u(t)$ is
the control parameter that drives the system. This second-order differential
equation is readily transformed into a set of two first order differential
equation with the variable $x$ and $y=\dot{x}$. The control problem we focus
on consists in designing the time-dependent control parameter $u(t)$ to
transfer the system from an initial condition $(x_{i}=0,~{}y_{i}=0)$ to an
arbitrary target $(x_{f}=r\cos\varphi,~{}y_{f}=r\omega_{0}\sin\varphi)$ in
phase space in a finite amount of time $t_{f}$. For $y_{f}\neq 0$, such a
control problem is referred to as shuttling Erikion ; Qi2021 . The particle
acquires the desired velocity at a given position after application of a
proper driving. For a vanishing final velocity $y_{f}=0$, this problem is
nothing but a transport problem for which the particle initially at rest is
displaced from $x=0$ to $x=x_{f}$ and remains at rest afterwards DGO08 ;
Erikion .

Figure 1: The energy cost $E_{1}$ as a function of the protocol duration
$t_{f}$ in units of $\omega_{0}^{-1}$ for transport ($\varphi=0$, solid black
line) and shuttling ($\varphi=\pi/2$, solid gray line). Parameters: $d=1$,
$\omega_{0}t_{f}=10$.

Using the previous formalism for the minimization of the standard cost
function $E_{1}=\int_{0}^{t_{f}}u^{2}dt$, the 2-dimensional Lagrangian
multiplier $\bm{\mu}_{0}$ can be readily derived from the boundary conditions
yielding the following expression for the control parameter

$u_{1}(t)\\!=\\!\frac{2r\omega_{0}^{2}[\varphi_{f}\sin(\varphi\\!+\\!\varphi_{f}\\!-\\!\omega_{0}t)\\!-\\!\sin(\omega_{0}t\\!+\\!\varphi)\sin(\varphi_{f})]}{\varphi_{f}^{2}\\!-\\!\sin^{2}(\varphi_{f})},$

with $\varphi_{f}=\omega_{0}t_{f}$. The corresponding minimal energy cost
$E_{1}$ is

$E_{1}=2r^{2}\omega_{0}^{3}\frac{\varphi_{f}+\cos(2\varphi+\varphi_{f})\sin(\varphi_{f})}{\varphi_{f}^{2}-\sin^{2}(\varphi_{f})}.$
(2)

This energy cost provides the lower bound to reach the target
($r\cos\varphi,\,r\omega_{0}\sin\varphi$) in an amount of time $t_{f}$. In
Fig. 1 (a), we have depicted this energy as a function of the final time for
both a transport problem (target with $\varphi=0$, black solid line) and a
shuttling process (target with $\varphi=\pi/2$, gray solid line). As
intuitively expected, the shorter the driving time $t_{f}$, the larger the
energetic cost $E_{1}$, and $E_{1}$ is a decreasing function of
$\omega_{0}t_{f}$. Interestingly, this result reveals the existence of
plateaus for which the change of $t_{f}$ has a negligible impact on the energy
cost, meaning that the corresponding solution are intrinsically robust against
a small variation of the protocol duration $t_{f}$ for a succession of windows
of final time. The energetic cost (2) is the lowest possible as it is derived
without any constraints. We denote by $E_{1}^{\rm(min)}$ this lower bound.

Next, we add extra constraints and derive an inverse engineering Lagrangian
approach to build protocols robust against the exact knowledge of the trap
frequency $\omega_{0}$. The energetic cost $E_{1}$ associated to such a
requirement will be quantified and compared to $E_{1}^{\rm(min)}$. As a first
example, we elaborate the transport problem, i.e. a target $r=d$ and
$\varphi=0$ to be reached in an amount of time $t_{f}$. For an arbitrary
driving $x_{0}(t)=u(t)/\omega_{0}$, the final excitation energy reads
davidmugapra

$\Delta
E(\omega_{0};t_{f})=d^{2}\omega_{0}^{2}\left|\int_{0}^{t_{f}}\ddot{x}_{0}(t)e^{-i\omega_{0}t}dt\right|^{2}\equiv\left|g(\omega_{0};t_{f})\right|^{2}.$
(3)

A perfect transport (without residual ossilations) requires $\Delta
E(\omega_{0};t_{f})=0$. A transport robust against the exact value of
$\omega_{0}$ requires to smoothen $\Delta E$ as a function of $\omega$ about
$\omega_{0}$ davidmugapra :

$\frac{\partial^{p}\Delta
E(\omega;t_{f})}{\partial\omega^{p}}\bigg{|}_{\omega_{0}}=0,\;\;\mbox{with}\;\;p\geq
1.$ (4)

The calculation associated to this robustness issue are made simpler by using
the complex variable $z(t)=\dot{x}/\omega_{0}+ix$. The equations of motion are
then $\dot{z}(t)=i\omega_{0}z(t)+u(t)/\omega_{0}$. Under the driving $u(t)$,
the complex solution is
$z(t)=e^{i\omega_{0}t}z(0)+\tilde{u}(\omega_{0})/\omega_{0}$ where
$\tilde{u}(\omega_{0})=\int^{t_{f}}_{0}u(t)e^{-i\omega_{0}t}dt\equiv a_{0}$.
For an arbitrary target $z(t_{f})=r(\sin\varphi+i\cos\varphi)$,
$a_{0}=r\omega_{0}\left(\sin\varphi+i\cos\varphi\right)e^{-i\omega_{0}t_{f}}$.

For transport, the driving function obeys the boundary conditions: $u(0)=0$,
$u(t_{f})=\omega_{0}d$, and $\dot{u}(0)=\dot{u}(t_{f})=0$. To fulfill the
robustness requirements Eq. (4), we impose the nullity of the successive
derivative of $g(\omega;t_{f})$ (defined in Eq. (3)) with respect to $\omega$
about $\omega_{0}$. We infer the following set of extra constraints on the
Fourier function $\tilde{u}(\omega)$ and its successive derivatives:

$\omega_{0}\tilde{u}^{(p)}(\omega_{0})+2p\tilde{u}^{(p-1)}(\omega_{0})+\frac{p(p-1)}{\omega_{0}}\tilde{u}^{(p-2)}(\omega_{0})=\lambda,$
(5)

where $\lambda=d(\varphi_{f}+ip)(-it_{f})^{p-1}e^{-i\varphi_{f}}$. We find
$\tilde{u}^{(m)}(\omega_{0})=a_{m}$ with
$a_{1}=d(\varphi_{f}-i)e^{-i\varphi_{f}}/\omega_{0}$,
$a_{2}=d(-i\varphi_{f}^{2}-2\varphi_{f}+2i)e^{-i\varphi_{f}}/\omega_{0}^{2}$
and so on.

More generally, the class of problem to be solved here requires the
minimization of $E=\int_{\cal D}f(t)f^{*}(t)\mathrm{d}t$ on a compact support
${\cal D}$ under a discrete number of constraints on the Fourier transform of
$f$ and its derivatives for a finite set of angular frequencies $\omega_{j}$:
$\tilde{f}^{(m)}(\omega_{j})=\int_{\cal
D}(-it)^{m}f(t)e^{-i\omega_{j}t}\mathrm{d}t=\tilde{f}^{(m)}_{j}$. For a single
constraint on the Fourier transform, such a problem has been solved in the
context of superoscillations detailed in aharonove ; berry1 ; berry2 ;
SuperOsc . We generalize here this approach to enforce the robustness
requirement. The solution relies on the Euler-Lagrange formalism:

$\frac{\delta E}{\delta
f}+\sum_{j}\sum_{m=0}^{p}\nu_{j,m}\frac{\delta\tilde{f}^{(m)}(\omega_{j})}{\delta
f}=0,$ (6)

where $\nu_{j,m}$ are Lagrange multipliers. Interestingly, this set of
equations can be directly solved by a simple polynomial ansatz

$f(t)=\sum_{j}\bigg{(}\nu_{j0}^{*}+\sum_{m}(it)^{m}~{}\nu_{j,m}^{*}\bigg{)}e^{i\omega_{j}t}.$
(7)

The Lagrange multipliers, $\nu_{j,m}$, are subsequently determined by the
constraints. The inverse engineering procedure amounts here to fulfilling the
boundaries condition to set up the driving. One should be careful to keep
$u(t)$ real as it coincides with the position of the trap in real space. To
get a real driving function, we simply exploit the linearity and define the
control function as $u(t)=(f(t)+f^{*}(t))/2$.

Figure 2: In the case of transport ($\varphi=0$), (a) trajectories in phase
space, and (b) the final excitation energy as a function of $\omega$ are
compared between (i) Fourier Euler-Lagrange method with no extra robust
constraint protocol (black solid line), $p=1$ robust control (blue dot-dashed
line) and $p=2$ robust one (red double-dot-dashed line) and (ii) and Fourier
method for $N=1$ one frequency robust protocol (gray dot-dashed thin line) and
$N=2$ two frequencies robust one (gray double-dot-dashed thin line).
Parameters: $d=1$, $\omega_{0}t_{f}=10$.

In Fig. 2, we plot the excess of energy after transport in a harmonic
potential as a function of its angular frequency $\omega$ for protocols
optimized for a given $\omega_{0}$ and for a transport duration
$\omega_{0}t_{f}=10$. We have computed the solutions under the robustness
requirements on the first (blue dot-dashed line) and second (red double-dot-
dashed line) derivative of the Fourier transform of the driving function. We
clearly observe a strong improvement associated to the local flatness about
$\omega_{0}$ compared to the solution obtained by optimal control without
constraints (black solid line). We have also compared the solutions, see the
gray lines in Fig. 2, obtained by the so-called Fourier method davidmugapra
that has been developed to address similar requirements.

In fact, the formalism developed here is well-adapted to the optimization of
other class of transport problems, since solving this classical problem
amounts to solving its quantum counterpart as explained in Ref. davidmugapra .
The transport problem in the presence of dissipation requires to solve the
equation of motion for the density matrix. Regarding transport with losses,
the dynamics of coherent state $\left.|\alpha\right\rangle$ is given by the
following master equation:

$\frac{d\rho}{dt}=-\frac{i}{\hbar}[H,\rho]+\frac{\Gamma}{2}(2a\rho
a^{{\dagger}}-a^{{\dagger}}a\rho-\rho a^{{\dagger}}a),$ (8)

where $\Gamma$ is the Lindblad operator and $H$ the Hamiltonian,
$H=p^{2}/2m+m\omega^{2}[\hat{x}-x_{0}(t)]^{2}/2$, with $x_{0}(t)$ being the
center of moving potential Erikion ; DGO08 . This problem can be solved
exactly whatever is the driving using a time-dependent coherent state
$\rho(t)=\left.|\alpha\right\rangle\left\langle\alpha\right|$ with
$\alpha(t)=(i\omega/\sqrt{2}a_{0})e^{-(i\omega+\Gamma/2)t}\int_{0}^{t}e^{(i\omega+\Gamma/2)t^{\prime}}x_{0}(t^{\prime})dt^{\prime}$.
This problem belongs to the same class of problems as previously studied since
it can be expressed as a linear system on the variables
$x_{1}=\left\langle\hat{x}\right\rangle$,
$x_{2}=d\left\langle\hat{x}\right\rangle/dt$: $\dot{x}_{1}=x_{2}$ and
$\dot{x}_{2}=-(\omega_{0}^{2}+\Gamma^{2}/4)x_{1}-\Gamma
x_{2}+\omega_{0}^{2}x_{0}(t)$.

Figure 3:  Shuttling protocols (target $\varphi=\pi/2$). (a) Phase space
trajectories associated to the different protocols performed over a
$t_{f}=10/\omega_{0}$ time interval. (b) Energy cost of the different
protocols as a function of the protocol duration. Same notations as Fig. 2.

Moreover, the transport of a mesoscopic object in stochastic thermodynamics
can also be captured and optimized by our formalism. For this purpose, it
turns out to be more convenient to work with the distribution function rather
than with the Langevin-like stochastic differential equations. In the
overdamped regime, the distribution function $\rho(x,t)$ obeys the following
Fokker-Planck equation:

$\gamma\partial_{t}\rho(x,t)=\partial_{x}[m\omega_{0}^{2}\big{(}x-x_{0}(t)\big{)}\rho]+\gamma
D\partial_{xx}^{2}\rho,$ (9)

where $\gamma$ is the friction coefficient and $D=k_{B}T/\gamma$ the diffusion
coefficient. The solution of the FP equation is provided by a Gaussian
$\rho(x,t)=\sqrt{\frac{\alpha}{\pi}}\exp[-\alpha(x-x_{c}(t))^{2}]$ with
$\alpha=2\gamma D/(m\omega_{0}^{2})$. The time-dependence of $\rho(x,t)$ is
encapsulated in the variable $x_{c}(t)$ governed by the linear system
$\dot{x}_{c}+(m\omega_{0}^{2}/\gamma)(x_{c}-x_{0})=0$ for which the previous
formalism applies. These results can be readily generalized to the underdamped
regime. For this purpose, we introduce the phase space distribution
$\rho(x,v,t)$ that obeys the Kramers equations Cunuder16 . A Gaussian solution
can be readily worked out
$\rho(x,v,t)=\mathcal{N}\exp[-\alpha(x-x_{c}(t))^{2}-\beta(v-v_{c}(t))^{2}]$
with $\alpha=m\omega_{0}^{2}/2k_{B}T$ and $\beta=m/2k_{B}T\gamma$, and with a
time variation for the parameters $x_{c}$ and $v_{c}$ dictated by the
following set of linear equations: $\dot{x}_{c}=v_{c}$,
$\dot{v}_{c}=-\omega_{0}^{2}x_{c}-\gamma v_{c}/m+\omega_{0}^{2}x_{0}(t)$.

Besides targets in phase space located on the position axis, our protocol
encompasses the shuttling problem for which the particle acquire a finite
velocity. The results obtained for such shuttling protocols are summarized in
Fig. 3. As a general result, the inverse engineering Lagrangian method
developed to ensure robustness provides systematically a much lower $E_{1}$
energy for a time duration of the protocol larger than $t_{f}>10/\omega_{0}$
compared to the Fourier method davidmugapra . The latter turns out to be more
efficient for extremely short shuttling protocol duration but at the expense
of a large energy cost $E_{1}$. In Fig. 3, we also provides for comparison the
minimum energy obtained by optimal control theory for the minimization of the
energy cost in the absence of any other constraints (black solid line). We
therefore clearly identify the minimum cost of the robustness requirement.

Figure 4: Sensing the frequency between $\omega_{2}=1.05\omega_{1}$ by
separating the corresponding target states as follows:
$(r_{1}=\nu,\varphi_{1}=\pi/2)$ and $(r_{2}=d,\varphi_{2}=0)$. The comparison
is between the following two protocols: Fourier Euler-Lagrange approach (red
solid line) and Fourier method (blue dashed line). Parameter: $t_{f}=10$ and
$\omega_{1}=1$.

Finally, the inverse engineering Lagrangian formalism enables one to design a
driven protocol for quantum sensing, with the objective to measure the trap
frequency or its drift as a function of time. To this end, a protocol is
designed to assign different targets to two close but different trapping
frequencies $\omega_{1}$ and $\omega_{2}=(1+\epsilon)\omega_{1}$
($\epsilon\neq 0$). Such a strategy is reminiscent of the fingerprinting
method used in NMR Ma13 ; Ansel17 . We therefore consider here the same
control function that drives two independent set of equations:
$\dot{x}_{j}=y_{j}$ and $\dot{y}_{j}=-\omega_{j}^{2}y_{j}+u(t)$. To transfer
the system from the initial state
$\big{(}x_{j}(0),\dot{x}_{j}(0)\big{)}=(0,0)$ to distinguishable target states
$\big{(}x_{j}(t_{f}),\dot{x}_{j}(t_{f})\big{)}=(r_{j}\cos\varphi_{j},r_{j}\omega_{j}\sin\varphi_{j})$.
The Fourier Transform of control $u(t)$ should be equal to
$\tilde{u}(\omega_{j})=\omega_{j}r_{j}\left(\sin\varphi_{j}+i\cos\varphi_{j}\right)e^{-i\omega_{j}t_{f}}$.
Such conditions are readily fulfilled with an ansatz of the form
$u(t)=\sum_{j}[a_{j}\cos(\omega_{j}t)-b_{j}\sin(\omega_{j}t)]$. For example,
the two frequencies case refers to $j=1,2$ and the target can be
$r_{1}=\nu,\varphi_{1}=\pi/2$, $r_{2}=d$ and $\varphi_{2}=0$. Such conditions
are readily fulfilled with an ansatz of the form
$u(t)=\sum_{j=1}^{2}[a_{m}\cos(\omega_{m}t)-b_{m}\sin(\omega_{m}t)]$ . As
shown in Fig. 4, we plot in red the evolution of the trajectories for the two
different $\omega_{i}$ towards their target, along with the final state in
phase space $(x_{f},\dot{x}_{f})$ for different trapping angular frequencies
$\omega\in(\omega_{1},\omega_{2})$. It demonstrates the possibility to magnify
the sensitivity. Indeed, in the absence of driving, the information on the
trapping angular frequency cannot be recovered as the system remains at rest.
This sensing protocol is particularly relevant when the time duration of the
protocol $t_{f}$ remains smaller than the natural timescale
$\pi/(2\varepsilon\omega_{1})$ which would generate an angle separation in
phase space equal to $\pi/4$ for a step driving (sudden switch of the driving
function to a constant value). This is the reason why in Fig. 4, we have
chosen the parameters to fulfill the inequality
$\omega_{1}t_{f}\ll\pi/(2\varepsilon)$. The curves in blue dashed lines have
been obtained using de Fourier method davidmugapra . The same conclusion holds
concerning the sensing protocol. It is worth noticing that the Lagrange method
limits considerably the extension of the trajectory to obtain a similar
sensitivity. This can be a strong asset to limit the role of anharmonicities
in such a protocol an2015 ; QiJPB . The fingerprint sensing method discussed
here can be readily generalized to a larger number of $\omega_{i}$.

To conclude, we have developed a formalism inspired by Lagrangian minimization
to tackle optimal control of linear systems taking into account the different
constraints, such as energy-like minimization and robustness criterion. Our
inverse engineering Lagrangian formalism generalizes the one for driving
quantum systems with superoscillations to enforce the robustness or
sensitivity requirements. This has been adopted to both non-dissipative and
dissipative linear systems, describing the transport of ions or neutral atoms
with/without loss, and a mesoscopic object in stochastic thermodynamics. Our
results pave the way for the design of optimal shortcuts to adiabaticity in a
class of general linear systems with broad applications in quantum
technologies.

This work has been financially supported by NSFC (12075145), STCSM
(2019SHZDZX01-ZX04), the Agence Nationale de la Recherche through grant
ANR-18-CE30-0013, EU FET Open Grant EPIQUS (899368), QUANTEK project
(KK-2021/00070). X.C. acknowledges the Ramón y Cajal program (RYC-2017-22482).

## References

  * (1) S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, Eur. Phys. J. D 69, 279 (2015).
  * (2) C. P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhel, arXiv:2205.12110.
  * (3) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga Rev. Mod. Phys. 91, 045001 (2019).
  * (4) A. Walther, F. Ziesel, T. Ruster, S. T. Dawkins, K. Ott, M. Hettrich, K. Singer, F. Schmidt-Kaler, U. Poschinger, Phys. Rev. Lett. 109, 080501, (2012).
  * (5) R. Bowler, J. Gaebler, Y. Lin, T. R. Tan, D. Hanneke, J. D. Jost, J. P. Home, D. Leibfried, D. Wineland, Phys. Rev. Lett. 109, 080502 (2012).
  * (6) L. Qi, J. Chiaverini, H. Espinós, M. Palmero, and J. G. Muga EPL 134, 23001 (2021).
  * (7) A. Couvert, T. Kawalec, G. Reinaudi, D. Guéry-Odelin, Euro. Phys. Lett. 83, 13001 (2008).
  * (8) D. Guéry-Odelin and J. G. Muga, Phys. Rev. A 90, 063425 (2014).
  * (9) E. Torrontegui, S. Ibáñez, X. Chen, A. Ruschhaupt, D. Guéry-Odelin, J. Muga, Phys. Rev. A 83, 013415 (2011).
  * (10) E. Torrontegui, X. Chen, M. Modugno, S. Schmidt, A. Ruschhaupt, J. G. Muga, New J. Phys.14, 013031 (2012).
  * (11) R. Corgier, S. Amri, W. Herr, H. Ahlers, J. Rudolph, D. Guéry-Odelin, E. M. Rasel, E. Charron, N. Gaaloul, New J. Phys., 20, 055002 (2018).
  * (12) D. Becker, M. D.Lachmann, S. T. Seidel, H. Ahlers, A. N. Dinkelaker, J. Grosse, O. Hellmig, H. Müntinga, V. Schkolnik, T. Wendrich, Nature 562, 391 (2018).
  * (13) S. Amri, R. Corgier, D. Sugny, E. M. Rasel, N. Gaaloul, E. Charron, Scientific Reports, 9, 1, (2019).
  * (14) M. Dupont-Nivet, C. I. Westbrook, S. Schwartz, New J. Phys. 18, 113012, (2016).
  * (15) X. Chen, E. Torrontegui, D. Stefanatos, J.-S. Li, J. G. Muga, J., Phys. Rev. A 84, 043415 (2011).
  * (16) X.-J. Lu, J. G. Muga, X. Chen, U. Poschinger, F. Schmidt-Kaler, A. Ruschhaupt, Phys. Rev. A, 89, 063414 (2014).
  * (17) Q. Zhang, X. Chen, and D. Guéry-Odelin, Phys. Rev. A 92, 043410 (2015).
  * (18) Q. Zhang, J. G. Muga, D. Guéry-Odelin, X. Chen, J. Phys. B: Atomic, Molecular and Optical Physics, 49, 125503 (2016).
  * (19) X.-J. Lu, A. Ruschhaupt, and J. G. Muga, Phys. Rev. A 97, 053402 (2018).
  * (20) Q. Zhang, X. Chen, and D. Guéry-Odelin, Entropy 23, 84 (2021).
  * (21) H. Espinós, J. Echanobe, X.-J. Lu, and J. G. Muga, arXiv:2201.07555.
  * (22) S. González-Resines, D. D. Guéry-Odelin, A. Tobalina, I. Lizuain, E. Torrontegui, and J. G. Muga, Phys. Rev. Applied 8, 054008 (2017).
  * (23) E. Torrontegui, I. Lizuain, S. González-Resines, A. Tobalina, A. Ruschhaupt, R. Kosloff, and J. G. Muga, Phys. Rev. A 96, 02213 (2017).
  * (24) V. Martikyan, D. Guéry-Odelin, and D. Sugny Phys. Rev. A 101, 013423 (2020).
  * (25) V. Martikyan, A. Devra, D. Guéry-Odelin, S. J. Glaser, and D. Sugny, Phys. Rev. A 102, 053104 (2020).
  * (26) A. Le Cunuder, I. Martinez, A. Petrosyan, D. Guéry-Odelin, E. Trizac and S. Ciliberto, Appl. Phys. Lett. 109, 113502 (2016).
  * (27) G. Li, H. T. Quan, and Z. C. Tu, Phys. Rev. E 96, 012144 (2017).
  * (28) M. Chupeau, S. Ciliberto, D. Guéry-Odelin, and E. Trizac, New J. Phys. 20, 075003 (2018).
  * (29) J. A. C. Albay, S. R. Wulaningrum, C. Kwon, P.-Y. Lai, and Y. Jun, Phys. Rev. Research 1, 033122 (2019).
  * (30) The controllability Gramian matrix is defined by $\bm{W}(t_{f})=\int_{t_{0}}^{t_{f}}e^{\bm{A}(t_{f}-t^{\prime})}\bm{B}\bm{B}^{T}e^{\bm{A}^{T}(t_{f}-t^{\prime})}dt^{\prime}$.
  * (31) Y. Aharonov, J. Anandan, S. Popescu, and L. Vaidman, Phys. Rev. Lett. 64, 2965–2968 (1990).
  * (32) M. Berry, J. Phys. A: Math. Gen. 27(11), L391 (1994).
  * (33) M. Berry, Quantum Coherence and Reality (World Scientific, 1994), pp. 55–65.
  * (34) A. Kempf and A. Prain, J. Math. Phys. 58, 082101 (2017).
  * (35) D. Ma, V. Gulani, N. Seiberlich, K. Liu, J. L. Sunshine, J. L. Duerk and M. A Griswold, Nature, 495, 187 (2013).
  * (36) Q. Ansel, M. Tesch, S. J. Glaser, and D. Sugny, Phys. Rev. A 96, 053419 (2017).