#  Birth, Death, and Horizontal Flight: Malthusian flocks with an easy plane
in three dimensions

John Toner<EMAIL_ADDRESS>Department of Physics and Institute for
Fundamental Science, University of Oregon, Eugene, OR $97403$

###### Abstract

I formulate the theory of three dimensional “Malthusian flocks” – i.e.,
coherently moving collections of self-propelled entities (such as living
creatures) which are being “born” and “dying” during their motion – whose
constituents all have a preference for having their velocity vectors lie
parallel to the same two-dimensional plane. I determine the universal scaling
exponents characterizing such systems exactly, finding that the dynamical
exponent $z=3/2$, the “anisotropy” exponent $\zeta=3/4$, and the “roughness”
exponent $\chi=-1/2$. I also give the scaling laws implied by these exponents.

In addition to its obvious importance to biology, the study of “Active matter”
Active1 ; Active2 ; Active3 ; Active4 ; tissue ; PS ; Nematics  has been of
immense interest to physicists because of the novel collective phenomena that
occur in such systems. These include long-ranged orientational order in
spatial dimension $d=2$ Vicsek ; TT1 ; Chate1 ; Chate2 , and the breakdown of
linearized hydrodynamics that occurs in polar ordered active fluids ; aka
“flocks”) TT1 ; TT3 ; birdrev .

Reproduction (along with its inevitable companion, death) is another important
non-equilibrium biological phenomenon. One class of systems that exhibit both
spontaneous motion and birth and death is so-called “Malthusian flocks”Toner
(2012); Malthus_1789 ; maltd>3prl ; maltd>3pre ; that is, polar ordered active
fluids in which the number of active particles is not conserved. Such systems
are readily experimentally realizable in experiments on a, e.g., growing
bacteria colonies and cell tissuestissue , and “treadmilling” molecular motor
propelled biological macromolecules in a variety of intracellular structures,
including the cytoskletonProst_cyto , and mitotic spindlesBrugues , in which
molecules are being created and destroyed as they move.

This prior work has treated Malthusian flocks whose underlying dynamics are
completely isotropic. However, in many biological systems, anisotropy plays an
important role. For example, birds are much more inclined to fly horizontally
than vertically (for obvious evolutionary reasons!). Many oceanic creatures,
such as planktonplankton  preferentially swim horizontally as well.

In this paper, I consider the effects of such anisotropy on three-dimensional
Malthusian flocks. That is, I consider a flock that fills a three dimensional
space, whose constituents have a preference for moving parallel to an “easy
plane”, which I take to be the $x$-$y$ plane. The “hard“ direction orthogonal
to this axis I will callnot z  $r_{h}$. The two-component velocity field of
the flock $\mathbf{v}(\mathbf{r},t)$ is a function of the three-dimensional
Cartesian co-ordinate $\mathbf{r}=(x,y,r_{h})$ as well as time $t$.

I focus on what is known in the field of active matter as a “dry“ system; that
is, one without momentum conservation. In a three dimensional system as I am
considering here, this would correspond to flockers moving through a fixed,
frictional medium (e.g., worms crawling through soil).

Similarly to the isotropic theory of dry Malthusian flocks treated in Toner
(2012); Malthus_1789 ; maltd>3prl ; maltd>3pre , the hydrodynamic theory I
develop here may also be viewed as a generic non-equilibrium $3$-dimensional
$2$-component spin modelSRspin  (i.e., a non-equilibrium “XY“ model), with
space-spin couplings that are invariant under simultaneous rotations of space
and spin between the two components of a spin $\mathbf{s}(\mathbf{r},t)$ and
the $(x,y)$ coordinates of the three dimensional space. In such a spin model,
fluctuations in the system propagate spatially in a spin-direction-dependent
manner, but the spins themselves are not moving. Like the Malthusian flock,
which has no appreciable density fluctuations because the Malthusian balance
between birth and death suppresses them, such a spin system clearly has no
density fluctuations either, in this case because there is no motion.
Therefore, the only hydrodynamic variable in such a spin system is the spin
field $\mathbf{s}(\mathbf{r},t)$. Furthermore, since both the symmetries and
conservation laws of such a spin system are the same as those of the easy-
plane Malthusian flock I consider here, the hydrodynamic equation of motion
for $\mathbf{s}(\mathbf{r},t)$ is exactly the same as the one we derive here
for a Malthusian flock, with spin playing the role of the velocity field.

I focus here on the “polar ordered” state of the Malthusian flock; that is, a
state of an active polar system in which the spatially averaged velocity
vector $\langle\mathbf{v}(\mathbf{r},t)\rangle$ is non-zero in the
thermodynamic limit of $N\to\infty$, where $N$ is the number of active agents
(which I will hereafter call “flockers“). I will chose my coordinate system so
that the direction of the mean velocity is $x$; hence, $y$ is the direction
perpendicular to the mean flock motion in the easy plane.

The large-distance, long-time scaling behavior of the easy plane Malthusian
flock is characterized by three universal exponents: a “dynamical“ exponent
$z$, which relates time scales to length scales along $y$; an “anisotropy“
exponent $\zeta$, which relates length scales along both $y$ and $r_{h}$ to
length scales along $y$, and a “roughness“ exponent $\chi$, which relates the
scale of velocity fluctuations to length scales along $y$.

More precisely, the anisotropy exponent $\zeta$, relates length scales along
both $y$ and $r_{h}$ in a “boosted“ frame moving along the direction of flock
motion at a non-universal speed $\gamma$ which is not, in general, the average
speed $v_{0}$ of the flock itself, although it is likely to be of the same
order of magnitude.

That is, the exponents are defined via the “typical“ time $t(y)$, “typical“
perpendicular distance $R_{{}_{\perp}}(y)$, and “typical“ velocity fluctuation
$v_{y}(y)$

$t(y)\propto|y|^{z}\ \ \ ,\ \ \ R_{{}_{\perp}}(y)\propto|y|^{\zeta}\ \ \ ,\ \
\ v_{y}(y)\propto|y|^{\chi}\,.$ (1)

where

$R_{{}_{\perp}}=\sqrt{(x-\gamma t)^{2}+\alpha r_{h}^{2}}$ (2)

with $\alpha$ a dimensionless, non-universal (i.e., system-dependent)
parameter of $O(1)$.

These exponents can be determined from the correlations of the velocity field
$\mathbf{v}(\mathbf{r},t)$. I find that the two point velocity correlation
functions take the scaling form

$\displaystyle\langle\mathbf{v}(\mathbf{r}+\mathbf{r}^{\prime},t+T)\cdot\mathbf{v}({\mathbf{r}^{\prime},T})\rangle=|y|^{2\chi}F\left({(R_{{}_{\perp}}/a_{{}_{\perp}})\over(|y|/a_{y})^{\zeta}},{(|t|/\tau)\over(|y|/a_{y}^{z})}\right)\propto$
$\displaystyle\left\\{\begin{array}[]{ll}|y|^{2\chi}\,,&\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)\ll\left({|y|\over
a_{y}}\right)^{\zeta},{|t|\over\tau}\ll\left({|y|\over a_{y}}\right)^{z}\\\
R_{{}_{\perp}}^{2\chi\over\zeta},&\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)\gg\left({|y|\over
a_{y}}\right)^{\zeta},{|t|\over\tau}\ll\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)^{z/\zeta}\\\ \ |t|^{2\chi\over
z},&{|t|\over\tau}\gg\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)^{z/\zeta},{|t|\over\tau}\gg\left({|y|\over
a_{y}}\right)^{z}\,,\end{array}\right.$ (6)

where $\mathbf{r}=(x,y,r_{h})$,
$F\left({(R_{{}_{\perp}}/a_{{}_{\perp}})\over(|y|/a_{y})^{\zeta}},{(|t|/\tau)\over(|y|/a_{y}^{z})}\right)$
is a scaling function, which is universal up to a non-universal overall
multiplicative constant, $a_{y,\perp}$ are non-universal microscopic lengths,
and $\tau$ an equally non-universal microscopic time.

Despite the fact that this problem proves to be fundamentally non-linear, and
in contrast to both “immortal“ flocksrean , and isotropic Malthusian
flocksToner (2012); Malthus_1789 ; maltd>3prl ; maltd>3pre , I can, thanks to
the symmetries of this anisotropic problem, and the structure of its equation
of motion, determine these scaling exponents exactlyincomp . I find:

$z=3/2\ \ \ ,\ \ \ \zeta=3/4\ \ \ ,\ \ \ \chi=-1/2\,.$ (7)

Inserting these results into (6) gives

$\displaystyle\langle\mathbf{v}(\mathbf{r}+\mathbf{r}^{\prime},t+T)\cdot\mathbf{v}({\mathbf{r}^{\prime},T})\rangle=|y|^{-1}F\left({(R_{{}_{\perp}}/a_{{}_{\perp}})\over(|y|/a_{y})^{3/4}},{(|t|/\tau)\over(|y|/a_{y}^{3/2})}\right)\propto$
$\displaystyle\left\\{\begin{array}[]{ll}y^{-1}\,,&\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)\ll\left({y\over
a_{y}}\right)^{3/4},{|t|\over\tau}\ll\left({y\over a_{y}}\right)^{3/2}\\\
R_{{}_{\perp}}^{-(4/3)},&\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)^{3/4}\gg{y\over
a_{y}},{|t|\over\tau}\ll\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)^{2}\\\ \
|t|^{-(2/3)},&{|t|\over\tau}\gg\left({R_{{}_{\perp}}\over
a_{{}_{\perp}}}\right)^{2},{|t|\over\tau}\gg\left({y\over
a_{y}}\right)^{3/4}\,,\end{array}\right.$ (11)

I will now derive these results.

Equation of motion. The hydrodynamic equation of motion for a Malthusian flock
in an isotropic space was derived in Toner (2012). The hydrodynamic equation
of motion I derive here is essentially the same, except for the confinement of
the velocity field $\mathbf{v}(\mathbf{r},t)$ to lie in the $(x,y)$ plane.

Here I consider this confinement to be absolutely rigorous: in my approach,
the only non-zero components of $\mathbf{v}(\mathbf{r},t)$ are the $x$ and $y$
components $v_{x}$ and $v_{y}$.

One could imagine relaxing this constraint, and considering anisotropic flocks
in which the flockers were allowed to move in the hard direction, but with
motion along that direction suppressed by, e.g., decaying with a finite
lifetime. However, it is straightforward to show that such a suppression will
only grow upon renormalization, leading to complete suppression of the hard
component of the velocity as the system is renormalized. Physically, this
simply amounts to the statement that, if we look at the system on time scales
much longer than the relaxation time of the hard component of the velocity,
that component will be effectively zero. I therefore chose to restrict the
velocity field $\mathbf{v}$ in the hydrodynamic theory to lie in the $(x,y)$
plane.

Similar considerations apply to the magnitude $v\equiv|\mathbf{v}|$ of the
velocity $|\mathbf{v}|$. That is, I could choose to allow this magnitude to
fluctuate, as was done in, e.g., the “Toner-Tu“ equationsTT1 ; TT3 . However,
since we expect any flocking system to have a preferred speed $v_{0}$, with
fluctuations away from this speed decaying with a finite lifetime, we can
instead choose to work from the outset with a velocity which is fixed in
magnitude at $v_{0}$.polar

In addition to the above considerations, I now apply the general principles of
hydrodynamics:

I) The hydrodynamic variables are the “slow“ variables: that is, variables
that evolve slowly at long wavelengths. There are generically two reasons a
variable will be slow in this sense: first, if it is a “Goldstone mode“
associated with a spontaneously broken continuous symmetry, and second, if it
is a conserved variable.

If we define our $x$-axis to be the direction of the partially averaged
velocity $\langle\mathbf{v}\rangle$ of the flock, then the component
$v_{y}(\mathbf{r},t)$ of velocity perpendicular to that direction is a
Goldstone mode, and, hence, will be one of our hydrodynamic variables.

In flocks without birth and death, in which, therefore, the number of flockers
is conserved, the local density $\rho(\mathbf{r},t)$ of the flockers is a slow
variable for this reason. This additional variable complicates the
hydrodynamics considerably. This greater complexity makes it impossible to
determine the exact scaling exponents for such “immortal“ flocksrean ; sc ; cf
; comment ; mpp ; KPZ ; discodyn . Here, this difficulty does not arise, since
$\rho(\mathbf{r},t)$ is not conserved, and, hence, not a hydrodynamic mode.

Therefore, the only hydrodynamic variable is the velocity field
$\mathbf{v}(\mathbf{r},t)$. Its equation of motion is constrained by
rotational invariance in the $x$-$y$ plane to be:

$\displaystyle\partial_{t}v_{k}(\mathbf{r},t)=T_{ki}\bigg{(}-\lambda(\mathbf{v}\cdot{\bf\nabla})v_{i}+\mu_{{}_{BE}}\partial^{E}_{i}({\bf\nabla}\cdot\mathbf{v})+\mu_{{}_{TE}}\nabla_{{}_{E}}^{2}v_{i}+\mu_{h}\partial_{h}^{2}v_{i}+\mu_{{}_{A}}(\mathbf{v}\cdot{\bf\nabla})^{2}v_{i}+f_{i}\bigg{)}$
(12)

where the indices $i$ and $j$ run over $(x,y)$, and I have defined the
transverse projection operator

$T_{ki}\equiv\delta^{E}_{ki}-v_{k}(\mathbf{r},t)v_{i}(\mathbf{r},t)/v_{0}^{2}$
(13)

which projects any vector orthogonal to $\mathbf{v}(\mathbf{r},t)$. Its
presence in (12) insures that the fixed length condition
$|\mathbf{v}(\mathbf{r},t)|=v_{0}$ on $\mathbf{v}(\mathbf{r},t)$ is preserved.

The $\mathbf{f}$ term in (12) is a random Gaussian white noise, reflecting
errors made by the flockers, with correlations:

$\displaystyle\langle
f_{i}(\mathbf{r},t)f_{j}(\mathbf{r}^{\prime},t^{\prime})\rangle=2D\delta^{E}_{ij}\delta^{d}(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$
(14)

where the noise strength $D$ is a constant hydrodynamic parameter and $i,j$
label vector components.

In the ordered state (defined as the state in which
$\left<\mathbf{v}(\mathbf{r},t)\right>=v_{0}{\bf\hat{x}}\neq{\bf 0}$, where
I‘ve chosen the spontaneously picked direction of mean flock motion as our
$x$-axis), we can write the velocity as:

$\displaystyle\mathbf{v}(\mathbf{r},t)=\sqrt{v_{0}^{2}-v_{y}^{2}(\mathbf{r},t)}\,{\bf\hat{x}}+v_{y}(\mathbf{r},t){\bf\hat{y}}\,.$
(15)

Taking the $y$-component of the equation of motion (12), changing co-ordinates
to a new Galilean frame $\mathbf{r}^{\prime}$ moving with respect to our
original frame in the direction ${\bf\hat{x}}$ of mean flock motion at speed
$\gamma=\lambda v_{0}$ – i.e., $\mathbf{r}^{\prime}\equiv\mathbf{r}-\gamma
t{\bf\hat{x}}$, and, finally, expanding to cubic order in $v_{y}$ and dropping
terms of order $v_{y}^{2}$ with more than one spatial derivative, gives

$\displaystyle\partial_{t}v_{y}+{\lambda\over
2}\partial_{y}(v_{y}^{2})+\lambda_{3}\partial_{x}(v_{y}^{3})$ $\displaystyle=$
$\displaystyle\bigg{(}\mu_{x}\partial_{x}^{2}+\mu_{y}\partial_{y}^{2}+\mu_{h}\partial_{h}^{2}\bigg{)}v_{y}+f_{y}\
,$ (16)

where I have defined $\lambda_{3}\equiv-{\lambda\over 6v_{0}^{2}}$, and the
diffusion constants $\mu_{x,y}$ are given by

$\mu_{x}=\mu_{{}_{TE}}+\mu_{{}_{A}}v_{0}^{2}\ \ \ ,\ \ \
\mu_{y}=\mu_{{}_{BE}}+\mu_{{}_{TE}}\,.$ (17)

Stability of the homogeneous ordered state requires that $\mu_{x,y,r_{h}}$ are
all positive.

I will begin by studying the equation obtained from (16) by dropping the cubic
$\lambda_{3}$ term, which reads:

$\displaystyle\partial_{t}v_{y}+{\lambda\over 2}\partial_{y}(v_{y}^{2})$
$\displaystyle=$
$\displaystyle\bigg{(}\mu_{x}\partial_{x}^{2}+\mu_{y}\partial_{y}^{2}+\mu_{h}^{2}\partial_{h}^{2}\bigg{)}v_{y}+f_{y}\,.$

Dynamic renormalization group (DRG) analysis.

I now perform a DRG analysis FNS  on the equation of motion (LABEL:nocub).
Specifically, I first average over short wavelength degrees of freedom, and
then perform the following rescalinghard rescale :

$(x,r_{h})\mapsto{\rm e}^{\zeta\ell}(x,r_{h})\ ,\ y\mapsto{\rm e}^{\ell}y\ ,\
t\mapsto{\rm e}^{z\ell}t\ ,\ v_{y}\mapsto{\rm e}^{\chi\ell}v_{y}\ .$ (19)

There are a number of features of the structure of the reduced equation of
motion (LABEL:nocub) that prevent certain parameters appearing in
(LABEL:nocub) from being renormalized by the first step of the DRG process
(that is, the averaging over the short-wavelength degrees of freedom). In the
jargon of the field, these parameters get no “graphical” renormalization.

The first such parameter is the coefficient $\lambda$ of the only nonlinear
term in the equation of motion, $\lambda\partial_{y}(v_{y}^{2})$. This gets no
graphical renormalization because of the inherent pseudo-Galilean symmetry of
(LABEL:nocub), i.e., the invariance of the equation of motion under the
simultaneous replacements: $y\mapsto y+t\lambda w$ and $v_{y}\mapsto v_{y}+w$
for any arbitrary constant $w$.

The parameters $\mu_{x,h}$ also get no renormalization, because that sole non-
linearity $\lambda\partial_{y}(v_{y}^{2})$ is a total $y$-derivative, and,
hence, can only renormalize terms in the equation of motion that are also
total $y$-derivatives, which the $\mu_{x,h}$ are clearly not.

Finally, the noise strength $D$ gets no graphical renormalization for
essentially the same reason: it also involves no $y$-derivatives.

In light of these considerations, the recursion relations for the
aforementioned parameters can be obtained purely from the rescaling step of
the DRG, and are easily shown to read:

$\displaystyle\frac{1}{D}\frac{{\rm d}D}{{\rm d}\ell}$ $\displaystyle=$
$\displaystyle z-2\chi-d+1-\zeta\,,$ (20)
$\displaystyle\frac{1}{\lambda}\frac{{\rm d}\lambda}{{\rm d}\ell}$
$\displaystyle=$ $\displaystyle z+\chi-1\,,$ (21)
$\displaystyle\frac{1}{\mu_{x,z}}\frac{{\rm d}\mu_{x,z}}{{\rm d}\ell}$
$\displaystyle=$ $\displaystyle z-2\zeta\,,$ (22)

exactly.

To obtain a DRG fixed point, all three of these derivatives must vanish. This
gives us three independent linear equations in the three unknown exponents
$z$, $\zeta$, and $\chi$, whose solutions are equation (7).

With these exponents in hand, a straightforward and entirely standard DRG
analysis implies that the two-point velocity-velocity correlation function has
the form (11).

Proving the irrelvance of $\lambda_{3}$. Since both the $\lambda$ and the
$\lambda_{3}$ vertex in equation blah actually originate from the same
rotation-invariant $(\mathbf{v}\cdot\nabla)\mathbf{v}$ term in our original
equation of motion blah, rotational invariance implies that the ratio of these
two term remains fixed at $-6v_{0}^{2}$. Since our system has long-ranged
orientational order, which is equivalent to a non-zero average velocity,
$v_{0}$ can get no significant renormalization at large distances. Thus, the
ratio of $\lambda$ to $\lambda_{3}$ must remain unchanged by the graphical
corrections to them. Therefore, the DRG recursion relations for these two
parameters must read:

$\displaystyle\frac{d\ln\lambda}{d\ell}$
$\displaystyle=z+\chi-1+f(\lambda_{3},...)\,,$
$\displaystyle\frac{d\ln\lambda_{3}}{d\ell}$
$\displaystyle=z+2\chi-\zeta+f(\lambda_{3},...)\,,$ (23)

where $f(\lambda_{3},...)$ denotes the graphical corrections.AS

My earlier argument that $\lambda$ gets no graphical correction when
$\lambda_{3}=0$ implies that $f(\lambda_{3},...)$ vanishes as $\lambda_{3}\to
0$. Therefore, it is at most $O(\lambda_{3})$ for small $\lambda_{3}$. Hence
the recursion relation (23) for $\lambda_{3}$ implies

$\frac{d\lambda_{3}}{d\ell}=(z+2\chi-\zeta)\lambda_{3}+O(\lambda_{3}^{2})=-{1\over
4}\lambda_{3}+O(\lambda_{3}^{2})\,.$ (24)

Thus I see that a sufficiently small $\lambda_{3}$ is always irrelevant, with
RG eigenvalue $-1/4$. This justifies dropping the $\lambda_{3}$ term, as I did
earlier.

Summary. I have formulated the hydrodynamic theory of three dimensional
“Malthusian flocks” moving with an easy plane. This problem proves to be one
of the rare problemsincomp  for which, despite strongly relevant non-
linearities in the equation of motion, it is possible to determine the
universal scaling exponents characterizing it exactly, finding that the
dynamical exponent $z=3/2$, the “anisotropy” exponent $\zeta=3/4$, and the
“roughness” exponent $\chi=-1/2$. I‘ve also determined the scaling form of the
velocity-velocity two point correlation function implied by these exponents.

###### Acknowledgements.

I am grateful to Sriram Ramaswamy and Ananyo Maitra for extremely valuable
discussions, without which I would not have been able to demonstrate the
irrelevance of the cubic vertex $\lambda_{3}$, on which my entire argument for
exact exponents depends. I also thank The Higgs Centre for Theoretical Physics
at the University of Edinburgh, where this work was done, for their
hospitality and support.

## References

  * (1) J. Toner, The Physics of Flocking: Birth, Death, and Flight in Active Matter, Cambridge University Press, (2024).
  * (2) S. Ramaswamy, The mechanics and statics of active matter. Ann. Rev. Condens. Matt. Phys. 1, 323-345 (2010).
  * (3) M.C. Marchetti, J.F. Joanny, S. Ramaswamy, T.B. Liverpool, J. Prost, M. Rao, and R.A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143-1188 (2013).
  * (4) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Active particles in complex and crowded environments, Rev. Mod. Phys. 88, 045006 (2016).
  * (5) F. Schweitzer, Brownian Agents and Active Particles: Collective Dynamics in the Natural and Social Sciences. Springer Series in Synergetics (Springer, New York, 2003).
  * (6) J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. Jülicher, Fluidization of Tissues by Cell Division and Apoptosis, Proc. Natl. Acad. Sci. USA 107, 20863 (2010); J. Ranft, J. Prost, F. Jülicher, and J.-F. Joanny, Tissue Dynamics with Permeation, Eur. Phys. J. E 35, 46 (2012).
  * (7) M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur, Arrested phase separation in reproducing bacteria creates a generic route to pattern formation, Proc. Natl. Acad. Sci. USA 107(26), 11715 (2010).
  * (8) D. Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. C. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo and W. C. K. Poon, A growing bacterial colony in two dimensions as an active nematic, Nature Communications 9, 4190 (2018); Z. You, D. J. G. Pearce, A. Sengupta, and L. Giomi, Phys. Rev. X 8, 031065 (2018); A. Doostmohammadi, S. P. Thampi, and J. M. Yeomans, Defect-Mediated Morphologies in Growing Cell Colonies, Phys. Rev. Lett. 117, 048102 (2016).
  * (9) T. Vicsek, A. Czirók, E. Ben-jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-Driven particles. Phys. Rev. Lett. 75, 1226 (1995).
  * (10) J. Toner, and Y. Tu, Long-range order in a two-dimensional dynamical XY model: how birds fly together. Phys. Rev. Lett. 75, 4326 (1995).
  * (11) G. Grégoire and H. Chaté, Onset of Collective and Cohesive Motion. Phys. Rev. Lett. 92, 025702 (2004).
  * (12) H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Collective motion of self-propelled particles interacting without cohesion. Phys. Rev. E 77, 046113 (2008).
  * (13) J. Toner, and Y. Tu, Flocks, herds, and schools: a quantitative theory of flocking. Phys. Rev. E 58, 4828(1998).
  * (14) J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks. Ann. Phys. 318, 170 (2005).
  * Toner (2012) J. Toner, Birth, Death, and Flight: A Theory of Malthusian Flocks. Physical Review Letters 108, 088102 (2012), URl http://dx.doi.org/10.1103/PhysRevLett.108.088102. Sadly, the argument in this paper‘s claim for exact exponents proves to be incorrect. I thank A. Maitra and S. Ramaswamy for pointing the flaw in this argument out to me.
  * (16) T. R. Malthus, An Essay on the Principle of Population, edited by J. Johnson (St. Paul’s Churchyard, London, 1798).
  * (17) L. Chen, C. F. Lee, and J. Toner, Moving, reproducing, and dying beyond Flatland: Malthusian flocks in dimensions $d>2$, Phys. Rev. Lett. 125, 098003 (2020).
  * (18) L. Chen, C. F. Lee, and J. Toner, A novel nonequilibrium state of matter: a $d=4-\epsilon$ expansion study of Malthusian flocks, Phys. Rev. E 102, 022610 (2020).
  * (19) K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Generic Theory of Active Polar Gels: A Paradigm for Cytoskeletal Dynamics, Eur. Phys. J. E 16, 5 (2005).
  * (20) J. Brugués and D. Needleman, Physical basis of spindle self-organization, PNAS 111 (52) 18496 (2014).
  * (21) R. Benzi, D. R. Nelson, S, Shankar, F. Toschi, and X. Zhu, (2022). Spatial population genetics with fluid flow. Reports on Progress in Physics, 85(9), 096601 (2022).
  * (22) I refrain from the the obvious choice of calling this coordinate “$z$“ because I reserve the symbol $z$ for the dynamical exponent.
  * (23) This conclusion was independently reached for isotropic systems by L. P. Dadhichi, J. Kethapelli, R. Chajwa, S. Ramaswamy, and A. Maitra, Nonmutual torques and the unimportance of motility for long-range order in two-dimensional flocks, Phys. Rev. E 101, 052601.
  * (24) J. Toner, Reanalysis of the hydrodynamic theory of fluid, polar-ordered flocks. Phys. Rev. E 86, 031918-1-031918-9 (2012).
  * (25) It has proven possible to determine exact exponents for a number of incompressible systems. See, e.g., L. Chen, C. F. Lee, and J. Toner, Squeezed in three dimensions, moving in two: Hydrodynamic theory of 3D incompressible easy-plane polar active fluids, Phys. Rev. E 98, 040602(R) (2018); L. Chen, C. F. Lee, and J. Toner, L. Chen, Incompressible polar active fluids in the moving phase in dimensions $d>2$, arXiv:1806.02548, New J. Phys. 20, 113035 (2018); L. Chen, C.F. Lee, A. Maitra, and J. Toner, Packed swarms on dirt: two dimensional incompressible flocks with quenched and annealed disorder, arXiv:2202.02865, and Phys. Rev. Lett. 129, 188004 (2022); L. Chen, C.F. Lee, A. Maitra, and J. Toner, Incompressible polar active fluids with quenched disorder in two dimensions, arXiv:2205.14384, and Phys. Rev. E 106, 044608 (2022); L. Chen, C.F. Lee, A. Maitra, and J. Toner, Incompressible polar active fluids with quenched disorder in dimensions $d>2$, arXiv:2203.01892, and Phys. Rev. Lett. 129, 198001 (2022). However, since the systems considered in all of the above work are incompressible, they correspond to an extreme - indeed, a singular, and unstable under renormalization, limit of one of the parameters: namely, the limit of inverse compressibility going to infinity. In contrast, all of the parameters in the theory Malthusian flocks are finite, except for those that are driven to infinity by renormalization.
  * (26) Alternatively, one could instead choose to take as the hydrodynamic variable in the problem the unit vector $\mathbf{\hat{p}}\equiv\mathbf{v}/|\mathbf{v}|$ along the velocity. This variable, often called “the polarization“, is, indeed, often used in the field of Active Matter. I choose not to do so here. The difference between these choices is irrelevant in both the renormalization group and the conventional english sense.
  * (27) H. Chate and A. Solon, Dynamic Scaling of Two-Dimensional Polar Flocks, Dynamic Scaling of Two-Dimensional Polar Flocks, Phys. Rev. Lett. 132, 268302 (2024).
  * (28) C. F. Lee and P. Jentsch, A new universality class describes Vicsek‘s flocking phase in physical dimensions, arXiv:2402.01316.
  * (29) Both references sc  and cf  claim to obtain exact exponents for the immortal flocks problem; the former only in $d=2$; the latter in all $d$. The answers claimed in these two references are slightly different. I believe that sc  is incorrect, since it assumes (incorrectly, in my view) that the time derivative of the Goldstone mode must be a total divergence, which is not the case. Numerous counter-examples to this claim exist: the KPZ equationKPZ  for one, and the hydrodynamics of crystalline solidsmpp , smecticsmpp , and discotic phasesdiscodyn  for three others. cf  is more promising, but still depends on a simplification of the full hydrodynamic equations that has not yet been completely justified. It also has difficulties in spatial dimensions $d=2$.
  * (30) M. Kardar, G. Parisi and Y.-C. Zhang, “Dynamic scaling of growing interfaces,” Physical Review Letters 56, 889-892 (1986).
  * (31) P.C. Martin, O. Parodi, and P.S. Pershan, Unified hydrodynamic theory for crystals, liquid crystals, and normal fluids, Phys. Rev. A 6, 2401 (1972).
  * (32) J. Prost and N. A. Clark, in Proceedings of the International Conference on Liquid Crystals, Bangalore (1979), edited by S. Chandrasekhar, Heyden, Philadelphia (1980); in Liquid Crystals of One- and Two-Dimensional Order, Garmisch-Partenkirchen, Fanvar (1980).
  * (33) D. Forster, D. R. Nelson and M. J. Stephen, Large-distance and long-time properties of a randomly stirred fluid. Phys. Rev. A 16 732 (1977).
  * (34) More generally, I could have chosen different anisotropy exponents $\zeta_{x,h}$ for the $x$ and hard directions respectively. However, doing so, one quickly finds that a fixed point can only be obtained by choosing $\zeta_{x}=\zeta_{h}$.
  * (35) This aproach to determining the relevance of $\lambda_{3}$, was first obtained, for the somewhat different problem of a two-dimensional Malthusian flock, by Sriram Ramaswamy and Ananyo Maitra. I thank them for patiently explaining this argument to me.


# Quantum control with a multi-dimensional Gaussian quantum invariant

Selwyn Simsek  Physics Department, Blackett Laboratory, Imperial College
London, Prince Consort Road, SW7 2BW, United Kingdom  Florian Mintert  Physics
Department, Blackett Laboratory, Imperial College London, Prince Consort Road,
SW7 2BW, United Kingdom

(March 8, 2021)

###### Abstract

The framework of quantum invariants is an elegant generalization of adiabatic
quantum control to control fields that do not need to change slowly. Due to
the unavailability of invariants for systems with more than one spatial
dimension, the benefits of this framework have not yet been exploited in
multi-dimensional systems. We construct a multi-dimensional Gaussian quantum
invariant that permits the design of time-dependent potentials that let the
ground state of an initial potential evolve towards the ground state of a
final potential. The scope of this framework is demonstrated with the task of
shuttling an ion around a corner which is a paradigmatic control problem in
achieving scalability of trapped ion quantum information technology.

##  1 Introduction

The development of hardware for quantum information processing has reached a
stage in which tens of qubits can be accurately controlled [Friis  _et
al._(2018)Friis, Marty, Maier, Hempel, Holzäpfel, Jurcevic, Plenio, Huber,
Roos, Blatt, and Lanyon, Wei  _et al._(2020)Wei, Lauer, Srinivasan,
Sundaresan, McClure, Toyli, McKay, Gambetta, and Sheldon]. In order to control
such a large register of qubits, it is advisable to arrange them into several
smaller units that are coherently interconnected. Such interconnections can be
realised in a multitude of different ways [Monroe  _et al._(2014)Monroe,
Raussendorf, Ruthven, Brown, Maunz, Duan, and Kim, Lekitsch  _et
al._(2017)Lekitsch, Weidt, Fowler, Mølmer, Devitt, Wunderlich, and Hensinger].
In the case of trapped ions, one of the most prominent envisioned
architectures is based on segmented traps, in which ions can be moved, i.e.
shuttled, between different segments [Kielpinski  _et al._(2002)Kielpinski,
Monroe, and Wineland, Lekitsch  _et al._(2017)Lekitsch, Weidt, Fowler, Mølmer,
Devitt, Wunderlich, and Hensinger] of a segmented trap.

Since trapped ion quantum logic usually requires the motional states of ions
to be close to their quantum mechanical ground state [Mølmer and
Sørensen(1999)], it is desirable that any such shuttling process should
transfer ions to their motional ground state with high fidelity. Performing
all operations of a quantum algorithm within the time-scale imposed by the
system’s coherence times, requires fast shuttling [Kielpinski  _et
al._(2002)Kielpinski, Monroe, and Wineland] implying that the transient
motional states during shuttling are far away from the motional ground state.
Any shuttling process thus needs to end with a stage of deceleration in which
an originally rapidly moving ion is being transferred to its quantum
mechanical motional ground state.

In practice, shuttling is realised in terms of time-dependent trapping
potentials resulting from voltages applied to trap electrodes [Rowe  _et
al._(2002)Rowe, Ben-Kish, Demarco, Leibfried, Meyer, Beall, Britton, Hughes,
Itano, Jelenković, Langer, Rosenband, and Wineland, Hensinger  _et
al._(2006)Hensinger, Olmschenk, Stick, Hucul, Yeo, Acton, Deslauriers, Monroe,
and Rabchuk]. Describing the quantum dynamics of a particle through such a
time-dependent potential landscape is generally a formidable numerical
challenge, but in the present context the problem can be substantially
simplified. The initial ground state is a Gaussian wavepacket characterised by
the expectation values of position and momentum and the corresponding
covariances of those operators. For realistic trapping potentials, the
potential can be approximated to be harmonic in the spatial domain occupied by
such a wavepacket [Leibfried  _et al._(2003)Leibfried, Blatt, Monroe, and
Wineland]. Within this approximation, an initial Gaussian wavepacket remains
Gaussian under time evolution, and thus the state is characterised in terms of
expected position in phase space and covariances, i.e. a small number of
parameters [Heller(1975a), Heller(1975b), Heller(1976), Heller(1977),
Heller(1981)].

Identifying time-dependent gate voltages that realise a shuttling protocol
that transfers one or several ions from their ground state in the initial
potential to the ground state of their final potential, typically requires
testing a large number of different potential solutions. As such, the
reduction of numerical effort resulting from the Gaussian approximation is
invaluable.

An equally important reduction in effort can be achieved with the framework of
quantum invariants [Lewis(1968)], since those allow one to directly construct
shuttling protocols that end with the ions in their motional ground state
[Chen  _et al._(2011)Chen, Torrontegui, and Muga]. The framework of quantum
invariants and its use for control of wave-packets via invariant-based inverse
engineering is well-established for one-dimensional problems [Chen  _et
al._(2011)Chen, Torrontegui, and Muga, Torrontegui  _et al._(2011)Torrontegui,
Ibáñez, Chen, Ruschhaupt, Guéry-Odelin, and Muga, Levy  _et al._(2018)Levy,
Kiely, Muga, Kosloff, and Torrontegui, Torrontegui  _et al._(2013)Torrontegui,
Ibáñez, Martínez-Garaot, Modugno, del Campo, Guéry-Odelin, Ruschhaupt, Chen,
and Muga, Tobalina  _et al._(2017)Tobalina, Palmero, Martínez-Garaot, and
Muga].

Invariant-based inverse engineering has been also suggested to realise the
transport of cold atoms in optical lattices [Lu  _et al._(2020)Lu, Ruschhaupt,
Martínez-Garaot, and Muga] or optical tweezers [Zhang  _et al._(2016)Zhang,
Muga, Guéry-Odelin, and Chen] as well as atom clouds with a corresponding
experimental implementation [Ness  _et al._(2018)Ness, Shkedrov, Florshaim,
and Sagi]. Invariant-based inverse engineering has also motivated the study of
fast transport of spin-orbit-coupled Bose-Einstein condensates [Chen  _et
al._(2018)Chen, Jiang, Li, Ban, and Sherman], and quantum invariants have been
used as a theoretical tool to compute topological phases in planar waveguides
[Rosu and Romero(1999)] and minispace quantum cosmologies [Espinoza(2000)].
They may also be constructed for light beam propagation in nonlinear
inhomogeneous media [Goncharenko  _et al._(1991)Goncharenko, Logvin, Samson,
Shapovalov, and Turovets].

Processes in which an ion is being shuttled across a junction in a more
complicated trap geometry, however, cannot yet be approached with invariant-
based inverse engineering, because the development of invariants for systems
with more than one translational degree of freedom [Tobalina  _et
al._(2020)Tobalina, Torrontegui, Lizuain, Palmero, and Muga]is still in its
infancy. Here, we develop a framework of quantum invariants that is suitable
for optimal control of Gaussian wave-packets in any number of spatial
dimensions via invariant-based inverse engineering.

In Sec. 2, invariant-based inverse engineering and Gaussian states are
summarised. Sec. 3 introduces the invariant that is at the core of this paper
and discusses how this invariant can be used to design time-dependent
potentials that result in ground-state-to-ground-state shuttling. This section
is followed by explicit examples of such protocols discussed in Sec. 4, while
the actual derivation of the invariant is in Sec. 5, and Sec. 6 contains a
discussion of how well-established one-dimensional invariants are contained as
special case in the present framework. A reader who is more interested in the
underlying framework than in its application to a specific control problem,
can skip Sec. 4 or read Secs. 5 and 6 in any order.

##  2 Invariant-based control

In this section several elements of Gaussian quantum states and quantum
invariants upon which the present work builds will be reviewed.

In order to ease notation, from hereon in the momentum variables $p_{i}$ will
include a scaling factor depending on the mass of each particle, such that the
kinetic energy reads $\sum_{i}\hat{p}_{i}^{2}/(2m)$ instead of the usual
expression $\sum_{i}\hat{p}_{i}^{2}/(2m_{i})$ that explicitly reflects the
possibility of having different ions with different masses. A corresponding
scaling of the position variables $x_{i}$ ensures that $\hat{x}_{i}$ and
$\hat{p}_{i}$, as well as their classical counterparts, are conjugate to each
other.

###  2.1 Gaussian dynamics

The framework uses the Gaussian approximation of wavepackets that simplifies
the description in terms of general quantum states to a description in terms
of classical trajectories and covariances. If $\hat{X}$ denotes the
$2d$-dimensional vector comprised of the phase space operators of a
$d$-dimensional system, i.e. $d$ position operators and $d$ momentum operators

$\hat{X}=\begin{pmatrix}\hat{r}_{1},\hat{r}_{2},...,\hat{r}_{n},\hat{p}_{1},\hat{p}_{2},...,\hat{p}_{n}\end{pmatrix},$
(1)

then any Gaussian wavepacket is completely characterised in terms of the
expectation value of $\hat{X}$, denoted in the following by the vector $X$,
and by the $2d\times 2d\ -$ dimensional covariance matrix $\Sigma$ with the
elements

$\Sigma_{ij}=\frac{1}{2}\langle\hat{X}_{i}\hat{X}_{j}+\hat{X}_{j}\hat{X}_{i}\rangle-\langle\hat{X}_{i}\rangle\langle\hat{X}_{j}\rangle\
,$ (2)

where the symbol $\langle\ \rangle$ denotes the expectation value with respect
to the Gaussian quantum state.

In a potential that is not always harmonic, an initially Gaussian quantum
state will lose its Gaussian character, but if the wavepacket is sufficiently
well localised in real-space, then the potential can be approximated by its
second-order Taylor expansion around the center of the wavepacket, and the
quantum state remains Gaussian [Heller(1975b)].

The center of the wavepacket, $X=\langle\hat{X}\rangle$, follows the classical
trajectory, i.e. it satisfies the equation of motion

$\dot{X}={\cal S}\vec{V}\ ,$ (3)

where

$\vec{V}=\nabla H_{c}$ (4)

is the vector containing the first derivatives of the classical Hamiltonian
$H_{c}$ with respect to the phase space variables, evaluated at the centre of
the wavepacket $X$, and the symplectic matrix

${\cal S}=\left[\begin{array}[]{rr}\mathbbm{O}&\mathbbm{1}\\\
-\mathbbm{1}&\mathbbm{O}\\\ \end{array}\right]\ ,$ (5)

is defined in terms of the $d\times d-$dimensional identity matrix
$\mathbbm{1}$ and zero-matrix $\mathbbm{O}$.

The covariance matrix satisfies the equation of motion [Mintert and
Heller(2009)]

$\dot{\Sigma}=\Sigma\Omega{\cal S}-{\cal S}\Omega\Sigma\\\ $ (6)

in terms of the matrix $\Omega$ containing the second derivatives

$\Omega_{ij}=\frac{\partial^{2}H_{c}}{\partial x_{i}\partial x_{j}}$ (7)

of the classical Hamiltonian $H_{c}$ with respect to the phase space variables
evaluated at $X$.

From hereon in it is assumed that the anharmonic component of $H$, which is to
say the Taylor expansion of $H$ beyond the second order, is negligible, in
which case the approximation holds and the system Hamiltonian can be expressed
as

$H=\frac{1}{2}\hat{X}^{T}\Omega\hat{X}+\vec{V}\hat{X}\ .$ (8)

###  2.2 Quantum invariants

The problem at hand lies in finding $\Omega$ and $\vec{V}$ such that $X(T)$
matches the desired phase space position at the final instance in time $T$ of
the shuttling protocol, and such that $\Sigma(T)$ corresponds to the
covariances of the quantum mechanical ground state of $H(T)$. Quantum
invariants allow one to derive straightforwardly expressions for $\Omega$ and
$\vec{V}$ that ensure that the boundary conditions at $t=0$ and $t=T$ are
satisfied [Tobalina  _et al._(2017)Tobalina, Palmero, Martínez-Garaot, and
Muga].

A quantum invariant is any operator ${\cal I}(t)$ satisfying the equation of
motion [Lewis(1967)]

$\frac{\partial{\cal I}(t)}{\partial t}=i[{\cal I}(t),H(t)]\ .$ (9)

Crucially, the instantaneous eigenstates of an invariant, are solutions of the
Schrödinger equation with the time-dependent Hamiltonian $H(t)$ [Lewis(1968)].
A quantum invariant with non-degenerate spectrum that commutes with the
Hamiltonian at the beginning and the end of a time interval of interest, i.e.
$[{\cal I}(0),H(0)]=0$ and $[{\cal I}(T),H(T)]=0$, can thus be used to
identify shuttling protocols that ensure that the system ends up in an
eigenstate of the final Hamiltonian $H(T)$, if it is initialised in an
eigenstate of $H(0)$ [Chen  _et al._(2011)Chen, Torrontegui, and Muga]. If the
system is initialised in the ground state of $H(0)$, then it will be
transferred to the ground state of $H(T)$ since the only eigenstate of the
quadratic operator $H(T)$ that is Gaussian is the ground state. This is
similar to the idea of adiabatic transitions, but in contrast, the invariant
based method is exact and holds also under fast changes in the Hamiltonian,
and the system state remains an eigenstate of the invariant ${\cal I}(t)$
rather than an eigenstate of the system Hamiltonian $H(t)$.

###  2.3 Control with quadratic invariants

In the case of quadratic, time-dependent, Hamiltonians, quantum invariants can
be constructed from which one can deduce quadratic Hamiltonians $H(t)$ that
correspond to those invariants. This, however, is not very useful in practice
since, the resulting Hamiltonians contain in general bilinear terms of phase
space operators; in particular, the part corresponding to the kinetic energy
will be of the form $\sum_{ij}h_{ij}(t)\hat{p}_{i}\hat{p}_{j}$, which is not
necessarily equal to the desired term $\sum_{i}\hat{p}_{i}^{2}/2m$.

Since, in practice, shuttling must be realised in terms of modulation of only
the trapping potential, it is crucial to ensure that Hamiltonians resulting
from the invariant-based framework are of the form

$H(t)=\sum_{i}\frac{\hat{p}_{i}^{2}}{2m}+\frac{1}{2}m\sum_{ij}M_{ij}\hat{x}_{i}\hat{x}_{j}-\sum_{i}\vec{F}_{i}\hat{x}_{i}\
,$ (10)

with tuneable parameters $\vec{F}$ and $M$, and no Lorentz-type terms of the
form $\hat{x}_{i}\hat{p}_{j}+\hat{x}_{j}\hat{p}_{i}$, or terms linear in the
momenta $\hat{p}_{i}$.

In one-dimensional systems, there does exist a class of quantum invariants
[Chen  _et al._(2011)Chen, Torrontegui, and Muga] that ensures the
Hamiltonians are of this form and gives a way to deduce $M$ and $\vec{F}$ from
the invariant. Despite the utility and applicability of these one-dimensional
invariants, generalisations to higher dimensional systems have not been found
to date.

##  3 Optimal control with high-dimensional quadratic invariants

As the present work deals with quadratic Hamiltonians, it is most useful to
work with an invariant that contains the phase space operators $\hat{X}$ only
up to second order. Any such invariant can be written as

${\cal I}(t)=\frac{1}{2}\hat{X}^{T}\Gamma\hat{X}+\vec{W}\cdot\hat{X}+\theta\
,$ (11)

parametrised in terms of a time-dependent, Hermitian matrix $\Gamma(t)$, a
real vector $\vec{W}(t)$, and a real scalar $\theta(t)$.

${\cal I}(t)$ is an invariant if and only if Eq. (9) holds, which is the case
if the equations

$\displaystyle\frac{d\Gamma}{dt}$ $\displaystyle=\Omega{\cal
S}\Gamma-\Gamma{\cal S}\Omega,$ (12) $\displaystyle\frac{d\vec{W}}{dt}$
$\displaystyle=\Omega{\cal S}\vec{W}-\Gamma{\cal S}\vec{V},$ (13)
$\displaystyle\frac{d\theta}{dt}$ $\displaystyle=\vec{V}^{T}{\cal S}\vec{W},$
(14)

with $\Omega$ and $\vec{V}$ as defined in Eqs.(3) and (7) are satisfied. For
any choice of $\Gamma$ and $\vec{W}$, one may deduce choices of $\Omega$,
$\vec{V}$ and $\theta$ which satisfy Eqs. (12), (13) and (14).

This approach would permit us to find a quadratic Hamiltonian for an invariant
specified by functions $\Gamma(t)$ and $\vec{W}(t)$ that one would be free to
choose. For most choices of such invariants, however, any compatible
Hamiltonian will not necessarily be of the form given in Eq. (10), i.e. it
would not be realizable in practice. It is thus essential to find instances of
$\Gamma$ and $\vec{W}$, such that $\Omega$ is of the form

$\Omega=\left[\begin{array}[]{cc}mM&\mathbbm{O}\\\
\mathbbm{O}&\frac{1}{m}\mathbbm{1}\end{array}\right]\ ,$ (15)

and such that $\vec{V}$ is of the form

$\vec{V}=\left[\begin{array}[]{c}-\vec{F}\\\ 0\end{array}\right]\ .$ (16)

###  3.1 A suitable invariant

While a general invariant of the form specified in Eq. (11) does not
correspond to a Hamiltonian of desired form, the following explicit
parametrization forces the Hamiltonian to be of the form of Eq. (10), as shown
in Sec.5.

The quadratic part $\Gamma$ of the invariant is given by

$\begin{array}[]{l}\Gamma=\\\ \left[\begin{array}[]{cc}\displaystyle
m\left(\dot{R}^{2}+\Re\left([\dot{R},R]_{A}-RA^{2}R\right)\right)&\displaystyle\frac{{\cal
J}-\\{R,\dot{R}\\}}{2}\\\\[6.0pt] \displaystyle\frac{-{\cal
J}-\\{R,\dot{R}\\}}{2}&\displaystyle\frac{R^{2}}{m}\end{array}\right]\end{array}$

where the symbol $[x,y]_{z}=xzy-yzx$ denotes a generalized commutator and
$\\{x,y\\}=xy+yx$ is the anti-commutator.

$\Gamma$ is parametrised in terms of a $d\times d$ real, positive and semi-
definite square matrix $R$ satisfying $\dot{R}(0)=0$, and the matrices $A$ and
${\cal J}$ are determined uniquely from the relations

$\displaystyle A$ $\displaystyle=$ $\displaystyle
iR^{-2}+\frac{1}{2}[R^{-1},\dot{R}]+\frac{1}{2}R^{-1}{{\cal J}}R^{-1}\ ,$
$\displaystyle\\{{{\cal J}},R^{-2}\\}$ $\displaystyle=$
$\displaystyle[\dot{R},R^{-1}]+[R,R^{-2}]_{\dot{R}}\ .$ (17)

The linear part $\vec{W}$ of the invariant ${\cal I}$ is given by

$\vec{W}=-\Gamma\left[\begin{array}[]{c}\vec{L}\\\
m\dot{\vec{L}}\end{array}\right]\ ,$ (18)

parametrized by a d-dimensional real vector $\vec{L}$,

and the scalar part is given by

$\theta(t)=\frac{1}{2}\vec{W}^{T}\Gamma^{-1}\vec{W}\ .$ (19)

The Hamiltonian corresponding to an invariant is determined by the equations
of motion Eqs. (12) and (13). As shown in the following, the quantities $M$
and $\vec{F}$ (defined in Eq. (15) and (16)) used to parametrize and
accessible Hamiltonians are determined by Eqs. (20) and (21).

The quadratic component $M$ is determined by

$\\{R^{2},M\\}=2[\dot{R},R]_{A}-\\{\ddot{R},R\\}-2RA^{2}R$ (20)

in terms of $R$ following Eq. (50). This is still not an explicit solution for
$M$, but it can be obtained in terms of an expansion into the eigenstates of
$R$ as stated in Eq. (53).

Once $M$ is obtained, the linear component $\vec{F}$ is obtained through the
relation

$\vec{F}=m\left(\ddot{\vec{L}}+M\vec{L}\right)\,$ (21)

following Eq. (61).

$R$ and $\vec{L}$ may be chosen freely subject to a set of boundary conditions
discussed in the next section.

###  3.2 Boundary conditions

As stated in Sec. 2.2, a prerequisite for invariant-based control is that the
invariant and Hamiltonian commute at initial and final times, i.e.

$[{\cal I}(0),H(0)]=[{\cal I}(T),H(T)]=0\ .$ (22)

Since the initial and final Hamiltonians $H(0)$ and $H(T)$ are determined by
the choice of problem at hand, this commutativity must be ensured only by
imposing boundary conditions on the invariant ${\cal I}(t)$, and not $H(t)$.
Since the invariant ${\cal I}(t)$ is parametrised in terms of $\vec{L}$ and
$R$ via Eqs. (11), (3.1) and (18), this implies boundary conditions for
$\vec{L}$ and $R$.

The choices

$\begin{array}[]{rclcrcl}R(t)&=&M(t)^{-\frac{1}{4}}\ ,&&\dot{R}(t)&=&0\ ,\\\
\vec{L}(t)&=&(mM(t))^{-1}\vec{F}(t)\
,&\mbox{and}&\dot{\vec{L}}(t)&=&0\end{array}$ (23)

directly result in commutativity of ${\cal I}(t)$ and $H(t)$. Imposing the
conditions of Eq. (23) for $t=0$ and $t=T$ thus imposes the desired boundary
conditions Eq. (22).

An additional set of conditions arise from the requirement that ${\cal I}$ be
an invariant. The equations

$\displaystyle\\{\ddot{R},R\\}+\\{R^{2},M\\}$ $\displaystyle=$ $\displaystyle
2[\dot{R},R]_{A}-2RA^{2}R$ $\displaystyle\ddot{\vec{L}}+M\vec{L}$
$\displaystyle=$ $\displaystyle\frac{\vec{F}}{m},$ (24)

which are derived in Sec. 5 as Eqs. (50) and (61), must be satisfied at all
times. Substituting Eqs. (23) into these equations gives an additional set of
boundary conditions

$\ddot{\vec{L}}(t)=0\ ,\ \mbox{ and }\ \ddot{R}(t)=0.$ (25)

For any time other than $t=0$ and $t=T$, the conditions of Eqs.(23) and (25)
do not need to be satisfied, and any time-dependent form of $\vec{L}$ and $R$
that satisfies the boundary conditions may be employed.

In practice, however, some choices of $\vec{L}(t)$ and $R(t)$ can result in
trajectories $X(t)$ or covariances $\Sigma(t)$ of the quantum state that are
conflicting with practical requirements, such as a trapped ion colliding with
a trap electrode. It is thus desirable to relate the variables $\vec{L}(t)$
and $R(t)$ characterizing the invariant ${\cal I}$ to the variables $X(t)$ and
$\Sigma(t)$ characterizing the quantum state.

Since the invariant ${\cal I}(t)$ is quadratic, this is most elegantly done
via the fact that the time-dependent quantum state is the instantaneous ground
state of ${\cal I}(t)$, which is characterized by covariance matrix

$\Sigma_{g}=\frac{1}{2}\Gamma^{-1}$ (26)

and phase-space coordinate

$X_{g}=-\Gamma^{-1}\vec{W}\ .$ (27)

Since $\Gamma$ (given in Eq. (3.1)) is symplectic, its inverse can be
expressed as $\Gamma^{-1}=-{\cal S}\Gamma{\cal S}$, resulting in the explicit
form

$\begin{array}[]{l}\Sigma_{g}=\\\
\displaystyle\frac{1}{2}\left[\begin{array}[]{cc}\displaystyle\frac{R^{2}}{m}&\displaystyle\frac{\\{R,\dot{R}\\}-{{\cal
J}}}{2}\\\\[6.0pt] \displaystyle\frac{\\{R,\dot{R}\\}+{{\cal
J}}}{2}&m\left(\dot{R}^{2}+\Re([\dot{R},R]_{A}-RA^{2}R)\right)\end{array}\right]\end{array}$

for the covariance matrix.

With the form of $\vec{W}$ given in Eq. (18), the explicit solution (Eq. (27))
for $X_{g}$ reads

$X_{g}=\left[\begin{array}[]{c}\vec{L}\\\ m\dot{\vec{L}}\end{array}\right]\ .$
(28)

That is, $X_{g}$ is determined straightforwardly in terms of $\vec{L}$, which
means that the trajectory can be easily chosen at will.

It is of course desirable to attribute an equally clear physical
interpretation to the quantity $R$ as found for $\vec{L}$ in the sense of the
trajectory. This, however, seems feasible only in the adiabatic limit, in
which the time-derivative $\dot{R}$ becomes negligible. ${\cal J}$ in Eq. (17)
thus also becomes negligible, and Eq. (20) reduces to

$R=M^{-\frac{1}{4}}\ .$ (29)

$R$ thus directly determines the curvature of the trapping potential. Apart
from this limiting case, however, the relation between $R$ and the potential
in the system Hamiltonian is substantially more complicated, consistent with
the fact that diabatic dynamics is generally more complicated than adiabatic
dynamics.

##  4 Shuttling protocol

With the general framework for optimal control based on high-dimensional
invariants established, its practical use can be demonstrated with the problem
of shuttling an ion through a two-dimensional potential landscape, such as the
transfer from an initial trapping potential

$V_{i}=\frac{1}{2}m\left(\omega_{t}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}{2}}}x^{2}+\omega_{r}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}{2}}}(y-r)^{2}\right)$
(30)

centered around the initial position $[0,r]$ to a final potential

$V_{f}=\frac{1}{2}m\left(\omega_{r}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}{2}}}(x-r)^{2}+\omega_{t}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}{2}}}y^{2}\right)$
(31)

centered around the final position $[r,0]$. In addition to the change in
position, in particular, also the shape of the potential will change during
the shuttling process with the (typically strong) radial confinement
(characterized by the frequency $\omega_{r}$) initially in the $y$-direction
and finally in the $x$-direction, and the (typically weak) axial confinement
(characterized by the frequency $\omega_{t}$) initially in the $x$-direction
and finally in the $y$-direction.

As discussed above in the context of Eq. (28), the trajectory of the ion can
be chosen at will, as long as it satisfies the required boundary conditions.
We are thus free to chose an arc

$\vec{L}(t)=r\left[\begin{array}[]{c}\displaystyle\sin\left(\frac{\pi}{2}p(\tau)\right)\vspace{.2cm}\\\
\displaystyle\cos\left(\frac{\pi}{2}p(\tau)\right)\end{array}\right]$ (32)

with $\tau=t/T$, radius $r$, and a function $p(\tau)$ satisfying the boundary
conditions $p(0)=0$ and $p(1)=1$. Motivated by its successful use in one-
dimensional problems [Torrontegui  _et al._(2011)Torrontegui, Ibáñez, Chen,
Ruschhaupt, Guéry-Odelin, and Muga, Palmero  _et al._(2013)Palmero,
Torrontegui, Guéry-Odelin, and Muga, Lu  _et al._(2014)Lu, Muga, Chen,
Poschinger, Schmidt-Kaler, and Ruschhaupt] we will choose

$p(\tau)=10\tau^{3}-15\tau^{4}+6\tau^{5}\ ,$ (33)

which ensures that the boundary conditions

$\dot{\vec{L}}(0)=\dot{\vec{L}}(T)=\ddot{\vec{L}}(0)=\ddot{\vec{L}}(T)=0$ (34)

are satisfied.

For the quadratic component $R(t)$ of the invariant, the boundary conditions
read

$R(0)=\left[\begin{array}[]{cc}\sqrt{\omega_{r}}&0\\\ 0&\sqrt{\omega_{t}}\\\
\end{array}\right]=:R_{i}\ ,$ (35)

and

$R(T)=\left[\begin{array}[]{cc}\sqrt{\omega_{t}}&0\\\ 0&\sqrt{\omega_{r}}\\\
\end{array}\right]=:R_{f}\ ,$ (36)

as well as $\dot{R}=0$ and $\ddot{R}=0$ at $t=0$ and $t=T$.

A suitable time-dependent matrix $R(t)$ that satisfies the boundary conditions
can be defined with the ansatz

$R(t)=(1-p(\tau))R_{i}+p(\tau)R_{f}+\tau^{3}\left(1-\tau\right)^{3}R_{c}\ ,$
(37)

The matrix $R_{c}$ that will determine how the trapping potential rotates
along the trajectory of the ion can still be chosen at will, and the following
examples are based on the choice

$R_{c}=(\omega_{t}\omega_{r})^{\frac{1}{4}}\left[\begin{array}[]{cc}1&1\\\
1&1\\\ \end{array}\right]\ .$ (38)

With the given specific choices for $R(t)$ and $\vec{L}(t)$ one can solve
Eqs.(12) and (13) for $\Omega$ and $\vec{V}$, which yields the desired objects
$M$ and $\vec{F}$ via Eqs.(15) and (16).

The explicit time-dependent solution for $M$ depends on the system parameters
$\omega_{r}$ and $\omega_{t}$ and on the duration $T$ of the shuttling
protocol. Due to scale-invariance, the shuttling protocol depends only on two
independent parameters, and in the following, we will consider the variation
of $\frac{\omega_{r}}{\omega_{t}}$ and $\omega_{t}T$.

The explicit time-dependent solution for $\vec{F}$ depends additionally also
on the mass $m$ of the ion and the radius $r$ of the trajectory. Discussing
the solution of the time-dependent trapping potential in terms of its center
$\vec{C}=\frac{1}{m}M^{-1}\vec{F}$ instead of the force $\vec{F}$ yields a
solution that is independent of $m$, and that simply scales linearly with $r$.

As such, one can exemplify shuttling protocols for all values of $\omega_{r}$,
$\omega_{t}$, $m$, $r$ and $T$ in terms of the two independent parameters
$\frac{\omega_{r}}{\omega_{t}}$ and $\omega_{t}T$.

Fig. 1 depicts some explicit solutions for time-dependent trapping potentials
resulting in ground-state to ground-state shuttling protocols. Insets (a) and
(d) correspond to a very fast shuttling protocol with $\omega_{t}T=3$, whereas
insets (b) and (e) depict instances of a slower shuttling protocol with
$\omega_{t}T=5$, and insets (c) and (f) corresponds to a shuttling protocol
with $\omega_{t}T=10$ that is close to adiabatic. The insets (a)-(c)
correspond to a highly anisotropic trap with $\omega_{r}=10\omega_{t}$,
whereas insets (d)-(f) correspond to a less anisotropic trap with
$\omega_{r}=2\omega_{t}$.

(a)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=10\ $, $\ \omega_{t}T=3$

(b)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=10\ $, $\omega_{t}T=5$

(c)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=10\ $, $\omega_{t}T=10$

(d)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=2\ $, $\omega_{t}T=3$

(e)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=2\ $, $\omega_{t}T=5$

(f)1.5cm1.3cm

$\displaystyle\frac{\omega_{r}}{\omega_{t}}=2\ $, $\omega_{t}T=10$

Figure 1:  Trajectories for center of trapping potential (blue) and trap
frequencies (yellow ellipses) for shuttling protocols with increasing duration
from the left insets to the right insets. The purple line denotes the
trajectory of the ion, which traces out the same path in each instance. Fast
protocols (left) result in substantial deviations between the trajectories of
trap center and ion, and this effect is even more pronounced for the more
isotropic trapping potentials.

The arc-shaped trajectory $\vec{L}(t)$ of the ion is depicted in purple, and
the center of the trapping potential is depicted in blue. In the close-to-
adiabatic cases (c) and (f) the trajectories of ion and center of trapping
potential nearly perfectly coincide, which is consistent with the fact that
the ion remains in the instantaneous ground state of the trapping potential in
the adiabatic limit. Deviations between the two trajectories are clearly
discernible for the faster shuttling protocols. In insets (b) and (e) the ion
trajectory remains close to the trajectory of the trapping potential,
indicating that this protocol is close to the regime of validity of the
adiabatic approximation. In insets (a) and (d), however, the trajectory of the
center of trapping potential shares no similarity with the arc-shaped
trajectory of the ion. The ion will thus be substantially displaced from the
center of trapping potential, but the carefully chosen dynamics of the trap
center ensures that the ion comes to rest at the final instance of the
protocol.

Comparison between the instances of the strongly anisotropic trapping
potentials featured in (a)-(c) with the instances of less anisotropic trapping
potentials (d)-(f) indicate that the isotropy has only a minor impact on the
dynamics of the trapping potential in the slow protocols, but that anisotropy
tends to alter dramatically the trajectory of the center of the trapping
potential.

The dynamics of the quadratic component of the trapping potential is
represented by yellow ellipses indicating equipotential curves, and these
ellipses are depicted for $t=jT/8$ with $j=0,\ldots,8$. Even in the fast
shuttling protocols, the shape of the trapping potential rotates without
discernible rapid or non-monotonic dynamics, which is in quite some contrast
to the seemingly wild trajectories of the center of trapping potential.

In order to gauge the experimental resources required for the present
shuttling protocol, it is instructive to specify the maximum of the electric
field $\max_{t\in[0,T]}\lVert\vec{F}(t)/e\rVert$ that is required to keep the
ion on its trajectory. For an Ytterbium-171 ion, with turning radius
$r=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ a trap frequency
$\omega_{t}=2\pi\times$1\text{\,}\mathrm{MHz}$$ this results in an electric
field strength

$\lVert\vec{E}\rVert=2.10\times
10^{5}$\text{\,}\mathrm{V}$$\text{\,}\mathrm{m}$^{-1}$ (39)

for the cases depicted in insets (a) to (c) of Fig. 1, and an electric field
strength

$\lVert\vec{E}\rVert=8.36\times
10^{3}$\text{\,}\mathrm{V}$$\text{\,}\mathrm{m}$^{-1}$ (40)

for the cases depicted in insets (d) to (f). The electric field strength
required in these examples thus depend only on the isotropy of the trapping
potential, but not on the duration of the shuttling protocol. This seems
surprising, but also encouraging since it suggests that increasing the speed
of shuttling does not require increasingly strong control fields.

##  5 Construction of the invariant

Since the derivation of the invariant introduced in Sec. 3.1 is rather
involved, it is deferred to this Section, with some technical steps discussed
in Appendices A and B.

The goal is not to find an invariant for a given Hamiltonian; instead, the
time-dependent invariant is the starting point, and the quest is for a
Hamiltonian such that Eqs. (12), (13) and (14) are satisfied for the given
invariant ${\cal I}(t)$. This can be done in subsequent steps, first
determining the quadratic component of the invariant, parametrized in terms of
$\Gamma$, and subsequently determining the linear component of the invariant,
parametrized in terms of $\vec{W}$ in Eq. (11).

###  5.1 A suitable invariant – the quadratic component

As stated in Eq. (15), $\Omega$ must take the form

$\Omega=\left[\begin{array}[]{cc}mM&\mathbbm{O}\\\
\mathbbm{O}&\frac{1}{m}\mathbbm{1}\end{array}\right]\ ,$ (41)

where $M$ is a real symmetric matrix. The lower right sub block
$\frac{1}{m}\mathbbm{1}$ is associated with the kinetic energy while the term
$M$ is associated with the potential energy, and is the unknown to be deduced
from any proposed invariant.

For a general matrix $\Gamma$, the resultant quadratic part $\Omega$ of the
system Hamiltonian will not be of the desired form given in Eq. (15), but the
parametrization

$\Gamma=\Re\left[\begin{array}[]{cc}m\dot{P}^{\dagger}\dot{P}&-\dot{P}^{\dagger}P\\\
-P^{\dagger}\dot{P}&\frac{1}{m}P^{\dagger}P\end{array}\right]\ ,$ (42)

with four sub-blocks defined in terms of a $d\times d$-dimensional, complex
matrix $P$ will result in the desired form. This can be seen by expressing the
equation of motion for $\Gamma$ (Eq. (12)) in terms of $P$, resulting in

$\left[\begin{array}[]{cc}m\Re(\dot{P}^{\dagger}{\cal D}+{\cal
D}^{\dagger}\dot{P})&-\Re({\cal D}^{\dagger}P)\\\ -\Re(P^{\dagger}{\cal
D})&0\end{array}\right]=0\ $ (43)

where ${\cal D}=\ddot{P}+PM$. Any choice of $P$ satisfying the differential
equation

$\ddot{P}+PM=0\,$ (44)

thus results in a quadratic term $\Gamma$ for the invariant ${\cal I}$ that is
consistent with the initial ansatz in Eq. (15).

With the current parametrisation of $\Gamma$ in terms of $P$, one could
construct $M$ as

$M=-P^{-1}\ddot{P}\ .$ (45)

Usually, however, this does not result in the identification of a suitable
matrix $M$, because $M$ defined in this fashion would not necessarily be
symmetric or real. It is thus necessary to find additional restrictions on $P$
that guarantee the symmetry and reality of $M$.

To this end, it is helpful to express $P$ in terms of its polar decomposition

$P=UR\ ,$ (46)

where $U$ is unitary and $R$ is Hermitian and positive semi-definite. It is
possible to leave the factor $R$ as a free parameter and determine $U$ through
the requirement that $M$ be real-symmetric.

In practice, this task can be formulated in terms of the anti-Hermitian
operator

$A=U^{\dagger}\dot{U}\ .$ (47)

Since $A$ satisfies the differential equation $\dot{U}=UA$, the unitary factor
$U$ is uniquely determined in terms of $A(t)$ and the initial condition
$U(0)=\mathbbm{1}$.

As derived explicitly in Sec.A, requiring $M$ to be Hermitian requires $A$ to
be of the form

$A=iR^{-2}+\frac{1}{2}[R^{-1},\dot{R}]+\frac{1}{2}R^{-1}{{\cal J}}R^{-1}\ ,$
(48)

with

${{\cal J}}=\int_{0}^{t}d\tau\ [M(\tau),R(\tau)^{2}]\ .$ (49)

With this specification of $A$ in terms of $\dot{R}$, $R$ and $M$, one obtains
the differential equation

$\\{\ddot{R},R\\}+\\{R^{2},M\\}=2[\dot{R},R]_{A}-2RA^{2}R$ (50)

for $R$.

So far, it is only imposed that $M$ ought to be Hermitian, and not that it
should be real as well. As shown in Sec. B, however, $M$ is indeed real-
symmetric for all times, if $R$ is chosen such that it is real and positive at
all times and that

$\dot{R}(0)=0\ .$ (51)

Even though Eq.(50), together with Eqs.(48) and (49) may be used to determine
$M(t)$, they are unwieldy in practice, since they relate $M(t)$ to the
integral ${{\cal J}}$ that depends on the entire history of $M$. Since $M$ is
real-symmetric, however, one can – as discussed in detail in Appendix A
arriving at Eq. (78) – specify ${{\cal J}}$ explicitly as

${{\cal
J}}=\sum_{ij}\frac{(\lambda_{i}-\lambda_{j})(\lambda_{i}+\lambda_{j})^{2}}{\lambda_{i}^{2}+\lambda_{j}^{2}}\langle\Phi_{i}|\dot{R}|\Phi_{j}\rangle\ket{\Phi_{i}}\bra{\Phi_{j}}$
(52)

in terms of the eigenvalues $\lambda_{i}$ and eigenvectors $\ket{\Phi_{i}}$ of
$R$. With this explicit form of ${{\cal J}}$, in terms of $R$ and $\dot{R}$,
but without dependence on $M$, one can finally solve Eq. (50) for $M$, which
yields

$M=\sum_{ij}\frac{\langle\Phi_{i}|B|\Phi_{j}\rangle}{\lambda_{i}^{2}+\lambda_{j}^{2}}\ket{\Phi_{i}}\bra{\Phi_{j}}\
,$ (53)

with

$B=2[\dot{R},R]_{A}-2RA^{2}R-\\{\ddot{R},R\\}\ .$ (54)

Any real positive matrix $R$ that satisfies the initial condition
$\dot{R}(0)=0$ (Eq. (51)) thus yields the quadratic term $\Gamma$ of the
invariant via Eq. (42) that is consistent with the desired form for $\Omega$
specified in Eq. (15) with a real-symmetric potential matrix $M$.

Substituting for $P=UR$ and $A=U^{\dagger}\dot{U}$ in the Ansatz for $\Gamma$
(Eq. (42)), results in the form of $\Gamma$ that is quoted in Eq. (3.1) in
Sec. 3.1.

###  5.2 A suitable invariant – the linear component

The goal of this section is to show that the ansatz

$\vec{W}=-\Gamma\left[\begin{array}[]{c}\vec{L}\\\
m\dot{\vec{L}}\end{array}\right]\ ,$ (55)

of Eq.(18) results in a linear part $\vec{V}$ of the Hamiltonian that is of
the form

$\vec{V}=\left[\begin{array}[]{c}-\vec{F}\\\ 0\end{array}\right]\ ,$ (56)

specified in Eq. (16).

$\vec{V}$ can be expressed in terms of $\vec{W}$ as

$\vec{V}={{\cal S}}^{-1}\Gamma^{-1}\left(\Omega{{\cal
S}}\vec{W}-\dot{\vec{W}}\right)$ (57)

via the equation of motion for $\vec{W}$ given in Eq. (13). Rather than
expressing the right-hand-side in terms of $\vec{W}$, it will turn out helpful
to express it in terms of $\vec{Z}$ defined via the relation
$\vec{W}=-\Gamma\vec{Z}$. Together with the equation of motion for $\Gamma$
given in Eq. (12), this yields

$\vec{V}={{\cal S}}^{-1}\dot{\vec{Z}}\ -\Omega\vec{Z}.$ (58)

With the explicit form of the quadratic term $\Omega$ of the Hamiltonian
specified in Eq. (15) this results in the explicit form

$\vec{V}=\left[\begin{array}[]{c}-mM\vec{Z}_{x}-\dot{\vec{Z}}_{p}\\\
\dot{\vec{Z}}_{x}-\frac{1}{m}\vec{Z}_{p}\end{array}\right]$ (59)

for $\vec{V}$, where the vector $\vec{Z}$ is decomposed into spatial and
momentum parts

$\vec{Z}=\left[\begin{array}[]{c}\vec{Z}_{x}\\\ \vec{Z}_{p}\end{array}\right]\
.$ (60)

Finally, choosing $\vec{Z}_{p}=m\dot{\vec{Z}}_{x}$ ensures that $\vec{V}$ is
of the form given in Eq. (16) as desired, and the identification of $\vec{L}$
with $\vec{Z}_{x}$ results in the sought-after form of $\vec{W}$ (Eq. (18)).

The resulting equation of motion

$\ddot{\vec{L}}+M\vec{L}=\frac{\vec{F}}{m},$ (61)

permits the determination of $\vec{F}$ in terms of $\vec{L}$ and $R$ (through
$M$).

###  5.3 A suitable invariant – the scalar component

The scalar part $\theta(t)$ of the invariant ${\cal I}$ is determined by Eq.
(14). Strictly speaking, it does not need to be constructed, since the
eigenstates of the invariant that are the relevant objects are independent of
this phase term. For the sake of completeness, however, the solution

$\theta(t)=\frac{1}{2}\vec{Z}^{T}\Gamma\vec{Z}=-\frac{1}{2}\vec{Z}^{T}\vec{W}=\frac{1}{2}\vec{W}^{T}\Gamma^{-1}\vec{W}\
,$ (62)

which satisfies Eq. (14) is stated.

##  6 The one-dimensional case

The invariant derived in Sec. 5 does not seem to share many similarities with
the invariant that are routinely used in the one-dimensional case [Lewis and
Leach(1982)]. The regular one-dimensional framework is, however, naturally
contained as special case in the present framework, as we will show in the
following.

Since all matrices such as $M$ and $R$ reduce to scalars $M_{s}$ and $R_{s}$
in the one-dimensional case, the commutators and generalized commutators
vanish in the one-dimensional case. In particular, the object ${\cal J}$
defined in terms of a commutator in Eq. (49) vanishes, so that Eq. (48)
simplifies to

$A_{s}=iR_{s}^{-2}.$ (63)

Substituting Eq. (63) into Eq. (50) and exploiting the fact that all objects
commute, gives that

$\ddot{R}_{s}R_{s}+M_{s}R_{s}^{2}=R_{s}^{-2}\ .$ (64)

The quadratic part $\Gamma$ defined in Eq. (3.1) simplifies to

$\displaystyle\Gamma_{s}=\left[\begin{array}[]{cc}\displaystyle
m\left(\dot{R}_{s}^{2}+R_{s}^{-2}\right)&\displaystyle-R_{s}\dot{R}_{s}\\\
\displaystyle-R_{s}\dot{R}_{s}&\frac{1}{m}R_{s}^{2}\end{array}\right]\ ,$ (67)

in the case of commutativity. From Eqs. (18), (19) and (67), it is possible to
write the invariant ${\cal I}$ out explicitly in terms of $\hat{r}$ and
$\hat{p}$,

$\displaystyle{\cal I}_{s}$
$\displaystyle=\frac{1}{2m}\left(R_{s}(\hat{p}-m\dot{L})-m\dot{R}_{s}(\hat{r}-L)\right)^{2}$
$\displaystyle+\frac{1}{2}m\left(\frac{\hat{r}-L}{R_{s}}\right)^{2},$ (68)

which is the Ermakov-Lewis invariant [Lewis(1967)] with additional terms
linear in position and momentum [Lewis and Leach(1982)]. Dividing both sides
of Eq. (64) by $R_{s}$ gives

$\ddot{R}_{s}+R_{s}M_{s}=R_{s}^{-3},$ (69)

which is the Ermakov equation that appears in the original treatment
[Lewis(1967)], where the symbol $\rho$ is used instead of $R_{s}$, and the
symbols $\omega^{2}$ and $\alpha$ are used instead of $M_{s}$ and $\vec{L}$.

##  7 Outlook

The extension of the framework of quantum invariants to systems with quadratic
potentials of more than one translational degree of freedom opens up a wide
perspective for the control of quantum systems. While the present work is
motivated by control of trapped ions, none of the features of the present
framework are specific to this system, and the only precondition for the use
of the presented invariant is the validity of the Gaussian approximation.
Since also shuttling of cold neutral atoms has been experimentally realised in
one-dimensional optical lattice potentials [Lu  _et al._(2020)Lu, Ruschhaupt,
Martínez-Garaot, and Muga] with time-dependent potential designed with the
Ermakov-Lewis invariant, the presently derived framework shall open up
possibilities to extend such tasks to higher-dimensional optical lattices.

The ability to systematically devise time-dependent potentials that realise
ground-state to ground-state transfer permits to identify those potentials
that are optimized with respect to additional desirable properties. Typical
examples might include shuttling protocols with particularly weak forcing or
potentials that suit best the geometry of the electrodes that generate the
potentials. As such, the present framework promises to substantially advance
the shuttling of trapped ions, which is a central step towards achieving
scalability of quantum information technology with trapped ions.

The specific choice of invariant discussed here, does not necessarily need to
be seen as the unique generalization of one-dimensional invariants, but rather
as evidence that such generalizations are possible. It thus seems conceivable
that different generalization can be found, and that different restrictions of
the achievable Hamiltonians can be found. One may, for example, envision the
Lorentz force caused by a homogeneous magnetic field and aim at finding an
invariant that ensures that no inhomogeneities in the magnetic field arise.

In addition to the extension of experimental tasks such as noise-resilient
shuttling [Levy  _et al._(2018)Levy, Kiely, Muga, Kosloff, and Torrontegui, Lu
_et al._(2014)Lu, Muga, Chen, Poschinger, Schmidt-Kaler, and Ruschhaupt, Lu
_et al._(2018)Lu, Ruschhaupt, and Muga] to higher dimensional settings, the
present work can also initiate new steps in the conceptual developments of
quantum invariants.

## Acknowledgements

We are indebted to stimulating discussions with Adam Callison, Alexander
Paige, Pedro Taylor-Burdett, Sebastian Weidt and Winnie Hensinger. This work
was supported through a studentship in the Centre for Doctoral Training on
Controlled Quantum Dynamics at Imperial College London funded by
EPSRC(EP/L016524/1).

##  Appendix A Derivation of the invariant equations

This section contains the explicit derivation of Eqs. (48),(49) and (50).

The equation of motion $\ddot{P}+PM=0$ (Eq. (44)) with $P$ expressed in terms
of its polar decomposition $P=UR$ (Eq. (46)) results in

$K:=\dot{A}R+A^{2}R+2A\dot{R}+\ddot{R}+RM=0\ .$ (70)

Since $M$ is Hermitian, and the relation $RK=K^{\dagger}R$ holds because $K$
vanishes, this is equivalent to

$\displaystyle\frac{d}{dt}RAR=$
$\displaystyle\frac{1}{2}\left([\ddot{R},R]+[M,R^{2}]\right)\ .$ (71)

The right hand side of Eq. (71) can be rewritten in terms of a total
derivative to obtain

$\frac{d}{dt}RAR=\frac{1}{2}\frac{d}{dt}\left([\dot{R},R]+{{\cal J}}\right),$
(72)

with

${{\cal J}}=\int^{t}_{0}d\tau\ [M(\tau),R(\tau)^{2}]\ .$ (73)

This can be directly integrated, resulting in the explicit solution

$A=R^{-1}CR^{-1}+\frac{1}{2}[R^{-1},\dot{R}]+\frac{1}{2}R^{-1}{{\cal
J}}R^{-1},$ (74)

for $A$.

The matrix $C$ is the constant of integration that can be chosen freely. With
the choice $C=i\mathbbm{1}$ motivated by the Ermakov equation in the one-
dimensional case [Lewis(1967), Leach and Andriopoulos(2008)], this results in
Eq. (48).

Using again the fact that $K$ vanishes, we can form the relation
$RK+K^{\dagger}R=0$ and rearrange to derive Eq. (50).

##  Appendix B Reality of $M$

This section contains the proof that $M(t)$ determined by Eq.(50) is indeed
real and symmetric for all times, provided that the initial condition
$\dot{R}(0)=0$ is satisfied.

Eq. (50) reads explicitly

$\\{\ddot{R},R\\}+\\{R^{2},M\\}=2[\dot{R},R]_{A}-2RA^{2}R\ ,$ (75)

and there is indeed no evident reason why $M$ determined from this relation
together with Eqs. (48) and (49) should be real and symmetric.

In the following, explicit, closed-form expressions for ${{\cal J}}$, $A$ and
$M$ are derived, that hold under the assumption that $M$ is real and
symmetric. It is then demonstrated that such closed-form expressions form
indeed a solution of Eqs. (48), (49) and (50), which not only completes the
proof that $M$ is real symmetric, but gives a direct way to calculate $M$.

###  B.1 Necessary conditions for the reality of $M$

Under the assumption that $M$ is real and symmetric, the imaginary part of the
equation of motion for $R$ (Eq. (50)) reads

$\\{{{\cal J}},R^{-2}\\}=[\dot{R},R^{-1}]+[R,R^{-2}]_{\dot{R}}\ .$ (76)

This determines ${{\cal J}}$ in terms of $R$ and $\dot{R}$, and, in contrast
to the defining relation (Eq. (49)) for ${{\cal J}}$ it is time-local.

In order to solve Eq. (76) for ${{\cal J}}$, it is helpful to consider the
eigen-decomposition

$R=\sum_{j}\lambda_{j}|\Phi_{j}\rangle\langle\Phi_{j}|$ (77)

of the real and positive semi-definite matrix $R$.
In terms of matrix elements with respect to the eigenstates $\ket{\Phi_{j}}$
of $R$, Eq. (76) can directly be solved, yielding

${{\cal
J}}_{jk}=\frac{(\lambda_{j}-\lambda_{k})(\lambda_{j}+\lambda_{k})^{2}}{\lambda_{j}^{2}+\lambda_{k}^{2}}{\dot{R}}_{jk}\
,$ (78)

where

${O}_{jk}=\langle\Phi_{j}|O|\Phi_{k}\rangle$ (79)

is a short hand notation for the matrix elements of operator $O$, and

$\dot{O}_{jk}=\langle\Phi_{j}|\dot{O}|\Phi_{k}\rangle$ (80)

does generally not coincide with $\frac{\partial}{\partial t}O_{jk}$ because
of the time-dependence of the eigenstates $\ket{\Phi_{j}}$.

Similarly to Eq.(78), also Eq. (50) can be solved for $M_{jk}$, resulting in

${M}_{jk}=\frac{1}{\lambda_{j}^{2}+\lambda_{k}^{2}}\left(2\sum_{l}\left(\lambda_{k}{\dot{R}}_{jl}A_{lk}-\lambda_{j}A_{jl}\dot{R}_{lk}-\lambda_{j}\lambda_{k}A_{jl}A_{lk}\right)-{\ddot{R}}_{jk}\left(\lambda_{j}+\lambda_{k}\right)\right)\
.$ (81)

Eq. (48) expressed similarly in terms of matrix elements reads

${A}_{jk}=\frac{1}{2\lambda_{j}\lambda_{k}}\left(2i\delta_{jk}-(\lambda_{j}-\lambda_{k}){\dot{R}}_{jk}+{{\cal
J}}_{jk}\right)\ .$ (82)

Substituting Eq. (78) into Eq. (82) yields

$A_{jk}=\frac{i}{\lambda_{j}\lambda_{k}}\delta_{jk}+\frac{\lambda_{j}-\lambda_{k}}{\lambda_{j}^{2}+\lambda_{k}^{2}}\dot{R}_{jk}.$
(83)

Substituting Eq. (83) into Eq. (81) yields

${M}_{jk}=\frac{\delta_{jk}}{\lambda_{j}^{4}}-\frac{(\lambda_{j}+\lambda_{k})}{\lambda_{j}^{2}+\lambda_{k}^{2}}{\ddot{R}}_{jk}+2\sum_{l}{\dot{R}}_{jl}{\dot{R}}_{lk}p(\lambda_{j},\lambda_{k},\lambda_{l})\
,$ (84)

with

$p(\lambda_{j},\lambda_{k},\lambda_{l})=\frac{\lambda_{l}^{3}\left(\lambda_{j}+\lambda_{k}\right)-\lambda_{l}^{2}\left(\lambda_{j}^{2}+\lambda_{k}^{2}-\lambda_{j}\lambda_{k}\right)-\lambda_{j}^{2}\lambda_{k}^{2}}{(\lambda_{j}^{2}+\lambda_{k}^{2})(\lambda_{j}^{2}+\lambda_{l}^{2})(\lambda_{k}^{2}+\lambda_{l}^{2})}.$
(85)

The closed-form expressions for ${{\cal J}}$, $A$ and ultimately $M$, in terms
of their matrix elements, given in Eqs, (78), (83) and (84), hold if $M$ is
real symmetric.

Conversely, if Eqs. (78), (83) and (84) hold, then $M$ is guaranteed to be
real symmetric since Eq. (84) is manifestly symmetric and defined solely using
real quantities. As a consequence, $M$ is real symmetric if and only if these
expressions hold.

###  B.2 Direct proof of the reality of $M$

What remains is thus to verify directly that taking Eqs. (78), (83) and (84)
as defining relations for ${\cal J}$, $A$ and $M$, implies that Eqs. (48),
(49) and (50) are satisfied. Doing this completes the proof that $M$ is real
symmetric.

Henceforth it is assumed that Eqs. (78), (83) and (84) hold. As Eq. (83) is
derived using Eq. (82) which is Eq. (48) expressed in terms of the matrix
elements of the eigenstates of $R$, Eq. (48) holds by construction. Similarly,
Eq. (84) is derived from Eq. (81), which itself is Eq. (50) expressed in terms
of matrix elements, so Eq. (50) additionally holds.

It remains to show that Eq. (49), or its equivalent differential form

$\displaystyle{{\cal J}}(0)$ $\displaystyle=0\ ,\ \mbox{and}$ (86)
$\displaystyle\dot{{\cal J}}$ $\displaystyle=R^{2}M-MR^{2}\ .$ (87)

is satisfied Since $\dot{R}(0)=0$, Eq. (86) follows directly from inspection
of Eq. (78).

Eq.(87) reads in terms of matrix elements

$\dot{{\cal J}}_{jk}=(\lambda_{j}^{2}-\lambda_{k}^{2})M_{jk}.$ (88)

In order to complete the proof, it is sufficient to show that Eq. (88) is
satisfied. As ${\cal J}$ is defined in terms of its matrix elements ${\cal
J}_{jk}$ in Eq. (78), it is necessary to relate these to the $\dot{{\cal
J}_{jk}}$ that appear in Eq. (88).

The time-derivative of ${\cal J}_{jk}$ reads

$\frac{d}{dt}{\cal J}_{jk}=\sum_{l}W_{jl}{\cal J}_{lk}-\sum_{l}{\cal
J}_{jl}W_{lk}+\dot{{\cal J}}_{jk}$ (89)

with

$W_{jk}=\langle{{\dot{\Phi}}_{j}}|{\Phi_{k}}\rangle.$ (90)

Substituting this expression into Eq. (88) gives

$\frac{d}{dt}{{{{\cal J}}}}_{jk}=\sum_{l}W_{jl}{{\cal J}}_{lk}-\sum_{l}{{\cal
J}}_{jl}W_{lk}+(\lambda_{j}^{2}-\lambda_{k}^{2}){M}_{jk},$ (91)

which is equivalent to Eq. (88).

The next step is to differentiate Eq. (78) and use the result to show that Eq.
(91) is satisfied identically.

Defining, motivated by Eq. (91),

$\Xi=\frac{d}{dt}{{\cal J}}_{jk}-\sum_{l}W_{jl}{{\cal J}}_{lk}+\sum_{l}{{\cal
J}}_{jl}W_{lk}+(\lambda_{k}^{2}-\lambda_{j}^{2}){M}_{jk},$ (92)

the remaining task is to show that $\Xi=0$.

Substituting for ${{\cal J}}_{jk}$ and ${M}_{jk}$ using Eqs. (78) and (84) in
Eq. (92) gives

$\displaystyle\Xi$
$\displaystyle=\frac{d}{dt}\left(\frac{(\lambda_{j}^{2}-\lambda_{k}^{2})(\lambda_{j}+\lambda_{k})}{\lambda_{j}^{2}+\lambda_{k}^{2}}\right){\dot{R}}_{jk}+\frac{(\lambda_{j}^{2}-\lambda_{k}^{2})(\lambda_{j}+\lambda_{k})}{\lambda_{j}^{2}+\lambda_{k}^{2}}\frac{d}{dt}\left({{\dot{R}}}_{jk}\right)$
$\displaystyle-\sum_{l}\frac{(\lambda_{l}^{2}-\lambda_{k}^{2})(\lambda_{l}+\lambda_{k})}{\lambda_{l}^{2}+\lambda_{k}^{2}}W_{jl}{\dot{R}}_{lk}+\sum_{l}\frac{(\lambda_{j}^{2}-\lambda_{l}^{2})(\lambda_{j}+\lambda_{l})}{\lambda_{j}^{2}+\lambda_{l}^{2}}{\dot{R}}_{jl}W_{lk}$
$\displaystyle+\frac{(\lambda_{j}+\lambda_{k})(\lambda_{k}^{2}-\lambda_{j}^{2})}{\lambda_{j}^{2}+\lambda_{k}^{2}}{\ddot{R}}_{jk}-2\left(\lambda_{k}^{2}-\lambda_{j}^{2}\right)p(\lambda_{j},\lambda_{k},\lambda_{l}){\dot{R}}_{jl}{\dot{R}}_{lk}\
.$ (93)

Given that the states $\ket{\Phi_{i}}$ are the eigenvectors of $R$, the
explicit expressions for $R_{ij}$, $\dot{R}_{ij}$ and $\ddot{R}_{ij}$ read

$\displaystyle R_{jk}$ $\displaystyle=\lambda_{j}\delta_{jk}\ ,$ (94)
$\displaystyle\dot{R}_{jk}$
$\displaystyle=(\lambda_{j}-\lambda_{k})W_{jk}+\dot{\lambda}_{j}\delta_{jk}\
,$ (95) $\displaystyle{\ddot{R}}_{jk}$
$\displaystyle=\sum_{l}(\lambda_{j}+\lambda_{k}-2\lambda_{l})W_{jl}W_{lk}+(\lambda_{j}-\lambda_{k})\dot{W}_{jk}+2(\dot{\lambda}_{j}-\dot{\lambda}_{k})W_{jk}+\ddot{\lambda}_{j}\delta_{jk}\
.$ (96)

Substituting Eq. (95) and Eq. (96) into Eq. (B.2) gives

$\Xi=W_{jk}r(\lambda_{j},\lambda_{k},\dot{\lambda}_{j},\dot{\lambda}_{k})+\sum_{l}W_{jl}W_{lk}s(\lambda_{j},\lambda_{k},\lambda_{l}),$
(97)

where

$\displaystyle r(\lambda_{j},\lambda_{k},\dot{\lambda}_{j},\dot{\lambda}_{k})$
$\displaystyle=(\lambda_{j}-\lambda_{k})\frac{d}{dt}\left(\frac{(\lambda_{j}^{2}-\lambda_{k}^{2})(\lambda_{j}+\lambda_{k})}{\lambda_{j}^{2}+\lambda_{k}^{2}}\right)-\frac{(\lambda_{j}^{2}-\lambda_{k}^{2})(\lambda_{j}+\lambda_{k})(\dot{\lambda}_{j}-\dot{\lambda}_{k})}{\lambda_{j}^{2}+\lambda_{k}^{2}}$
$\displaystyle-2(\lambda_{j}-\lambda_{k})(\lambda_{k}^{2}-\lambda_{j}^{2})\left(\dot{\lambda}_{j}p(\lambda_{j},\lambda_{k},\lambda_{j})+\dot{\lambda}_{k}p(\lambda_{j},\lambda_{k},\lambda_{k})\right)\
,$ (98)

and

$\displaystyle s(\lambda_{j},\lambda_{k},\lambda_{l})$
$\displaystyle=-\frac{(\lambda_{l}^{2}-\lambda_{k}^{2})^{2}}{\lambda_{l}^{2}+\lambda_{k}^{2}}+\frac{(\lambda_{j}^{2}-\lambda_{l}^{2})^{2}}{\lambda_{j}^{2}+\lambda_{l}^{2}}+\frac{(\lambda_{j}+\lambda_{k})(\lambda_{k}^{2}-\lambda_{j}^{2})(\lambda_{j}+\lambda_{k}-2\lambda_{l})}{\lambda_{j}^{2}+\lambda_{k}^{2}}$
$\displaystyle-2(\lambda_{j}-\lambda_{l})(\lambda_{k}-\lambda_{l})(\lambda_{j}^{2}-\lambda_{k}^{2})p(\lambda_{j},\lambda_{k},\lambda_{l})\
.$ (99)

Both $r$ and $s$ are rational functions of their arguments, and they are
indeed identically vanishing. Therefore $\Xi=0$ which completes the proof.

## References

  * [Friis  _et al._(2018)Friis, Marty, Maier, Hempel, Holzäpfel, Jurcevic, Plenio, Huber, Roos, Blatt, and Lanyon] N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt, and B. Lanyon, Phys. Rev. X 8, 021012 (2018).
  * [Wei  _et al._(2020)Wei, Lauer, Srinivasan, Sundaresan, McClure, Toyli, McKay, Gambetta, and Sheldon] K. X. Wei, I. Lauer, S. Srinivasan, N. Sundaresan, D. T. McClure, D. Toyli, D. C. McKay, J. M. Gambetta, and S. Sheldon, Phys. Rev. A 101, 032343 (2020).
  * [Monroe  _et al._(2014)Monroe, Raussendorf, Ruthven, Brown, Maunz, Duan, and Kim] C. Monroe, R. Raussendorf, A. Ruthven, K. R. Brown, P. Maunz, L.-M. Duan, and J. Kim, Phys. Rev. A 89, 022317 (2014).
  * [Lekitsch  _et al._(2017)Lekitsch, Weidt, Fowler, Mølmer, Devitt, Wunderlich, and Hensinger] B. Lekitsch, S. Weidt, A. G. Fowler, K. Mølmer, S. J. Devitt, C. Wunderlich, and W. K. Hensinger, Science Advances 3(2), e1601540 (2017).
  * [Kielpinski  _et al._(2002)Kielpinski, Monroe, and Wineland] D. Kielpinski, C. Monroe, and D. J. Wineland, Nature 417, 709 (2002), iSBN: 0028-0836.
  * [Mølmer and Sørensen(1999)] K. Mølmer and A. Sørensen, Phys. Rev. Lett. 82, 1835 (1999).
  * [Rowe  _et al._(2002)Rowe, Ben-Kish, Demarco, Leibfried, Meyer, Beall, Britton, Hughes, Itano, Jelenković, Langer, Rosenband, and Wineland] M. A. Rowe, A. Ben-Kish, B. Demarco, D. Leibfried, V. Meyer, J. Beall, J. Britton, J. Hughes, W. M. Itano, B. Jelenković, C. Langer, T. Rosenband, and D. J. Wineland, Quant. Inf. Comp. 2(4), 257 (2002).
  * [Hensinger  _et al._(2006)Hensinger, Olmschenk, Stick, Hucul, Yeo, Acton, Deslauriers, Monroe, and Rabchuk] W. K. Hensinger, S. Olmschenk, D. Stick, D. Hucul, M. Yeo, M. Acton, L. Deslauriers, C. Monroe, and J. Rabchuk, Appl. Phys. Lett. 88, 034101 (2006).
  * [Leibfried  _et al._(2003)Leibfried, Blatt, Monroe, and Wineland] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  * [Heller(1975a)] E. J. Heller, J. Chem. Phys 62, 1544 (1975a).
  * [Heller(1975b)] E. J. Heller, Chem. Phys. Lett. 34(2), 321 (1975b).
  * [Heller(1976)] E. J. Heller, J. Chem. Phys 65, 4979 (1976).
  * [Heller(1977)] E. J. Heller, J. Chem. Phys 66, 5777 (1977).
  * [Heller(1981)] E. J. Heller, J. Chem. Phys 75, 2923 (1981).
  * [Lewis(1968)] H. R. Lewis, Phys. Rev. 172, 1313 (1968).
  * [Chen  _et al._(2011)Chen, Torrontegui, and Muga] X. Chen, E. Torrontegui, and J. G. Muga, Phys. Rev. A 83, 062116 (2011).
  * [Torrontegui  _et al._(2011)Torrontegui, Ibáñez, Chen, Ruschhaupt, Guéry-Odelin, and Muga] E. Torrontegui, S. Ibáñez, X. Chen, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. A 83, 013415 (2011).
  * [Levy  _et al._(2018)Levy, Kiely, Muga, Kosloff, and Torrontegui] A. Levy, A. Kiely, J. G. Muga, R. Kosloff, and E. Torrontegui, New J. Phys. 20, 025006 (2018).
  * [Torrontegui  _et al._(2013)Torrontegui, Ibáñez, Martínez-Garaot, Modugno, del Campo, Guéry-Odelin, Ruschhaupt, Chen, and Muga] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Adv. At. Mol. Opt. Phys. 62, 117 (2013).
  * [Tobalina  _et al._(2017)Tobalina, Palmero, Martínez-Garaot, and Muga] A. Tobalina, M. Palmero, S. Martínez-Garaot, and J. G. Muga, Sci. Rep. 7, 5753 (2017).
  * [Lu  _et al._(2020)Lu, Ruschhaupt, Martínez-Garaot, and Muga] X.-J. Lu, A. Ruschhaupt, S. Martínez-Garaot, and J. G. Muga, Entropy 22(3), 262 (2020).
  * [Zhang  _et al._(2016)Zhang, Muga, Guéry-Odelin, and Chen] Q. Zhang, J. G. Muga, D. Guéry-Odelin, and X. Chen, J. Phys. B 49(12), 125503 (2016).
  * [Ness  _et al._(2018)Ness, Shkedrov, Florshaim, and Sagi] G. Ness, C. Shkedrov, Y. Florshaim, and Y. Sagi, New J. Phys. 20, 095002 (2018).
  * [Chen  _et al._(2018)Chen, Jiang, Li, Ban, and Sherman] X. Chen, R.-L. Jiang, J. Li, Y. Ban, and E. Y. Sherman, Phys. Rev. A 97, 013631 (2018).
  * [Rosu and Romero(1999)] H. C. Rosu and J. L. Romero, Nuovo Cim. B 114, 569 (1999).
  * [Espinoza(2000)] P. B. Espinoza, arXiv:math-ph/0002005 (2000).
  * [Goncharenko  _et al._(1991)Goncharenko, Logvin, Samson, Shapovalov, and Turovets] A. M. Goncharenko, Y. A. Logvin, A. M. Samson, P. S. Shapovalov, and S. I. Turovets, Phys. Lett. A 160(2), 138 (1991).
  * [Tobalina  _et al._(2020)Tobalina, Torrontegui, Lizuain, Palmero, and Muga] A. Tobalina, E. Torrontegui, I. Lizuain, M. Palmero, and J. G. Muga, Phys. Rev. A 102, 063112 (2020).
  * [Mintert and Heller(2009)] F. Mintert and E. J. Heller, EPL 86(5), 50006 (2009).
  * [Lewis(1967)] H. R. Lewis, Phys. Rev. Lett. 18, 510 (1967).
  * [Palmero  _et al._(2013)Palmero, Torrontegui, Guéry-Odelin, and Muga] M. Palmero, E. Torrontegui, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. A 88, 053423 (2013).
  * [Lu  _et al._(2014)Lu, Muga, Chen, Poschinger, Schmidt-Kaler, and Ruschhaupt] X.-J. Lu, J. G. Muga, X. Chen, U. G. Poschinger, F. Schmidt-Kaler, and A. Ruschhaupt, Phys. Rev. A 89, 063414 (2014).
  * [Lewis and Leach(1982)] H. R. Lewis and P. G. L. Leach, J. Math. Phys 23, 2371 (1982).
  * [Lu  _et al._(2018)Lu, Ruschhaupt, and Muga] X.-J. Lu, A. Ruschhaupt, and J. G. Muga, Phys. Rev. A 97, 053402 (2018).
  * [Leach and Andriopoulos(2008)] P. Leach and S. Andriopoulos, Appl. Anal. Discrete Math. 2(2), 146 (2008).