# Controlled Quasi-Latitudinal Solutions for ultra-fast Spin-Torque
Precessional Magnetization Switching

Alessandro Fortunati, Massimiliano d’Aquino and Claudio Serpico

Department of Electrical Engineering and Information Technologies,
University of Naples Federico II, I-80125 Naples, Italy E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The aim of the paper is to present a novel class of time-dependent controls to
realize ultra-fast magnetization switching in nanomagnets driven by spin-
torques produced by spin-polarized electric currents. Magnetization dynamics
in such systems is governed by the Landau-Lifshitz-Slonczewski equation which
describes the precessional motion of (dimensionless) magnetization vector on
the unit-sphere. The relevant case of nanoparticles with uniaxial anisotropy
having in-plane easy and intermediate axes and out-of-plane hard axis is
considered. By exploiting the characteristic smallness of damping and spin-
torque intensity, the aforementioned controls are constructed via suitable
perturbative tools in a way to realise approximate _latitudinal solutions_
(i.e. motions on a sphere in which the out-of-plane magnetization component
stays constant) with the effect to fast “switch” the system from one
stationary state to another. The possibility to keep a (“small”) bounded value
of the out-of-plane coordinate throughout this process of “transfer”, turns
out to be advantageous in the applications as it sensibly reduces the post-
switching relaxation oscillations that may cause the failure of switching in
real samples. Further relevant quantitative results on the behaviour of the
solutions during the pre- and post- switching stages (termed “expulsion” and
“attraction”, respectively), are given as a byproduct. A selection of
validating numerical experiments is presented alongside the corresponding
theoretical results.

Keywords: Magnetisation Dynamics, Landau-Lifshitz-Slonczewski equation,
Spintronics, Perturbation Theory, Qualitative Methods.
2010 MSC. Primary: 78A25, 34D10, 34H05. Secondary: 37C50, 37C75.

Figure 1: Sketch of a MTJ used as STT-MRAM cell.

##  1 Introduction

The efficient and high-speed manipulation of magnetic nanoelements holds
immense significance in the framework of magnetization dynamics, particularly
in the context of magnetic storage nanodevices and spintronics. Over the past
decade, an extensive body of research has been dedicated to investigating
ultra-fast magnetization switching within spintronic devices[DPG+20], which
hold promise as potential candidates for advancing Magnetic Random Access
Memory (MRAM) technology [ADS16, BSH+17].

In the pursuit of achieving rapid magnetization switching, researchers have
turned their attention to the utilization of electric current pulses via the
spin-transfer torque (STT) effect [Slo96]. This approach has shown the
potential to facilitate sub-nanosecond switching in devices like MRAM cells
[BLS+10]. Nevertheless, it is worth noting that ballistic spin-torque
magnetization switching engendered by pulsed injected currents encounters a
challenge reminiscent of the circumstance faced in precessional switching via
transverse external magnetic fields [KR02].

Indeed, both of these switching methodologies demand meticulous
synchronization of the excitation pulse timing to ensure effective switching
occurs at the precise moment[dSS+04]. This critical timing requirement is
essential for avoiding unsuccessful switching and ensuring reliable
performance. Once the excitation pulse is appropriately timed, the
magnetization state transits from a high-energy to a low-energy configuration,
ultimately reaching an equilibrium magnetization state [BMSd03]. Intriguingly,
the underlying relaxation mechanism driving this transition is inherently
stochastic, even in scenarios where thermal fluctuations are relegated to a
secondary role [SdBM09].

This inherent stochasticity can be attributed to the intricate interplay
between the extreme sensitivity of the system to initial conditions and its
multistable nature [BSM13]. These factors collectively bestow upon the system
a probabilistic dimension, contributing to the inherent randomness in the
relaxation process. This intrinsic unpredictability presents both challenges
and opportunities for controlling and optimizing spin torque-induced
magnetization switching in nanoscale devices.

In this respect, there is great interest in developing strategies to achieve
fast, reliable and energy-efficient spin-torque switching. It has been shown
that the reliability of spin-torque ballistic switching in in-plane
ferromagnetic nanodots can be strongly improved by controlling the quasi-
randomness by using appropriate bias fields [dPS+15b]. On the other hand,
optimization of current pulse design targeted at minimization of the energy
cost has been recently proposed [VKL+22].

In this paper, we propose a full analytical treatment for precessional spin-
torque switching of in-plane magnetized nanomagnets that relies on the idea of
forcing magnetization to evolve on controlled quasi-latitudinal (CQL)
trajectories associated with lower values of the ferromagnet’s free energy,
which drive magnetization very close to the target reversed state. This is
achieved by exciting magnetization precession using injected current pulses of
suitable shape. This switching scheme minimizes ringing phenomena that are the
main source for failure of the switching process due to back-hopping into the
original magnetization state [MSB+09, DBB+20]. By developing suitable
perturbation theory, the rigorous conditions for the realization of CQL
magnetization switching as well as the analytical expression for the optimal
current pulse shape are derived as function of the system physical parameters
and summarised in dedicated Lemmata.

The paper is organized as follows. In section 2, a simple schematic of the
magnetic nanosystem studied as the archetypal for in-plane STT-MRAM cells is
introduced and the equation governing the magnetization dynamics driven by the
spin-transfer torque is recalled along with the relevant parameters of the
system. Then, after some preliminaries in sec. 3, the main result of the paper
is stated in sec. 4 as Theorem 4.1. The development of this result requires
the definition and the detailed study of three stages of the spin-torque
precessional switching, termed expulsion, transfer and attraction, which is
carried out in the sections 5-7. Finally, validation of the proposed approach
that demonstrated its effectiveness is presented.

##  2 Magnetization dynamics driven by Spin-Transfer Torque

We consider the magnetic nanosystem sketched in fig.1, representing a magnetic
tunnel junction (MTJ) sandwiched between two electrodes and subject to a
current of intensity $i$ flowing through the MTJ in the $z$ direction. A
typical MTJ structure [HYN+20] is composed by two ferromagnetic layers
separated by a non-magnetic (NM) (insulating) layer. The magnetization of one
layer, termed fixed (lower layer labelled as FM in fig.1), is artificially
pinned to a given orientation and acts as polarizer for electron spins flowing
through it. The second ferromagnetic layer, termed free (upper FM layer in
fig.1), is where magnetization dynamics can take place driven by the external
actions (injected current and applied magnetic field). The MTJ acts as a
single bit-cell of a MRAM cell array where the bit state $0,1$ is coded into
the mutual orientation of the free and fixed layers, namely parallel (P) or
anti-parallel (AP) magnetization, respectively. The switching of the bit is
triggered by the spin-transfer (STT) torque created by the electric current
produced when both the bit and word lines are simultaneously addressed [ADS16,
BSH+17], which is able to eventually switch the magnetization in the free
layer of the MTJ.

The mechanism governing the spin-torque switching can be conveniently
explained and quantified under the assumption of spatially-uniform
magnetization and negligible thermal fluctuations. In this situation,
magnetization dynamics is described by the Landau-Lifshitz-Slonczewski (LLS)
equation (written in dimensionless form)[Slo96, MBS09]:

$\frac{d\bm{m}}{dt}=-\bm{m}\times\bm{h}_{\mathrm{eff}}-\alpha\bm{m}\times(\bm{m}\times\bm{h}_{\mathrm{eff}})+\beta\bm{m}\times(\bm{m}\times\bm{e}_{p})\quad,$
(1)

where $\bm{m}$ is the magnetization unit-vector, $\bm{h}_{\mathrm{eff}}$ is
the effective magnetic field normalized by the saturation magnetization
$M_{s}$, $\alpha$ is the Gilbert damping, $\beta$ is the normalized current
which measures the strength of the spin-transfer torque, $\bm{e}_{p}$ is the
fixed layer unit-vector acting as polarizer for the traversing electric
current, and time is measured in units of $(\gamma M_{s})^{-1}$ ($\gamma$ is
the absolute value of the gyromagnetic ratio). The function $\beta(t)$ is
proportional to the injected current $i(t)$ through the relationship
$\beta(t)=b_{p}i(t)/(SJ_{p})$ where $b_{p}$ is a model-dependent parameter in
the order of unity, $S$ is the device cross sectional area, and
$J_{p}=\mu_{0}M_{s}^{2}|e|d/\hbar$ is a characteristic current density
($\mu_{0}$ is the vacuum permeability, $e$ is the electron charge, $d$ is the
thickness of the free layer, and $\hbar$ is the reduced Planck constant).

The effective field is expressed as:

$\bm{h}_{\mathrm{eff}}=-\frac{\partial
g_{L}}{\partial\bm{m}}=-D_{1}m_{1}\bm{e}_{1}-D_{2}m_{2}\bm{e}_{2}-D_{3}m_{3}\bm{e}_{3}+\bm{h}_{a}\quad,$
(2)

where $\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}$ are the cartesian unit-vectors along
the coordinate axes, $g_{L}$ is the magnetic free energy of the particle

$g_{L}(\bm{m},\bm{h}_{a})=\frac{1}{2}D_{1}m_{1}^{2}+\frac{1}{2}D_{2}m_{2}^{2}+\frac{1}{2}D_{3}m_{3}^{2}-\bm{h}_{a}\cdot\bm{m}\,\,,$
(3)

$\bm{h}_{a}$ is the external applied magnetic field and $D_{1},D_{2},D_{3}$
are effective demagnetizing factors taking into account shape and magneto-
crystalline anisotropy. In the sequel, we assume

$0<D_{1}<D_{2}<D_{3}\mbox{,}$ (4)

meaning that $\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}$ refer to the easy,
intermediate, hard axes, respectively.

In the absence of injected current $\beta=0$ and small enough bias applied
field $\bm{h}_{a}$, magnetization lies in one of the stable equilibria aligned
with the easy axis. We assume that this equilibrium is $s^{-}$ such that
$m_{1}\approx-1$. Under this assumption, it is apparent from eq.(1) that, in
order to have the maximum torque on magnetization when the current is turned
on, the polarizer $\bm{e}_{p}$ must be orthogonal to the easy axis. In this
paper we assume $\bm{e}_{p}=\bm{e}_{3}$ for reasons that will be clear in the
sequel.

##  3 Preliminaries and set-up

Leu us consider the following system of ODEs

$\begin{cases}\dot{u}_{1}=h_{2}u_{3}+D_{3,2}u_{2}u_{3}+u_{1}u_{3}\beta(t)+\alpha
u_{1}\left(D_{3,1}u_{3}^{2}+D_{2,1}u_{2}^{2}-h_{2}u_{2}\right)\\\
\dot{u}_{2}=-D_{3,1}u_{1}u_{3}+u_{2}u_{3}\beta(t)+\alpha\left(D_{3,2}u_{2}u_{3}^{2}-D_{2,1}u_{2}u_{1}^{2}+h_{2}u_{1}^{2}+h_{2}u_{3}^{2}\right)\\\
\dot{u}_{3}=-h_{2}u_{1}+D_{2,1}u_{1}u_{2}-\left(u_{1}^{2}+u_{2}^{2}\right)\beta(t)-\alpha
u_{3}\left(D_{3,2}y^{2}+D_{3,1}u_{1}^{2}+h_{2}u_{2}\right)\end{cases}\mbox{.}$
(5)

The latter is immediately obtained from (1) by setting $\bm{u}:=\bm{m}$ and
$\bm{h}_{a}=:(0,h_{2},0)$. Given $D_{i}$ satisfying (4), we have defined, for
all $1\leq j<i\leq 3$,

$D_{i,j}:=D_{i}-D_{j}\mbox{.}$

Clearly, these quantities are such that

$D_{3,1}>D_{3,2},\qquad D_{3,2}+D_{2,1}-D_{3,1}=0\mbox{.}$ (6)

The function $\beta(t):[0,+\infty)\rightarrow\mathbb{R}$ (_injected current_)
plays here the role of “control”.
As it is well known (see, e.g. [MBS09]), the function “distance from the
origin”

$\Psi(\bm{u}):=\sum_{j=1}^{3}u_{j}^{2}$ (7)

is conserved along the solutions of (5), namely $\Psi(\bm{u})=1$, meaning that
the motion occurs on the unit-sphere.
A switching process (see [MBS09] for a comprehensive treatment) can be briefly
described as the possibility to “move” the state of the system (5) from (a
neighbourhood of) the equilibrium $\bm{s}^{-}$ to $\bm{s}^{+}$, where
$\bm{s}^{\pm}:=(\pm\sqrt{1-(h_{2}/D_{2,1})^{2}},h_{2}/D_{2,1},0)$ when
$\beta(t)=0$ (straightforward check) and this can be achieved by applying a
suitable injected current $\beta(t)$. We shall say that the “switching” has
taken place, once the system has reached a state which is attracted by
$\bm{s}^{+}$ without the need to supply any further control via $\beta(t)$.
For reasons of topological nature, such a transition cannot exist if the
current is switched off all the time.
As already mentioned in the foreword, several approaches can be found in the
literature. Possibly the most paradigmatic one is the ballistic one, which
consists in applying a current $\beta(t)=\operatorname{const.}>0$ for a
prescribed period of time $t\in[0,T_{e}]$, then switching it off. The time
$T_{e}$ is meant to be be determined.
However, as already stressed, due to the highly involved phase space space
structure of the system at hand pointed out by specialised analyses (see,
e.g., [dPS+15a]) such as a successful switching is more properly intended as a
probabilistic feature of the described technique. In other terms, without any
supporting result of quantitative nature, this class of switching approaches
have an experimental validity only.
In a different spirit, this paper aims to provide a _deterministic_ argument
to successfully complete a switching process under mild assumptions. This has
been possible by recognising some intrinsic perturbative features of the
system at hand, and exploiting them via perturbation methods borrowed from the
Hamiltonian world. As a main achievement, this has led to the realisation of
the mentioned CQL solutions, via a suitably constructed controlling current
$\beta(t)$. As the name suggests, these solutions are characterised by a
quasi-constant (in the perturbative sense) $u_{3}(t)$ throughout the motion,
until a neighbourhood of the target equilibrium $\bm{s}^{+}$ is reached.
Interestingly enough, “small” values of $|u_{3}(t)|$ turn out to provide a
very adequate option, being the solution “rapidly” attracted by the
equilibrium (characterised by $u_{3}=0$), right after the current is switched
off. This constitutes another difference with the ballistic approach, where
the required amount of initial injected current leads to “higher” values of
$|u_{3}(t)|$, requiring in this way a “long time” to come back to zero, as per
effect of the “small” dissipation $\alpha$.

##  4 Perturbative setting and main result

The starting point consists in noticing that a class of realistic models
exhibit different scales amongst the involved parameters. More specifically,
whilst $D_{3,2},D_{3,1}$ are “of order one”, denoted with $O(1)$, the
remaining $h_{2},D_{2,1},\alpha$ are sensibly smaller, namely “of order
$\lambda$”, where $\lambda$ is typically $\sim 10^{-2}$. We choose
$\beta(t)=O(\lambda)$ as it is usually done, since the spin torque compensates
the damping and is of the same order of magnitude. The described feature of
the parameters at hand, naturally leads to a perturbative formulation for (5).
In particular, by defining $h_{2}=:\lambda\tilde{h}_{2}$, with
$\tilde{h}_{2}=O(1)$ and similarly for all the other mentioned $O(\lambda)$
quantities, the system (5) is immediately cast into the following form

$\begin{cases}\dot{u}_{1}=D_{3,2}u_{2}u_{3}+\lambda\left[\tilde{h}_{2}u_{3}+u_{1}u_{3}\tilde{\beta}(t)+\tilde{\alpha}u_{1}D_{3,1}u_{3}^{2}\right]+\tilde{\alpha}\lambda^{2}u_{1}u_{2}\left[\tilde{D}_{2,1}u_{2}-\tilde{h}_{2}\right]\\\
\dot{u}_{2}=-D_{3,1}u_{1}u_{3}+\lambda
u_{2}u_{3}\left[\tilde{\beta}(t)+\tilde{\alpha}D_{3,2}u_{3}\right]+\tilde{\alpha}\lambda^{2}\left[\tilde{h}_{2}\left(u_{1}^{2}+u_{3}^{2}\right)-\tilde{D}_{2,1}u_{1}^{2}u_{2}\right]\\\
\dot{u}_{3}=\lambda\left[-\tilde{h}_{2}u_{1}+u_{1}u_{2}\tilde{D}_{2,1}-\left(u_{1}^{2}+u_{2}^{2}\right)\tilde{\beta}(t)-\tilde{\alpha}u_{3}\left(D_{3,2}u_{2}^{2}+D_{3,1}u_{1}^{2}\right)\right]-\tilde{\alpha}\tilde{h}_{2}\lambda^{2}u_{2}u_{3}\end{cases}$
(8)

or, in a more compact notation,

$\dot{\bm{u}}={\bm{A}}(u_{3})\bm{u}+\lambda\bm{\mathcal{N}}(\bm{u},\lambda;\tilde{\beta}(t))\mbox{,}$
(9)

where

${\bm{A}}(\xi):=\begin{pmatrix}0&D_{3,2}\xi&0\\\ -D_{3,1}\xi&0&0\\\ 0&0&0\\\
\end{pmatrix}$

and the definition of $\bm{\mathcal{N}}$ is obvious.
It is clear that at the order zero in $\lambda$, the variable $u_{3}(t)$ is
constant for all $t$. As a consequence, the dynamics in the variables
$(u_{1},u_{2})$ is described by the equation of a harmonic oscillator with
(constant) frequency determined once and for all by the initial condition
$u_{3}(0)$. We shall restrict ourselves to the choice $u_{3}(0)=-K$, with
$K\in\mathbb{R}^{+}$. We remark that the latter will be achieved in the full
process via the “expulsion” stage, which will be described later on.
Hence, by defining

$\omega:=K\sqrt{D_{3,2}D_{3,1}}$ (10)

the energy of such an oscillator is given by
$E=(\omega^{2}/2)[u_{1}^{2}+(u_{2}/D_{3,1})^{2}]$, and hence a family of
ellipses in the plane $(u_{1},u_{2})$ parameterised by $E$. A motion starting
in $(u_{1}(0),u_{2}(0))=(-1,0)$ (i.e. $\lambda$-close to
$(s_{1}^{-},s_{2}^{-})$) has energy $\hat{E}=\omega^{2}/2$. As a consequence,
after a semi-period

$T_{tr}^{+}:=\pi/\omega\mbox{,}$ (11)

it will evolve in $(u_{1}(T_{tr}^{+}),u_{2}(T_{tr}^{+}))=(1,0)$, which is
$\lambda$-close to $(s_{1}^{+},s_{2}^{+})$. This is, essentially, the key
mechanism of transfer we will rely on.
Note that this represents an archetypal _latitudinal_ solution, as (at least
at the zero-th order in $\lambda$) it takes place along the same latitudinal
line of the unit sphere, determined once and for all by the choice of $K$.
However, should $K$ be “small” (for reasons that will be clarified later on),
the transfer process would take a “long” time, and then the $O(\lambda)$
contributions would become significant. This means that, in order to preserve
the “transfer” feature of the solutions, we need to find a way to “contain”
the variations of $\bm{u}$ from the unperturbed solution.
At this point, we are about to summarise the argument used to construct the
control that is apt to realise the desired switching process.
Let us now consider, for this purpose, the auxiliary system in a renamed set
of variables, say $v_{j}(t)$, which is obtained from (8) when
$\tilde{\beta}(t)$ is replaced by some function of $\bm{v}$, namely
$\hat{\beta}(\bm{v})$ and subject to the same initial conditions. Hence, the
key step to obtain a latitudinal solution consists in choosing $\hat{\beta}$
in such a way the r.h.s. of the third equation is zero. In fact, as a
consequence of this choice, we have $v_{3}(t)=-K$ for all $t$.
Hence, in order to select a class of latitudinal solutions we can choose

$\hat{\beta}_{lat}(\bm{v}):=(v_{1}^{2}+v_{2}^{2})^{-1}\left[v_{1}v_{2}\tilde{D}_{2,1}+\tilde{\alpha}K\left(D_{3,2}v_{2}^{2}+D_{3,1}v_{1}^{2}\right)+\tilde{\alpha}\tilde{h}_{2}\lambda
v_{2}K-\tilde{h}_{2}v_{1}\right]\mbox{,}$ (12)

this leads to

$\dot{\bm{v}}={\bm{A}}(v_{3})\bm{v}+\lambda\bm{\mathcal{N}}(\bm{v},\lambda;\hat{\beta}_{lat}(\bm{v}))=:{\bm{A}}(-K)\bm{v}+\lambda\bm{\mathcal{F}}(\bm{v};\bm{v})+\lambda^{2}\bm{\mathcal{R}}(\bm{v};\bm{v})\mbox{.}$
(13)

In other terms, $\lambda\bm{\mathcal{F}}(\bm{v};\bm{v})$ denotes the first
order terms (in $\lambda$) of the field $\bm{\mathcal{N}}(\bm{v},\cdot;\cdot)$
in which $\hat{\beta}$ has been substituted as a function of $\bm{v}$. It
might be useful, in order to clarify the notation, to denote as
$\lambda\bm{\mathcal{F}}(\bm{u};\bm{v})$ the first order terms arising in (9)
following to the substitution
$\tilde{\beta}\leftarrow\hat{\beta}_{lat}(\bm{v})$. The meaning of
$\bm{\mathcal{R}}$ is obvious as a consequence. Their explicit expression is
slightly more involved and is given in Appendix A, however, as for
$\bm{\mathcal{F}}(\bm{v};\bm{v})$ and $\bm{\mathcal{R}}(\bm{v};\bm{v})$, we
have

$\begin{aligned} \mathcal{F}_{1}(\bm{v};\bm{v})&=-\rho
Kv_{2}\left(\tilde{h}_{2}v_{2}+\tilde{D}_{2,1}v_{1}^{2}\right)\\\
\mathcal{F}_{2}(\bm{v};\bm{v})&=\rho
Kv_{1}\left(\tilde{h}_{2}v_{2}-\tilde{D}_{2,1}v_{2}^{2}\right)\\\
\mathcal{R}_{1}(\bm{v};\bm{v})&=\rho\tilde{\alpha}v_{1}v_{2}\left[K^{2}\tilde{D}_{2,1}v_{2}+v_{1}^{2}(v_{2}\tilde{D}_{2,1}-\tilde{h}_{2})-\tilde{h}_{2}(K^{2}+v_{2}^{2})+\tilde{D}_{2,1}v_{2}^{3}\right]\\\
\mathcal{R}_{2}(\bm{v};\bm{v})&=\rho\tilde{\alpha}v_{1}^{2}\left[-K^{2}\tilde{D}_{2,1}v_{2}-v_{1}^{3}(v_{2}\tilde{D}_{2,1}-\tilde{h}_{2})+\tilde{h}_{2}(K^{2}+v_{2}^{2})-\tilde{D}_{2,1}v_{2}^{3}\right]\end{aligned}\mbox{,}$
(14)

where, by virtue of (7), we have set $v_{1}^{2}+v_{2}^{2}=1-K^{2}=:1/\rho$.
Clearly, by construction

$\mathcal{F}_{3}(\bm{v};\bm{v})=\mathcal{R}_{3}(\bm{v};\bm{v})=0\mbox{.}$ (15)

Hence (13) is, by all means, a two dimensional system of ODEs.
Similarly to what has been observed for (8), (13) is clearly integrable at the
zero-th order in $\lambda$.
Let us imagine for a moment to be able to compute a solution $v_{1,2}(t)$ to
(13). The construction above implies that, by choosing the particular
_control_

$\tilde{\beta}(t):=\hat{\beta}_{lat}(v_{1}(t),v_{2}(t),-K)\mbox{,}$ (16)

and substituting it in (8), one has $\bm{u}(t)=(v_{1}(t),v_{2}(t),-K)$,
provided that (8) and (13) are subject to the same initial conditions. In
other terms, (exact) solutions to (13) yield latitudinal solutions (in the
sense defined above) to (8) via (16). Hence, it is reasonable to expect that
“approximated” solutions to (13) could provide a “quasi-constant” behaviour
for the third variable $u_{3}(t)$, in a sense to be made precise. We shall
refer to these as _controlled quasi-latitudinal_ (CQL) solutions to (8). Under
this setting, the main result states as follows:

######  Theorem 4.1.

There exist $\lambda_{0},r^{-},T_{e},T_{tr}>0$, a non-trivial set of
parameters and a control

$\tilde{\beta}(t):=\begin{cases}\tilde{\beta}_{e}&t\in[-T_{e},0)\\\
\tilde{\beta}_{tr}(t)&t\in[0,T_{tr}]\\\ 0&t\in(T_{tr},+\infty)\end{cases}$
(17)

such that the flow of the system (8) controlled via (17) and denoted with
$\tilde{\Phi}^{t}$, for any $r\leq r^{+}$ and initial condition
$\bm{u}_{0}:=\bm{u}(-T_{e})\in\mathfrak{B}_{r}({\bm{s}^{-}})$, satisfies

$\lim_{t\rightarrow+\infty}\tilde{\Phi}^{t}(\bm{u}_{0})\in\mathfrak{B}_{fr}({\bm{s}^{+}})\mbox{,}$
(18)

for some $f>0$. In particular, the system realises a CQL solution in the time
interval $[0,T_{tr}]$.

The rest of the paper is devoted to the proof of Thm. 4.1. This will be
achieved in three steps. The first one, called “expulsion”, consists in
choosing a suitable neighbourhood of the point $\bm{s}^{-}$ and injecting a
constant current of magnitude $\tilde{\beta}_{e}$ for a certain time $T_{e}$.
It is shown how this “pushes” the third variable $u_{3}(t)$ in a region in
which it is strictly negative, a key requirement in order to obtain the class
of switching motions we are interested in, even at the zero-th order, as
anticipated in sec. 4.
A second stage, referred to as “transfer”, represents the very heart of the
argument, i.e. where the QLS are constructed via a suitable controlling
current, with the property of being able to “deliver” the solutions starting
in the vicinity of the point reached once the first step is complete, in a
neighbourhood of the target point $\bm{s}^{+}$.
The last step, called “attraction”, consists in showing that $\bm{s}^{+}$ is
an attractive point when the current is switched off, attracting in this way
any points arriving in its vicinity. Although this is a very well known
property, quantitative information are needed for a precise formulation of the
result, and this is exactly the aim of this third stage. It is immediate to
realise that being able to carry out the described three steps, implies a
successful realisation of the switching process. The whole argument clearly
exploits the known group property of a dynamical system.
From a technical viewpoint, the “core” of the proof (consisting in the
“transfer” stage), relies on the possibility to construct a suitable control
for the system at hand, by approximating the non-linear flow via Hamiltonian
perturbative tools. This is possible via the well known possibility to
interpret any system of ODEs as a Hamiltonian system in a suitably extended
phase space, see e.g. [Ber09]. However, such a simple observation has been
profitably used in several cases, due to the potential of a full all-orders
generalisation of the perturbative setting. See, for instance, [FW16], [FW19].
This has the potential to increase the threshold $\lambda_{0}$, expanding in
this way the class of systems which can be dealt with via this approach.
The remaining stages, “expulsion” and “attraction”, are carried out by using
established tools from the theory of ODEs that will be specified later.

##  5 Expulsion

Let us consider the system (9) and set
$\tilde{h}_{2}=:-\tilde{D}_{2,1}\Omega$, with $\Omega>0$. In this setting, in
particular, one has $\bm{s}^{\pm}\equiv(\pm\gamma,-\Omega,0)$, where
$\gamma:=\sqrt{1-\Omega^{2}}$. Let us now perform a translation-rescaling of
$\hat{\bm{s}}^{-}$ at the origin of the new system of coordinates

$\lambda\bm{\xi}:=\bm{u}-\bm{s}^{-}\mbox{.}$ (19)

Throughout this section we shall consider the autonomous control
$\lambda\tilde{\beta}(t)=:\lambda\tilde{\beta}_{e}(t)\equiv{\beta}_{e}=\operatorname{const.}$
and the new time $\tau:=T_{e}+t$, for any $T_{e}>0$. Consequently,
$\Phi_{e}^{\tau}$ will denote the phase flow of (9) with respect to the new
time. Hence, by (19), (9) reads as

$\frac{d}{d\tau}\bm{\xi}={\bm{\mathcal{L}}}\bm{\xi}+\bm{f}+\lambda\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})\mbox{,}$
(20)

where,

${\bm{\mathcal{L}}}:=\begin{pmatrix}0&0&\bar{a}\\\ 0&0&\bar{b}\\\
2\beta_{e}\gamma&2\beta_{e}\Omega&0\\\
\end{pmatrix}\mbox{,}\qquad\bm{f}=(0,0,-\beta_{e}/\lambda)^{\top}\mbox{,}$
$\bar{a}:=-D_{3,2}\Omega-\beta_{e}\gamma,\qquad\bar{b}:=D_{3,2}\gamma-\beta_{e}\Omega\mbox{,}$
(21)

and $\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})$ is defined as a
consequence. Let us remark that we have used the relation
$D_{3,1}=D_{3,2}-\lambda\tilde{D}_{2,1}$, see (6), in order to simplify the
structure of ${\bm{\mathcal{L}}}$. As a result, the linear term
$-\tilde{D}_{2,1}\xi_{3}$ is included in
$\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})$. In particular, it is evident
that for sufficiently small $\lambda$ the dominant behaviour is “expulsive”
i.e. it drives $\xi_{3}$ towards negative values.
As it can be easily checked, by defining

${\bm{S}}:=\begin{pmatrix}1&1&1\\\
\bar{b}/\bar{a}&\bar{b}/\bar{a}&-\gamma/\Omega\\\
\sqrt{2}i\beta_{e}/\bar{a}&-\sqrt{2}i\beta_{e}/\bar{a}&0\\\ \end{pmatrix}$

one has
${\bm{S}}^{-1}{\bm{\mathcal{L}}}{\bm{S}}={\bm{\Gamma}}\equiv\operatorname{diag}(-\sqrt{2}i\beta_{e},\sqrt{2}i\beta_{e},0)$.
Hence, the solution to (9) reads as

$\bm{\xi}(\tau)=e^{{\bm{\mathcal{L}}}\tau}\bm{\xi}(0)+\int_{0}^{T_{e}}e^{{\bm{\mathcal{L}}}(\tau-s)}\bm{f}ds+\lambda\int_{0}^{T_{e}}e^{{\bm{\mathcal{L}}}(\tau-s)}\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})ds\mbox{,}$
(22)

where, by noticing that $\bar{a}\gamma+\bar{b}\Omega=-\beta_{e}$, one has

$e^{{\bm{\mathcal{L}}}\tau}=\frac{1}{\beta_{e}}\begin{pmatrix}-\bar{b}\Omega-\bar{a}\gamma\cos(\sqrt{2}\tau\beta_{e})&\bar{a}\Omega-\bar{a}\Omega\cos(\sqrt{2}\tau\beta_{e})&(\bar{a}/\sqrt{2})\sin(\sqrt{2}\tau\beta_{e})\\\
\bar{b}\gamma-\bar{b}\gamma\cos(\sqrt{2}\tau\beta_{e})&-\bar{b}\Omega\cos(\sqrt{2}\tau\beta_{e})-\bar{a}\gamma&(\bar{b}/\sqrt{2})\sin(\sqrt{2}\tau\beta_{e})\\\
\sqrt{2}\beta_{e}\gamma\sin(\sqrt{2}\tau\beta_{e})&\sqrt{2}\beta_{e}\Omega\sin(\sqrt{2}\tau\beta_{e})&\beta_{e}\sin(\sqrt{2}\tau\beta_{e})\\\
\end{pmatrix}\mbox{.}$

Furthermore, the norm of the latter is uniformly bounded in time, i.e.

$||e^{{\bm{\mathcal{L}}}\tau}||\leq\mathcal{M}_{e}\mbox{,}$ (23)

for some suitable
$\mathcal{M}_{e}=\mathcal{M}_{e}(\beta_{e},D_{3,2},\gamma,\Omega)\geq 1$ and
all $\tau\in\mathbb{R}$.
It is now easy to realise that, if $O(\lambda)$ are disregarded in (22), the
integral appearing in it can be computed immediately and this gives rise to a
function, namely $\bm{\xi}^{\prime}(\tau)$, which approximates the “full”
$\bm{\xi}(\tau)$ for “small” $\lambda$ (under suitable assumptions). Such a
function reads as

$\bm{\xi}^{\prime}(\tau)=e^{{\bm{\mathcal{L}}}\tau}\bm{\xi}^{\prime}(0)+(2\beta_{e}\lambda)^{-1}\left(\bar{a}(\cos(\sqrt{2}t\beta_{e})-1),\bar{b}(\cos(\sqrt{2}t\beta_{e})-1),-\sqrt{2}\beta_{e}\sin(\sqrt{2}t\beta_{e})\right)^{\top}\mbox{.}$
(24)

As the state of the system is supposed to be initially in a neighbourhood of
the origin, it is natural to set $\bm{\xi}^{\prime}(0)=\bm{0}$. Hence, in the
original set of variables, the third variable evolves as
$u_{3}(\tau)\sim-2^{-1/2}\sin(\sqrt{2}t\beta_{e})$ for initial data “close to”
$\bm{s}^{-}$ and “sufficiently small” $\lambda$. Once a target value $-K$ for
$u_{3}$ has been set, with

$0<\sqrt{2}K\leq 1\mbox{,}$ (25)

the expulsion stage is defined as the evolution of the system for
$\tau\in[0,T_{e}]$, where

$T_{e}:=(\sqrt{2}\beta_{e})^{-1}\arcsin(\sqrt{2}K)\mbox{.}$ (26)

From (24) one immediately gets

$\bm{\xi}^{\prime}(T_{e})=(2\beta_{e}\lambda)^{-1}(\bar{a}(-1+\sqrt{1-2K^{2}}),\bar{b}(-1+\sqrt{1-2K^{2}}),-2K\beta_{e})\mbox{.}$
(27)

######  Remark 5.1.

For realistic values of the parameters (for instance,
$\beta_{e},K^{2}=O(\lambda)$), the quantities
$\lambda\xi_{1,2}^{\prime}(T_{e})$ are $O(1)$, i.e. the point reached after
the expulsion stage is “far” from $\bm{s}^{-}$. For this reason, approximation
formula (27) will play a key role in the main proof. Furthermore, the latter
has some interest itself in ballistic switching processes, in which “long”
expulsions are typically considered.

The approximation features of (27), numerically validated for an example in
fig. 2, are stated in a quantitative form in the next

######  Lemma 5.1.

Let us set a target value $K$ satisfying (25) and $\beta_{e}>0$ arbitrarily
chosen. Then choose $T_{e}$ as in (26). Let us now define

$r^{*}:=(2\beta_{e}\lambda)^{-1}\left[4K^{2}+(D_{3,2}^{2}+\beta_{e}^{2})(\sqrt{1-2K^{2}}-1)^{2}\right]^{\frac{1}{2}}$
(28)

and let be $\rho_{e}>0$ arbitrarily chosen. Now set

$\lambda\mathcal{M}_{1}:=\max_{\bm{\xi}\in\mathfrak{B}_{r^{*}}(\bm{0})}|\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})|\mbox{,}\qquad\lambda\mathcal{M}_{2}:=\max_{\bm{\xi}\in\mathfrak{B}_{r^{*}+\rho_{e}}(\bm{0})}||{\bm{D}}\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})||_{\infty}\mbox{,}$
(29)

see Appendix B. for the explicit expressions of
${\bm{D}}\bm{\mathcal{V}}(\bm{\xi};\lambda,\beta_{e})$.
Then for all $\lambda\in(0,\lambda_{e}]$, with

$\lambda_{e}:=(T_{e}\mathcal{M}_{e})^{-1}\min\\{(4\mathcal{M}_{1})^{-1}\rho_{e},\mathcal{M}_{2}^{-1}\log
2\\}\mbox{,}$ (30)

the following property holds

$\Phi^{T_{e}}(\mathfrak{B}_{(4\mathcal{M}_{e})^{-1}\rho_{e}\lambda}(\bm{s}^{-}))\subseteq\mathfrak{B}_{\rho_{e}\lambda}(\bm{u}_{c}(T_{e}))\mbox{,}$
(31)

where $\bm{u}_{c}(\tau):=\lambda\bm{\xi}^{\prime}(\tau)+\bm{s}^{-}$. In other
terms, solutions starting within a sphere of radius $r^{-}$ around
$\bm{s}^{-}$ remain contained in a sphere of radius $\rho_{e}\lambda$ of the
(known) point $\lambda\bm{\xi}^{\prime}(T_{e})+\bm{s}^{-}$.

######  Remark 5.2.

As it can be easily noticed from (28), $r^{*}|_{K=0}=0$ and it grows
monotonically with $K$. The values of $\mathcal{M}_{1,2}$ will increase
accordingly. This implies that the allowed threshold for $\lambda_{e}$ worsens
as the target $K$ gets bigger.

###### Proof.

The proof relies on the fact that the solutions $\bm{\xi}(\tau)$, starting
from a suitable neighbourhood of the origin, do not escape the sphere of
radius $r^{*}+\rho_{e}$ centred at the origin, in such a way the bounds (29)
are justified. For this purpose, let us firstly check that
$|\bm{\xi}^{\prime}(T_{e})|\leq r^{*}$. This is immediate from (27),
assumption (28) (notice that $\bar{a}^{2}+\bar{b}^{2}=D_{3,2}^{2}+\Omega^{2}$)
and finally observing that $|\bm{\xi}^{\prime}(\tau)|$ is monotonically
increasing for all $\tau\in[0,T_{e}]$.
Let us now define $\bm{\delta}(\tau):=\bm{\xi}(\tau)-\bm{\xi}^{\prime}(\tau)$.
By substituting in (22) we get

$\bm{\delta}(\tau)=e^{{\bm{\mathcal{L}}}\tau}\bm{\delta}(0)+\lambda\int_{0}^{T_{e}}e^{{\bm{\mathcal{L}}}(\tau-s)}\bm{\mathcal{V}}(\bm{\xi}^{\prime}+\bm{\delta};\lambda,\beta_{e})ds\mbox{,}$
(32)

where $\bm{\delta}(0)\equiv\bm{\xi}(0)$. By taking the absolute values of
(32), then using the bound
$|\bm{\mathcal{V}}(\bm{\xi}^{\prime}+\bm{\delta};\lambda,\beta_{e})|\leq|\bm{\mathcal{V}}(\bm{\xi}^{\prime};\lambda,\beta_{e})|+\mathcal{M}_{2}|\bm{\delta}|$,
assumptions (29), and finally the classical Grönwall lemma, one obtains

$|\bm{\delta}(\tau)|\leq\mathcal{M}_{e}(|\bm{\delta}(0)|+\lambda\mathcal{M}_{1}\tau)\exp(\lambda\mathcal{M}_{e}\mathcal{M}_{2}\tau)\mbox{.}$

Hence, for all $\lambda\in(0,\lambda_{e}]$ and all $\tau\in(0,T_{e}]$ it is
sufficient to choose $|\bm{\delta}(0)|\leq(4\mathcal{M}_{e})^{-1}\rho_{e}$ as
suggested by the l.h.s. of (31), in order to get
$|\bm{\delta}(\tau)|\leq\rho_{e}$. This proves the r.h.s. of (31). ∎

\begin{overpic}[width=455.24408pt]{ExpulsionStereo.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(5.0,68.0){\footnotesize$w_{1}$}
\put(60.0,24.0){\footnotesize$w_{2}$} \end{overpic}

\begin{overpic}[width=455.24408pt]{ExpulsionErrors.pdf}
\put(0.0,78.0){\footnotesize(b)} \put(5.0,68.0){\footnotesize$|\bm{\delta}|$}
\put(62.0,24.0){\footnotesize$t$}
\put(70.0,71.0){\footnotesize$\lambda_{(1)}$}
\put(70.0,63.0){\footnotesize$\lambda_{(2)}$}
\put(70.0,55.0){\footnotesize$\lambda_{(3)}$} \end{overpic}

Figure 2: In panel (a) the trajectories $\bm{\xi}(\tau)$ (continuous line) and
$\bm{\xi}^{\prime}(\tau)$ (dashed line) are reported in the stereographic
coordinates defined via $w_{1,2}:=u_{1,2}/(1+u_{3})$ and transformation (19).
Panel (b) shows the behaviour of the error $|\bm{\delta}(t)|$ for
$t\in[-T_{e},0]$ with $T_{e}=2.3372$ and three different values of $\lambda$:
$\lambda_{(1)}=0.006$, $\lambda_{(2)}=0.004$ and $\lambda_{(3)}=0.001$. The
corresponding values for $\tilde{\beta}_{e}$ are $0.0018$, $0.0012$ and
$0.003$, respectively. The remaining parameters used are specified in Appendix
C.

##  6 Transfer

Let us firstly introduce the following notation. Given any vector
$\bm{x}\in\mathbb{R}^{3}$ we shall denote with $\bm{x}_{r}:=(x_{1},x_{2})$.
Vice-versa, we shall denote with $w_{1,2}$ the first two components of either
$\bm{w}_{r}$ or $\bm{w}$. It is immediate to notice from (14) and (15) that
$\bm{\mathcal{F}}(\bm{v};\bm{v})\equiv(\bm{\mathcal{F}}_{r}(\bm{v}_{r};\bm{v}_{r}),0)$
and similarly for $\bm{\mathcal{R}}$. Hence, it is meaningful to consider the
following system

$\dot{\bm{w}}_{r}={\bm{L}}\bm{w}_{r}+\lambda\bm{\mathcal{F}}_{r}(\bm{w}_{r},\bm{w}_{r})$
(33)

where

$\bm{L}:=\begin{pmatrix}0&-\sigma\omega\\\ \omega/\sigma&0\\\
\end{pmatrix},\qquad\sigma:=\sqrt{D_{3,2}/D_{3,1}}\mbox{,}$ (34)

and $\omega$ has been defined in (10). Clearly, (33) is nothing but the first
order truncation of (13), written in the reduced set of variables
$\bm{w}_{r}$. Moreover, it is immediate to check that

$\mathcal{G}_{tr}:=\sigma^{-1}w_{1}^{2}+\sigma w_{2}^{2}\mbox{,}$ (35)

is a prime integral for (33) if $\lambda=0$.
With the aim to construct CQL solutions, we ask ourselves if the system can be
solved, for instance, by means of perturbative tools. The answer is
affirmative, as stated in the following

######  Proposition 6.1.

It is possible to construct a function $\bm{w}_{r}^{[\leq 1]}(t)$ satisfying
(33) up to $O(\lambda)$.

######  Remark 6.1.

It is important to avoid any ambiguity about the meaning of the previous
statement. Solving (33) up to $O(\lambda)$ does not mean solving it exactly.
In fact one has $\dot{\bm{w}}_{r}^{[\leq 1]}-{\bm{L}}\bm{w}_{r}^{[\leq
1]}-\lambda\bm{\mathcal{F}}_{r}(\bm{w}_{r}^{[\leq 1]},\bm{w}_{r}^{[\leq
1]})=O(\lambda^{2})$. This is typical of perturbative arguments, which are
known to generate a remainder as a consequence of the expansions involved.

###### Proof.

Let us firstly cast the linear part of (33) into a diagonal form. It is
immediate to check that the required transformation is given by

$\bm{w}_{r}=\bm{C}\bm{x},\qquad\bm{C}:=\begin{pmatrix}\sigma&\sigma\\\ -i&i\\\
\end{pmatrix}\mbox{.}$ (36)

In fact,
$\bm{C}^{-1}\bm{L}\bm{C}=\operatorname{diag}(i\omega,-i\omega)=:\bm{\Lambda}$.
We can now cast system (33) in the new set of variables

$\dot{\bm{x}}=\bm{\Lambda}\bm{x}+\lambda\bm{C}^{-1}\bm{\mathcal{F}}_{r}(\bm{C}\bm{x})\mbox{,}$
(37)

into a Hamiltonian form via a phase space extension. More precisely, by
denoting with $y_{j}$ the momenta canonically conjugated to $x_{j}$, one has
that (37) is given by (part of) the canonical equations of

$H(\bm{y},\bm{x})=H_{0}(\bm{y},\bm{x})+\lambda H_{1}(\bm{y},\bm{x})\mbox{,}$

where

$\begin{aligned} H_{0}&:=i\omega(x_{1}y_{1}-x_{2}y_{2})\\\ H_{1}&:=K\rho
y_{1}\left[\frac{\tilde{h}_{2}}{2}\left(\sigma+\frac{1}{\sigma}\right)x_{1}^{2}+\frac{\tilde{h}_{2}}{2}\left(\frac{1}{\sigma}-\sigma\right)x_{2}^{2}-\frac{\tilde{h}_{2}}{\sigma}x_{1}x_{2}-i\tilde{D}_{2,1}\sigma
x_{1}x_{2}^{2}+i\tilde{D}_{2,1}\sigma x_{1}^{3}\right]\\\ &+K\rho
y_{2}\left[\frac{\tilde{h}_{2}}{2}\left(\sigma+\frac{1}{\sigma}\right)x_{2}^{2}+\frac{\tilde{h}_{2}}{2}\left(\frac{1}{\sigma}-\sigma\right)x_{1}^{2}-\frac{\tilde{h}_{2}}{\sigma}x_{1}x_{2}+i\tilde{D}_{2,1}\sigma
x_{1}^{2}x_{2}-i\tilde{D}_{2,1}\sigma x_{2}^{3}\right]\end{aligned}\mbox{.}$
(38)

As it is common in perturbation theory, we ask weather it is possible to find
a canonical, $\lambda-$close to the identity transformation of variables
$(\bm{x},\bm{y})=\mathcal{T}(\bm{X},\bm{Y})$ apt to “remove” the contribution
of $H_{1}$ i.e. such that

$H\circ\mathcal{T}=H_{0}(\bm{Y},\bm{X})+O(\lambda^{2})\mbox{.}$ (39)

For this purpose, by invoking the well known Gröbner Exchange Theorem, see
e.g. [Gio03], such a transformation will be determined by requiring that

$\exp(\mathcal{L}_{\lambda\chi})H=H_{0}+O(\lambda^{2})\mbox{,}$

i.e. the well known first order homological equation

$\\{\chi,H_{0}\\}=H_{1}\mbox{,}$ (40)

where the generating function $\chi$ will be sought of the form

$\chi(\bm{y},\bm{x})=\bm{y}\cdot\bm{\mathcal{C}}(\bm{x})\mbox{.}$ (41)

This choice is suggested by a general property when dealing with Hamiltonians
obtained via the above described phase space extension, see [FW16] for a
proof. In this particular case, we will look for $\mathcal{C}_{j}$ as non-
homogeneous polynomials of the form

$\mathcal{C}_{j}(\bm{x}):=\sum_{|\bm{\nu}|=2,3}c_{\bm{\nu}}^{(j)}\bm{x}^{\bm{\nu}}\mbox{,}$
(42)

where $\bm{\nu}\in\mathbb{N}^{2}$, $|\nu|:=\nu_{1}+\nu_{2}$,
$\bm{x}^{\bm{\nu}}:=x_{1}^{\nu_{1}}x_{2}^{\nu_{2}}$ and $c_{\bm{\nu}}^{(j)}$
are complex-valued unknown coefficients to be determined. By using (42) in
(41), then substituting in (40) one finds

$\displaystyle-c_{0,2}^{(1)}=c_{2,0}^{(2)}$
$\displaystyle=i\tilde{h}_{2}K\rho(\sigma^{2}-1)/(6\sigma\omega)$ (43)
$\displaystyle-c_{2,0}^{(1)}=c_{0,2}^{(2)}$
$\displaystyle=i\tilde{h}_{2}K\rho(\sigma^{2}+1)/(2\sigma\omega)$
$\displaystyle-c_{1,1}^{(1)}=c_{1,1}^{(2)}$
$\displaystyle=i\tilde{h}_{2}K\rho\sigma/(\sigma\omega)$

and

$\displaystyle c_{3,0}^{(1)}=c_{0,3}^{(2)}=c_{1,2}^{(1)}=c_{2,1}^{(1)}$
$\displaystyle=\tilde{D}_{2,1}K\rho\sigma/(2\omega)$ (44) $\displaystyle
c_{3,0}^{(1)}=c_{0,3}^{(2)}$ $\displaystyle=0$

The clear symmetry relations amongst these coefficients is related to the well
known property $x_{1}=\bar{x}_{2}$ (here $\bar{z}$ denotes the complex-
conjugated of $z\in\mathbb{C}$), which is typical of the coordinates maps
(36). In conclusion, if $O(\lambda^{2})$ terms are disregarded, the
normalising transformation reads as

$\bm{x}=\bm{\mathcal{N}}_{[\leq
1]}(\bm{X}):=\bm{X}+\lambda\bm{\mathcal{C}}(\bm{X})\mbox{.}$ (45)

If the above mentioned Gröbner Theorem is not used, checking that (45)
satisfies (39) is just a matter of patience.

######  Remark 6.2.

Interestingly enough, the variables $\bm{y}$ and $\bm{Y}$ do not appear in
(45), which consistently becomes a normalising transformation in the pair
$(\bm{x},\bm{X})$ alone for the original system (37).

It is a simple consequence of the construction above that if $\bm{X}(t)$ is a
solution to the normalised system (up to $O(\lambda)$), then
$\bm{x}(t):=\bm{\mathcal{N}}_{[\leq 1]}(\bm{X}(t))$ is a solution to the
original system (up to $O(\lambda)$).
On the other hand, the integration of the normalised system is immediate

$X_{1,2}(t)=(X_{a}\pm iX_{b})e^{\pm it\omega}\mbox{,}$

where

$X_{a}:=\Re X_{1}(0),\qquad X_{b}:=\Im X_{1}(0)\mbox{,}$ (46)

(recall that $X_{1}=\bar{X}_{2}$). Hence, by using (36), one gets

$\begin{aligned}
w_{1}(t)&=2\sigma(X_{a}\cos(t\omega)-X_{b}\sin(t\omega))+\lambda\omega^{-1}K\rho\left\\{\tilde{D}_{2,1}\sigma^{2}(X_{a}^{2}+X_{b}^{2})[X_{a}\cos(t\omega)-X_{b}\sin(t\omega)]\right.\\\
&+(2/3)\tilde{h}_{2}(\sigma^{2}+2)[2X_{a}X_{b}\cos(2t\omega)-(X_{b}^{2}-X_{a}^{2})\sin(2t\omega)]\\\
&\left.-\tilde{D}_{2,1}\sigma^{2}[X_{a}(3X_{b}^{2}-X_{a}^{2})\cos(3t\omega)-X_{b}(X_{b}^{2}-3X_{a}^{2})\sin(3t\omega)]\right\\}\\\
w_{2}(t)&=2(X_{b}\cos(t\omega)+X_{a}\sin(t\omega))+\lambda\omega^{-1}K\rho\left\\{-\tilde{D}_{2,1}\sigma(X_{a}^{2}+X_{b}^{2})[X_{b}\cos(t\omega)+X_{a}\sin(t\omega)]\right.\\\
&+(2/3)\tilde{h}_{2}(2\sigma^{2}+1)[(X_{b}^{2}-X_{a}^{2})\cos(2t\omega)+2X_{a}X_{b}\sin(2t\omega)]\\\
&\left.-\tilde{D}_{2,1}\sigma[X_{b}(X_{b}^{2}-3X_{a}^{2})\cos(3t\omega)+X_{a}(3X_{b}^{2}-X_{a}^{2})\sin(3t\omega)]-2\tilde{h}_{2}(X_{a}^{2}+X_{b}^{2})/\sigma\right\\}\\\
\end{aligned}\mbox{.}$ (47)

Clearly, the functions above provide the required approximated solution to a
given Cauchy problem for (33), once the quantities $X_{a,b}$ have been
determined from the initial condition $\bm{w}_{r}(0)$ (which will be specified
later on). More precisely, by (36) and (45) one has

$\bm{X}(0)=\bm{\mathcal{N}}_{[\leq 1]}^{-1}(\bm{C}^{-1}\bm{w}_{r}(0))\mbox{,}$
(48)

where we recall that, if $O(\lambda^{2})$ terms are disregarded, one simply
has $\bm{\mathcal{N}}_{[\leq
1]}^{-1}(\bm{x})=\bm{x}-\lambda\bm{\mathcal{C}}(\bm{x})$. The values of
$X_{a,b}$ are finally determined by recalling (46). These quantities determine
the required $\bm{w}_{r}^{[\leq 1]}(t)$ via (47). ∎

######  Remark 6.3.

Clearly, the normalisation order could be increased by generalising the above
described procedure, for instance by using the Lie Series or Lie Transform
methods, see e.g. [Gio03]. The Hamiltonian formulation plays a key role from
this point of view.
Knowingly, the computations get dramatically more and more involved as the
normalisation order is increased.

Let us now notice that at the zero-th order in $\lambda$, (48) yields

$\bm{X}^{[0]}(0)=2^{-1}(\sigma^{-1}w_{1}(0),w_{2}(0))\mbox{,}$

hence, from (47), one gets

$\bm{w}_{r}^{[0]}(t)=(w_{1}(0)\cos(\omega t)-w_{2}(0)\sigma\sin(\omega
t),w_{1}(0)\sigma^{-1}\sin(\omega t)+w_{2}(0)\cos(\omega t))\mbox{.}$

As $w_{1}(0)$ will be assumed as negative and bounded away from zero, the
first component can be written as $\mathcal{A}_{m}\cos(\omega t+\phi)$, where
$\phi:=\arctan(\sigma w_{2}(0)/w_{1}(0))$ and

$\mathcal{A}_{m}:=w_{1}(0)\sqrt{1+(\sigma w_{2}(0)/w_{1}(0))^{2}}<0\mbox{.}$
(49)

Hence, we shall define $T_{tr}$ as the time necessary for $w_{r,1}^{[0]}(t)$
to reach $-\mathcal{A}_{m}-K^{2}$. This yields

$T_{tr}:=\omega^{-1}\left[\arccos(-1-K^{2}/\mathcal{A}_{m})-\arctan(\sigma
w_{2}(0)/w_{1}(0))\right]\mbox{.}$ (50)

It is immediate to realise that this quantity is correctly defined as we we
shall only consider cases in which $K^{2}\ll|w_{1}(0)|\leq|\mathcal{A}_{m}|$.
Moreover, recalling (11), it is easy to check that

$T_{tr}\leq T_{tr}^{+}\mbox{.}$ (51)

In fact, $T_{tr}^{+}$ is attained in the (unrealistic as $D_{2,1}>0$) case
$\bm{w}_{r}(0)=(-1,0)$ i.e $\gamma=1$. However, as it is evident from (19) and
(22), the expulsion stage acts by increasing $w_{2}(t)$ (and hence decreasing
$w_{1}(t)$, and this clearly reduces $T_{tr}$.
Once the control has been constructed, it is important to obtain an estimate
of the measure of the initial conditions around
$\bm{w}(0)\equiv(\bm{w}_{r}(0),-K)$ which are “safely transported” in a
neighbourhood of the target point $\bm{w}(T_{tr})$. In order to achieve this,
we need to bound the difference between the solutions of two different
systems. The first one is given by (33) (i.e. the $O(\lambda)$ truncation of
(13)), whose solution is known by Prop. 6.1. By setting
$\bm{w}:=(\bm{w}_{r}^{[\leq 1]},-K)$ and disregarding $O(\lambda^{2})$, this
reads as

$\dot{\bm{w}}=:{\bm{A}}(-K)\bm{w}+\lambda\bm{\mathcal{F}}(\bm{w};\bm{w})\mbox{,}$
(52)

being subject to the prefixed initial condition $\bm{w}(0)$. The second system
is the full, controlled one, i.e. obtained from (8) by setting
$\tilde{\beta}(t):=\hat{\beta}_{lat}(\bm{w}^{[\leq 1]}(t))$. Its $O(\lambda)$
truncation can be written as

$\dot{\bm{u}}={\bm{A}}(u_{3})\bm{u}+\lambda\bm{\mathcal{F}}(\bm{u};\bm{w})\mbox{,}$
(53)

with an “uncertain” initial condition around $\bm{w}(0)$. Let us denote by
$\Phi_{tr}^{t}$ the flow of the controlled system (53). An example of the
constructed CQL solutions and the corresponding control is shown in fig. 3.

######  Lemma 6.1 ($O(\lambda)$ cut-off).

Let us define the following quantities

$\mathcal{M}_{tr}:=1+\max_{|\bm{u}|\leq
2}||{\bm{D_{\bm{u}}\mathcal{F}}}(\bm{u};\bm{w})||_{\infty},\qquad\mathcal{K}_{\bm{w}}:=1+\max_{t\in[0,T_{tr}]}|\bm{w}_{r}^{[\leq
1]}|\mbox{.}$

Then, for all $\lambda\in(0,\lambda_{tr}]$, with

$\lambda_{tr}:=(4\pi^{2}\mathcal{K}_{\bm{w}}\mathcal{M}_{tr})^{-1}\omega^{2}\mbox{,}$
(54)

and all $\rho_{tr}^{+}>0$, there exists
$\rho_{tr}^{-}=\rho_{tr}^{-}(\rho_{tr}^{+})$, such that following property
holds

$\Phi_{tr}^{T_{tr}}(\mathfrak{B}_{\lambda\rho_{tr}^{-}}(\bm{w}(0))\subset\mathfrak{B}_{\lambda\rho_{tr}^{+}}(\bm{w}(T_{tr}))\mbox{.}$
(55)

######  Remark 6.4.

The $O(\lambda^{2})$ terms in (53) can be bounded by using the Lie series
theory, and this requires some extra work. However, it is possible to check
that the information carried by them can be thought as “negligible”. More
precisely, by facing slightly more cumbersome estimates, one could repeat the
proof below by including these $O(\lambda^{2})$ contributions, and verify that
those terms would imply an additional $O(\lambda)$ term in (54). Hence, the
above statement still holds for the “non cut-off” problem at a price of a
possible (but “small”) further restriction of the threshold $\lambda_{tr}$.

###### Proof.

This proof has some similarities with the one of Lem. 5.1. Let us now define
$\bm{\eta}:=\bm{u}-\bm{w}$, hence, by (52) and (53) one gets

$\dot{\bm{\eta}}={\bm{A}}(-K)\bm{\eta}+{\bm{A}}(\eta_{3})(\bm{w}+\bm{\eta})+\lambda[\bm{\mathcal{F}}(\bm{w}+\bm{\eta};\bm{w})-\bm{\mathcal{F}}(\bm{w};\bm{w})]\mbox{.}$
(56)

It is evident from the latter and the structure of $\bm{A}$ that
$\dot{\eta}_{3}$ has an $O(\lambda)$ magnitude, hence it is appropriate to
consider the new variable $\eta_{3}=:\lambda\zeta_{3}$. By substituting it in
(56), one gets that the $O(1)$ evolution of $\eta_{1,2}$ is given by the
linear system $\dot{\bm{\eta}}_{r}={\bm{L}}\bm{\eta}_{r}$, where
$\bm{\eta}_{r}:=(\eta_{1},\eta_{2})$.
Let us now introduce $d_{r}^{\pm},d_{3}^{\pm}\in(0,1]$, which will play the
role of real parameters to be determined, then suppose for a moment that the
following bound holds

$|\zeta_{3}|\leq d_{3}^{+}\mbox{.}$ (57)

By multiplying both sides of the first two components of equation (56) by the
“integrating factor” $\exp(-{\bm{L}}t)$, integrating and finally considering
the absolute values, one obtains

$\displaystyle|\bm{\eta}_{r}(t)|$
$\displaystyle\leq|\bm{\eta}_{r}(0)|+\lambda\int_{0}^{t}\left\\{|\bm{A}(\zeta_{3})(\bm{w}+\bm{\eta})|+\mathcal{M}_{tr}|\bm{\eta}|\right\\}ds$
$\displaystyle\leq|\bm{\eta}_{r}(0)|+\lambda
d_{3}^{+}\mathcal{K}_{\bm{w}}t+\lambda(d_{3}+\mathcal{M}_{tr})\int_{0}^{t}|\bm{\eta}_{r}(s)|ds$

where we have used that $|\exp({\bm{L}}\cdot)\bm{z}|=|\bm{z}|$, for all
$\bm{z}\in\mathbb{R}^{3}$ as ${\bm{L}}$ has pure complex eigenvalues only, and
the bound

$|\mathcal{F}_{i}(\bm{w}+\bm{\eta};\bm{w})-\mathcal{F}_{i}(\bm{w};\bm{w})|\leq|\bm{\mathcal{F}}(\bm{w}+\bm{\eta};\bm{w})-\bm{\mathcal{F}}(\bm{w};\bm{w})|\leq\mathcal{M}_{tr}|\bm{\eta}|\mbox{.}$
(58)

Hence, by using the Grönwall lemma,

$|\bm{\eta}_{r}(t)|\leq\left[|\bm{\eta}_{r}(0)|+\lambda
d_{3}^{+}\mathcal{K}_{\bm{w}}t\right]e^{\lambda(d_{3}^{+}+\mathcal{M}_{tr})t}\mbox{.}$
(59)

Recalling that $T_{tr}\leq\pi/\omega$ by (51), by using assumption (54) one
gets $\lambda(1+\mathcal{M}_{tr})t\leq 2\lambda\mathcal{M}_{tr}t\leq\log(2)$
for all $t\leq T_{tr}$ as $\mathcal{M}_{tr}\geq 1$. Hence, in order to obtain

$|\bm{\eta}_{r}(t)|\leq\lambda d_{r}^{+}\mbox{,}$ (60)

it is sufficient to require $|\bm{\eta}_{r}(0)|\leq\lambda d_{r}^{-}$, and
then

$d_{r}^{-}\leq d_{r}^{+}/4,\qquad d_{3}^{+}\leq\omega
d_{r}^{+}/(\pi\mathcal{K}_{\bm{w}})\mbox{.}$ (61)

Let us now consider the third equation of (56), i.e.
$\dot{\zeta}_{3}=\mathcal{F}_{3}(\bm{w}+\bm{\eta};\bm{w})-\mathcal{F}_{3}(\bm{w};\bm{w})$.
The latter and bound (58), yield
$|\zeta_{3}(t)|\leq|\zeta_{3}(0)|+\int_{0}^{t}\mathcal{M}_{tr}\left(|\bm{\eta}_{r}(s)|+\lambda|\zeta_{3}(s)|\right)ds$.
Hence, (60) yields

$|\zeta_{3}(t)|\leq|\zeta_{3}(0)|+\lambda\mathcal{M}_{tr}d_{r}^{+}t+\lambda\mathcal{M}_{tr}\int_{0}^{t}|\zeta_{3}(s)|ds\mbox{.}$

By using the Grönwall lemma once again and assumption (54), which implies, _a
fortiori_ , $\exp(\lambda\mathcal{M}_{tr}t)<2$ for all $t\in[0,T_{tr}]$, we
obtain

$|\zeta_{3}(t)|\leq
2\left[|\zeta_{3}(0)|+\lambda\mathcal{M}_{tr}d_{r}^{+}\pi/\omega\right]\mbox{.}$

Hence, by setting $|\zeta_{3}(0)|\leq d_{3}^{-}$, a sufficient condition for
the property (57) to hold true, is

$d_{3}^{-}\leq d_{3}^{+}/4,\qquad\lambda d_{r}^{+}\leq\omega
d_{3}^{+}/(4\pi\mathcal{M}_{tr})\mbox{.}$ (62)

It is easy to realise that the second conditions appearing in (61) and in (62)
hold in the region of the plane $(d_{r}^{+},d_{3}^{+})$ given by
$d_{3}^{+}\in[4\omega^{-1}\pi\mathcal{M}_{tr}\lambda
d_{r}^{+},\pi^{-1}K_{\bm{w}}^{-1}\omega d_{r}^{+}]$, for all
$d_{r}^{+}\in[0,1]$. Such a region is correctly defined by virtue of (54).
Hence, for any $d_{r}^{+}$, we choose the medium point of the above mentioned
interval, that is,

$d_{3}^{+}:=\Theta d_{r}^{+}\mbox{,}$ (63)

where $\Theta:=2^{-1}(\omega^{-1}\lambda
4\pi\mathcal{M}_{tr}+\pi^{-1}K_{\bm{w}}^{-1}\omega)$. Note that $\Theta<1$
(as, in particular, $(2\pi)^{-1}K_{\bm{w}}^{-1}\omega<1$ by construction).
This implies $d_{3}^{+}<d_{r}^{+}$. Hence, as
$|\bm{\eta}(0)|\geq\max\\{|\bm{\eta}_{r}(0)|,|\eta_{3}(0)|\\}$, the first
conditions of (61) and (62) hold if one requires

$\lambda(d_{3}^{+}/4)\geq|\bm{\eta}(0)|\mbox{.}$ (64)

On the other hand, bounds (57) and (60) imply, for all $t\leq T_{tr}$,

$|\bm{\eta}(t)|\leq|\bm{\eta}_{r}(t)|+\lambda|\zeta_{3}(t)|\leq\lambda\left[1+\Theta\right]d_{r}^{+}\mbox{.}$
(65)

In conclusion, given $\rho_{tr}^{+}>0$, we have that (65) and (64), imply that
by defining

$\rho_{tr}^{-}:=[4(1+\Theta)]^{-1}\Theta\rho_{tr}^{+}\mbox{,}$ (66)

and choosing $|\bm{\eta}(0)|\leq\lambda\rho_{tr}^{-}$ as in (55), the proof is
complete. ∎

\begin{overpic}[width=455.24408pt]{TransferMzVariation.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(4.0,67.5){\footnotesize$u_{3}$}
\put(63.0,24.0){\footnotesize$t$}
\put(41.0,67.0){\footnotesize$\lambda_{(1)}$}
\put(50.0,61.0){\footnotesize$\lambda_{(2)}$}
\put(50.0,50.5){\footnotesize$\lambda_{(3)}$} \end{overpic}

\begin{overpic}[width=455.24408pt]{TransferCurrent.pdf}
\put(0.0,78.0){\footnotesize(b)} \put(4.0,67.5){\footnotesize$\tilde{\beta}$}
\put(63.0,24.0){\footnotesize$t$} \end{overpic}

\begin{overpic}[width=455.24408pt]{Transfer3D.pdf}
\put(0.0,78.0){\footnotesize(c)} \put(13.0,67.0){\footnotesize$\bm{s}^{-}$}
\put(27.0,71.0){\footnotesize$\bm{0}$}
\put(44.0,75.0){\footnotesize$\bm{s}^{+}$}
\put(2.0,58.0){\footnotesize$u_{3}$} \put(17.0,28.0){\footnotesize$u_{2}$}
\put(57.0,28.0){\footnotesize$u_{1}$} \end{overpic}

\begin{overpic}[width=455.24408pt]{TransferStereo.pdf}
\put(0.0,78.0){\footnotesize(d)} \put(5.0,66.0){\footnotesize$w_{1}$}
\put(57.0,24.0){\footnotesize$w_{2}$}
\put(15.0,69.0){\footnotesize$\bm{w}^{-}$}
\put(41.0,69.0){\footnotesize$\bm{0}$}
\put(62.0,69.0){\footnotesize$\bm{w}^{+}$} \end{overpic}

Figure 3: A case study to show the quasi-latitudinal features of the
constructed solutions during the “transfer” stage (the “expulsuion” stage is
not reported for the sake of clarity). In panel (a) the behaviour of
$u_{3}(t)$ for the values $\lambda_{(1)}=0.0055$, $\lambda_{(2)}=0.005$ and
$\lambda_{(3)}=0.001$. The remaining panels show the case with
$\lambda=\lambda_{(3)}$. Panel (b) shows the current realising the control,
whilst panels (c) and (d) represent the trajectory of the system controlled
with such a current in the space $\bm{u}$ and stereographic projection
$\bm{w}$, respectively. In particular, $\bm{w}^{\pm}$ represent the
stereographic projection of $\bm{s}^{\pm}$. The remaining parameters are
reported in Appendix C.

##  7 Attraction

As anticipated, the aim of this section is to study the dynamics in a
neighbourhood of the equilibrium $\bm{s}^{+}$ when the current is switched
off, i.e. $\beta=0$, and show that $\bm{s}^{+}$ is an attractive point for the
system, providing in addition an estimate for the basin of attraction. For
this purpose, we shall proceed in a non-perturbative fashion and consider
system (9). This is motivated by the fact that the attractive behaviour is
characterised by $\alpha$, whilst $\lambda$ does not play a particularly
relevant role in this case. As usual, we shall start by considering the
standard translation

$\bm{u}=\bm{U}+\bm{s}^{+}\mbox{,}$ (67)

and setting $\beta=0$ then $h_{2}=:-D_{2,1}\Omega$ with $\Omega>0$, and
$\gamma:=\sqrt{1-\Omega^{2}}$, in such a way (5) reads as

$\dot{\bm{U}}=\bm{\mathcal{G}}_{1}(\bm{U})+\alpha\bm{\mathcal{G}}_{2}(\bm{U})\mbox{,}$
(68)

where
$\bm{\mathcal{G}}_{1}(\bm{U}):=(D_{3,1}U_{3}(U_{2}-\Omega)-D_{3,2}U_{2}U_{3},-D_{3,1}U_{3}(\gamma+U_{1}),D_{2,1}U_{2}(\gamma+U_{1}))^{\top}$
and $\bm{\mathcal{G}}_{2}$ is defined as a consequence. We shall denote with
$\Phi_{a}^{t}$ the corresponding phase flow. Moreover, let us define
$\tilde{\Psi}(\bm{U}):=\Psi(\bm{u})|_{\bm{u}=\bm{U}+\bm{s}^{+}}$.
In this setting, we can state the following

######  Lemma 7.1.

Define the following quadratic function

$\mathfrak{W}(U_{2},U_{3}):=2^{-1}(D_{2,1}U_{2}^{2}+D_{3,1}U_{3}^{2})$ (69)

and suppose

$3\Omega^{2}\geq 2D_{2,1}/D_{3,1},\qquad
16\sqrt{D_{2,1}/D_{3,1}}\leq\gamma\mbox{.}$ (70)

Then, for all $\alpha>0$, and all $\mathcal{E}:=1+\delta_{a}$, with

$|\delta_{a}|\leq(\gamma/4)^{2}\mbox{,}$ (71)

all the solutions of (68) starting in the set

$\mathcal{B}_{\mathcal{E}}:=\\{\bm{U}\in\mathbb{R}^{3}:|U_{1}|\leq\gamma/4,\quad\mathfrak{W}(U_{2},U_{3})\leq(D_{2,1}\gamma)^{2}/(32D_{3,1}\Omega^{2})\\}\cap\\{\tilde{\Psi}(\bm{U})=\mathcal{E}\\}\mbox{,}$

satisfy the following property

$\bm{U}_{\infty}(\delta_{a}):=\lim_{t\rightarrow+\infty}\Phi_{a}^{t}\left(\mathcal{B}_{\mathcal{E}}\right)=(-\gamma+\sqrt{\gamma^{2}+\delta_{a}},0,0)\mbox{,}$
(72)

and

$|\bm{U}_{\infty}(\delta_{a})|\leq\delta_{a}/(2\gamma)\mbox{.}$ (73)

In particular, those trajectories with $\mathcal{E}=1$ are asymptotic to
$\bm{U}=\bm{0}$ (and hence $\bm{u}(t)\rightarrow\bm{s}^{+}$).
Furthermore, the projection of $\mathcal{B}_{\mathcal{E}}$ on the plane
$(U_{2},U_{3})$ contains the disk centred at the origin of radius

$r_{sm}:=(4|\Omega|D_{3,1})^{-1}\gamma D_{2,1}\mbox{.}$ (74)

######  Remark 7.1.

The importance of these last two lines of statement lies in the fact that, by
construction, the CQL solution which realises the transfer, ends up at a point
located in the vicinity of $\bm{U}=(0,0,-K)$. This will be made precise later
(see, for instance, (79)). However, it is important to stress for the moment
that this is necessary and we shall choose $2K\leq r_{sm}$ later on, in order
to ensure that such a point is suitably attracted once the current is switched
off.

A numerical validation of Lem. 7.1 is reported in fig. 4.

###### Proof.

The proof relies on the Theory of Lyapunov functions and related tools for the
basin of attraction estimation, see e.g. [Kha02] for a comprehensive
description. The proof will be carried out in the variables $U_{2},U_{3}$,
being the behaviour of the first variable a consequence, by the conservation
law $\Psi(\bm{u})=\operatorname{const.}$, see (7).
Let us start noticing that (69) satisfies $\mathfrak{W}(0,0)=0$ and
$\mathfrak{W}(U_{2},U_{3})>0$ for all $(U_{2},U_{3})\neq\bm{0}$. Furthermore,
it is immediate to check that its derivative along the solutions of the
“undamped” (68), i.e. with $\alpha=0$, satisfies

$\dot{\mathfrak{W}}|_{\alpha=0}:=\nabla_{(U_{2},U_{3})}\mathfrak{W}\cdot\bm{\mathcal{G}}_{1}(\bm{U})=0\mbox{,}$

hence $\mathfrak{W}$ is a (non-strict) Lyapunov function for the undamped
(partial) system, showing that it is stable. However, our aim is to show that
the damped system satisfies $\dot{\mathfrak{W}}<0$ on a suitable set, proving
in this way the (stronger) asymptotic stability.
As $D_{2,1}\ll D_{3,1}$ by assumption, the region we are considering to prove
the negative definiteness of $\dot{\mathfrak{W}}$ will be conveniently chosen
as non-isotropic. More precisely, we shall proceed by introducing in
$\dot{\mathfrak{W}}$ the following variables transformation

$\bm{U}=(\mu x,\varepsilon(\gamma D_{2,1})^{-1}\cos\theta,\varepsilon
D_{3,1}^{-1}\sin\theta)\mbox{,}$

with $(x,\theta)\in[0,1]\times[0,2\pi]$ and $\varepsilon,\mu>0$ to be
determined. The latter gives

$\displaystyle\dot{\tilde{\mathfrak{W}}}$
$\displaystyle=-\alpha\varepsilon^{2}\left[1+2\mu
x(\sin^{2}\theta+\gamma^{-1}\cos^{2}\theta)+(\mu
x)^{2}(\sin^{2}\theta+\gamma^{-2}\cos^{2}\theta)+\right.$ (75)
$\displaystyle\left.-2\Omega\varepsilon\gamma^{-1}(D_{2,1}^{-1}-D_{3,1}^{-1})\cos\theta\sin^{2}\theta+\varepsilon^{2}\gamma^{-2}(D_{2,1}D_{3,1})^{-2}(D_{3,1}-D_{2,1})^{2}\sin^{2}\theta\cos^{2}\theta\right]\mbox{,}$

where $\tilde{\mathfrak{W}}$ stands for $\mathfrak{W}$ in the new set of
variables. It is evident that for sufficiently small $\varepsilon,\mu$ and all
$\alpha>0$, the latter is strictly negative. Our goal is now to find a bound
for the thresholds $\varepsilon_{0},\mu_{0}$ in such a way this property
persists for all $\varepsilon\leq\varepsilon_{0}$ and $\mu\leq\mu_{0}$.
Let us firstly deal with $O(\varepsilon^{2})$ coefficients. For this purpose,
it is immediate to check that

$1+2\mu x(\sin^{2}\theta+\gamma^{-1}\cos^{2}\theta)+(\mu
x)^{2}(\sin^{2}\theta+\gamma^{-2}\cos^{2}\theta)\geq 1-2\mu/\gamma\mbox{.}$

Hence the latter will be, say, greater than $1/2$, for all
$\mu\leq\mu_{0}:=\gamma/4$.
As for the terms of $O(\varepsilon^{3})$, by recalling (6), one has

$(D_{2,1}^{-1}-D_{3,1}^{-1})\cos\theta\sin^{2}\theta\leq D_{2,1}^{-1}\mbox{.}$

Finally, $O(\varepsilon^{4})$ contributions are clearly strictly positive. In
conclusion one has

$\dot{\tilde{\mathfrak{W}}}<\alpha\varepsilon^{2}\left[2\Omega\varepsilon/(D_{2,1}\gamma)-1/2\right]\mbox{.}$

The latter implies that the derivative along the solutions
$\tilde{\mathfrak{W}}$ will be strictly negative by choosing

$\varepsilon_{0}:=D_{2,1}\gamma/(4\Omega)\mbox{.}$

This implies that, in the original set of variables $\bm{U}$, the cylinder
ellipse - shaped region on which $\dot{\mathfrak{W}}$ is strictly negative, is
given by

$\mathcal{C}_{el}:=\\{\bm{U}\in\mathbb{R}^{3}:|U_{1}|\leq\gamma/4,\quad
U_{2}^{2}+(U_{3}D_{3,1}/(D_{2,1}\gamma))^{2}=(4\Omega)^{-2}\\}\mbox{.}$

In order to obtain the required estimate for the basin of attraction, we need
to find the largest value $\mathfrak{W}^{*}$ for which the level curve
$\mathfrak{W}(U_{2},U_{3})=\mathfrak{W}^{*}$ is entirely contained in the
projection of $\mathcal{C}_{el}$ on the $(U_{2},U_{3})$ plane. This is a
simple problem of constrained optimisation in which $\mathfrak{W}$ has to be
minimised on the ellipse
$U_{2}^{2}+(U_{3}D_{3,1}/(D_{2,1}\gamma))^{2}=(4\Omega)^{-2}$. It is easy to
check that the required minimum is attained at

$(U_{2}^{*},U_{3}^{*}):=(0,\pm D_{2,1}\gamma/(4D_{3,1}\Omega))\mbox{,}$

implying
$\mathfrak{W}^{*}:=\mathfrak{W}(U_{2}^{*},U_{3}^{*})\equiv(D_{2,1}\gamma)^{2}/(32D_{3,1}\Omega^{2})$.
It is now necessary to check that $|U_{1}(t)|\leq\gamma/4$ for all $t$,
ensuring that $\bm{U}(t)\in\mathcal{C}_{el}$ and then $\dot{\mathfrak{W}}<0$
all $t$.
For this purpose, let us observe that the level
$\tilde{\Psi}(\bm{U})=\mathcal{E}$ is equivalent to

$(U_{1}+\gamma)^{2}+(U_{2}-\Omega)^{2}+U_{3}^{2}=1+\delta_{a}$

The latter defines a $\delta_{a}$-family of surfaces in which $U_{1}$ can be
written as a graph (surface) over the other two variables in a suitable
neighbourhood of the origin, as follows

$U_{1}=U_{1}(U_{2},U_{3};\delta_{a})\equiv-\gamma+\sqrt{1+\delta_{a}-(U_{2}-\Omega)^{2}-U_{3}^{2}}\mbox{,}$
(76)

and such that $U_{1}(0,0;0)=0$. The sign has been chosen according to the fact
that we are describing the portion of sphere located near $\bm{s}^{+}$. The
task consists in showing that, under the required condition for $\delta_{a}$,
any point $(U_{2},U_{3})\in\mathcal{B}_{\mathcal{E}}$, yields
$|U_{1}|\leq\gamma/4$. Given the structure of (76), it will be sufficient to
show $U_{1}\geq-\gamma/4$, being condition $U_{1}\leq\gamma/4$ a consequence.
From (76), one has $U_{1}\geq-\gamma+\sqrt{-|\delta_{a}|+\gamma^{2}-A_{a}}$,
where $A_{a}:=2\Omega|U_{2}|+U_{2}^{2}+U_{3}^{2}$. It is immediate to verify
that, by choosing $\delta_{a}$ as in (71), the desired condition is obtained
if

$A_{a}\leq(3/8)\gamma^{2}\mbox{.}$ (77)

On the other hand, in $\mathcal{B}_{\mathcal{E}}$ one has

$|U_{2}|\leq[\gamma/(4\Omega)]q_{a},\qquad
U_{2}^{2}+U_{3}^{2}\leq[\gamma/(4\Omega)]^{2}q_{a}^{2}\mbox{,}$

where $q_{a}:=\sqrt{D_{2,1}/D_{3,1}}$. This implies that condition (77) holds
if

$4q_{a}+2^{-1}\Omega^{-2}q_{a}^{2}\gamma\leq\gamma\mbox{.}$

A sufficient condition for the latter to hold is that, for instance,
$4q_{a}\leq\gamma/4$ and $2^{-1}\Omega^{-2}q_{a}^{2}\leq 3/4$. They are
equivalent to the first and the second of (70), respectively. Equation (72)
follows directly from (76) whilst bound (73) is immediate from (72) and a
Taylor estimate.
As for the last statement it is sufficient to observe that the ellipse
resulting from the projection of $\mathcal{B}_{\mathcal{E}}$ clearly contains
the disk whose radius is not larger than the minor semi-axis of the former,
which is given by the r.h.s. of (74). ∎

\begin{overpic}[width=455.24408pt]{Attractivity3D_E1.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(3.0,59.0){\footnotesize$U_{3}$}
\put(10.0,30.0){\footnotesize$U_{2}$} \put(60.0,26.0){\footnotesize$U_{1}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Attractivity3D_E1p1.pdf}
\put(0.0,78.0){\footnotesize(b)} \put(3.0,59.0){\footnotesize$U_{3}$}
\put(10.0,30.0){\footnotesize$U_{2}$} \put(60.0,26.0){\footnotesize$U_{1}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Attractivity2D.pdf}
\put(0.0,78.0){\footnotesize(c)} \put(3.0,60.0){\footnotesize$U_{3}$}
\put(51.0,24.0){\footnotesize$U_{2}$} \end{overpic}

\begin{overpic}[width=455.24408pt]{AttractivityLiap.pdf}
\put(0.0,78.0){\footnotesize(d)} \put(3.0,67.0){\footnotesize$\mathfrak{W}$}
\put(53.0,24.0){\footnotesize$t$} \end{overpic}

Figure 4: The field attractivity feature within the basin
$\mathcal{B}_{\mathcal{E}}$ when the injected current $\beta(t)$ is switched
off. The (black) ellipse in the first three panels represents the set
$\mathfrak{W}(U_{2},U_{3})=(D_{2,1}\gamma)^{2}/(32D_{3,1}\Omega^{2})$ (recall
that $\bm{U}$, defined via (67), represents a coordinate system in which the
equilibrium $\bm{s}^{+}$ corresponds to the origin). Panels (a) and (b) show a
trajectory starting on the boundary of $\mathcal{B}_{\mathcal{E}}$ for
$\delta_{a}=0$ and $\delta_{a}=0.1$, respectively. In panel (c) the projection
on the $(U_{2},U_{3})$ plane of the latter trajectory is depicted. The last
panel shows, as expected, the decreasing behaviour of the Lyapunov function
$\mathfrak{W}$ along the same trajectory. The full set of parameters is
reported in Appendix C.

#### Proof of Thm. 4.1

###### Proof.

Let us consider the class of expulsion targets and currents such that
$\beta_{e},K^{2}=O(\lambda)$. In order to use Lemmata 5.1, 6.1 and 7.1 we are
going to make the assumptions required by these results. First of all, given
$D_{2,1}=O(\lambda)$ and $D_{3,1}=O(1)$, let us choose $\Omega$ within the
limitation prescribed by the first of (70), this determines a range for
$h_{2}=-D_{2,1}\Omega$. It is easy to realise that $\Omega=O(\sqrt{\lambda})$.
As a consequence, $\gamma\sim 1-O(\lambda)$, hence we can assume without loss
of generality, that the last of (70) is satisfied as well. There is a certain
amount of freedom in the choice of $D_{3,2}$, provided that (6) holds. We
anticipate that the value of $\delta_{a}$ will be chosen, but only at the very
end of the proof, according to (71), hence we shall proceed by supposing that
such a condition is satisfied.
Once the possibility to use Lem. 7.1 has been guaranteed, let us examine the
setting of the remaining Lemmata. As for Lem. 5.1, the values of
$r^{*},\mathcal{M}_{1,2}$ are uniquely determined, this yields $\lambda_{e}$.
The same holds for the constants of $\mathcal{M}_{tr},\mathcal{K}_{\bm{w}}$
and the threshold $\lambda_{tr}$ of Lem. 6.1. In order to proceed with the
proof we shall assume $\lambda_{0}:=\min\\{\lambda_{e},\lambda_{tr}\\}$,
although practical applications of Lem. 6.1 may require a further restriction
of such a threshold, as already discussed in Rem. 6.4.
Let us now use the mentioned Lemmata to construct a suitable neighbourhood of
the equilibrium $\bm{s}^{-}$ which is transported in the vicinity of
$\bm{s}^{+}$ and then attracted by the point close to it on the corresponding
energy level. The proof will be complete once (18) is validated and a precise
value for $r^{-}$ is computed.
For this purpose let us start by observing that, from (27) and (19), at the
end of the expulsion stage started in $\bm{u}_{c}(0)\equiv\bm{s}^{-}$, we get

$\bm{u}_{c}(T_{e})=(-\gamma+\lambda\kappa\bar{a},-\Omega-\lambda\kappa\bar{b},-K)\mbox{,}$

with $\kappa:=(2\beta_{e})^{-1}(-1+\sqrt{1-2K^{2}})\sim
K^{2}/(2\beta_{e})=O(1)$ by assumption and $T_{e}$ given by (26). Note that
$T_{e}\sim K/\beta_{e}=O(\lambda^{-1/2})$. Furthermore,
$\Phi(\bm{u}_{c}(T_{e}))=1+O(\lambda)$.
Let us now set $\bm{w}_{r}(0):=\bm{u}_{c}(T_{e})$ in such a way to start the
transfer stage. This determines the value of $T_{r}$ via (50). By construction
we have
$w_{1}^{[0]}(T_{tr})=-\mathcal{A}_{m}-K^{2}=\gamma+O(\lambda)=:\hat{\gamma}$
and $w_{3}^{[0]}(T_{tr})=-K$. In order to determine $w_{2}^{[0]}(T_{tr})$, one
can use the conservation law (35), which reads as

$\sigma^{-1}(\lambda\kappa\bar{a}-\gamma)^{2}+\sigma(\Omega+\lambda\kappa\bar{b})^{2}=\sigma^{-1}\hat{\gamma}^{2}+\sigma(w_{2}^{[0]}(T_{r}))^{2}\mbox{.}$

The latter gives, recalling (21) and expanding in $\lambda$,

$w_{2}^{[0]}(T_{r})=-\Omega+O(\lambda)\mbox{.}$ (78)

The sign in the r.h.s. of the latter has been chosen by using the fact that
before and on the first occurrence of $w_{1}^{[0]}=-\mathcal{A}_{m}-K^{2}$,
the trajectory lies in the $w_{2}^{[0]}$-negative half-plane. Hence, by
recalling (47), we have

$\bm{w}(T_{e})=(\gamma,-\Omega,-K)+O(\lambda)\mbox{.}$ (79)

Once more, the usual conservation law is easily checked as
$\Phi(\bm{w}(T_{e}))=1+O(\lambda)$.
Now we need to make sure that (79) is actually contained in the basin of
attraction $\mathcal{B}_{\mathcal{E}}$ in such a the attractivity property
stated in Lem. 7.1 can be used. For this purpose, we recall the very last
statement of Lem. 7.1 and set $r_{sm}\geq(5/4)K$ (i.e. slightly bigger than
$K$ itself). Hence, (74) yields

$\Omega^{2}\leq(5KD_{3,1})^{2}(\gamma D_{2,1})^{2}\mbox{,}$

yielding, by recalling that $\gamma^{2}=1-\Omega^{2}$,
$\Omega^{2}\leq[1+25(D_{3,1}/D_{2,1}^{2})K^{2}]^{-1}$.
The latter, if compared with the first of (70) (which provides a lower bound
for $\Omega^{2}$), makes sense provided that $50K^{2}D_{3,1}\leq
3\gamma^{2}D_{2,1}$. One could choose, for instance,

$K=\bar{K}:=(\gamma/4)\sqrt{D_{2,1}/D_{3,1}}\mbox{,}$

then computing $\Omega^{2}$ as the mid-point given by the two conditions

$\Omega^{2}:=2^{-1}[(2/3)(D_{2,1}/D_{3,1})+\gamma^{2}/(25\bar{K}^{2})(D_{2,1}/D_{3,1})^{2}]=(52/75)(D_{2,1}/D_{3,1})\mbox{,}$

this determines $h_{2}$.
Hence, any suitable manifold $\mathcal{E}=1+\delta_{a}$ contained in the
cylinder $\\{\sqrt{U_{2}^{2}+U_{3}^{2}}\leq(5/4)K,\,\,|U_{1}|\leq\gamma/4\\}$
is contained in $\mathcal{B}_{\mathcal{E}}$. This implies that any
$\mathfrak{B}_{K/8}(\bm{w}(T_{e}))$ is contained in
$\mathcal{B}_{\mathcal{E}}$ as well. Hence, we can set
$\rho_{tr}^{+}:=K/(8\lambda)$ and
$\lambda\rho_{tr}^{-}:=[32(1+\Theta)]^{-1}\Theta K$ by (66), then use Lem.
6.1. The latter, allows us to conclude that the set
$\mathfrak{B}_{\lambda\rho_{tr}^{-}}(\bm{w}(T_{e}))$ will evolve inside
$\mathcal{B}_{\mathcal{E}}$ under the action of the controlled phase flow
$\Phi_{tr}^{t}$.
As $\omega=O(\lambda^{-1/2})$ because of the choice for $K$, we have
$\Theta=O(\sqrt{\lambda})$, hence $\lambda\rho_{tr}^{-}=O(\lambda)$. It is now
sufficient to recall Lem. 5.1 and set $\rho_{e}:=\rho_{tr}^{-}$ to determine
the value of $r^{-}$, which is easily found via (30) as

$r^{-}=[128\mathcal{M}_{e}(1+\Theta)]^{-1}\Theta K\mbox{.}$

Similarly, we notice that $r^{-}=O(\lambda)$.
As (73) holds for the family of points starting in $\mathcal{B}_{\mathcal{E}}$
and parameterised by $\delta_{a}$ in the admissible range (71), this property
will be true, _a fortiori_ , for the points of the set

$\mathcal{U}_{a}:=\Phi_{tr}^{T_{tr}}\left(\Phi_{e}^{T_{e}}(\mathfrak{B}_{r^{-}}(\bm{s}^{-}))\right)\subset\mathcal{B}_{\mathcal{E}}\mbox{.}$

Hence, by defining

$\bm{u}_{\infty}:=\lim_{t\rightarrow+\infty}\Phi_{a}^{t}\left(\bigcup_{\mathcal{E}=1+\delta_{a}}\mathcal{B}_{\mathcal{E}}\right)\mbox{,}$

its distance from $\bm{s}^{+}$ is readily bounded by $\delta_{a}$ via (73) as
follows

$|\bm{u}_{\infty}(\delta_{a})-\bm{s}^{+}|=|\bm{U}_{\infty}(\delta_{a})|\leq\delta_{a}/(2\gamma)\mbox{.}$

On the other hand, by the conservation of $\Psi(\bm{u})$, an upper bound for
$\delta_{a}$ is found by evaluating the maximum of $|\Psi(\bm{u})-1|$ for all
$\bm{u}:=\bm{s}^{-}+\bm{u}^{\prime}\in\mathfrak{B}_{r^{-}}(\bm{s}^{-})$ i.e.
for all $\bm{u}^{\prime}\in\mathfrak{B}_{r^{-}}(\bm{0})$. This is
straightforward, as $\Psi(\bm{s}^{-}+\bm{u}^{\prime})\leq 1+r^{-}$ by the
triangle inequality, hence $|\delta_{a}|\leq r^{-}$. This proves (18) with
$f:=(2\gamma)^{-1}$. The proof of Thm. 4.1 is now complete. ∎

\begin{overpic}[width=455.24408pt]{Fullnew_Fig5_u.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(25.0,60.0){\footnotesize$\bm{s}^{-}$}
\put(44.0,58.0){\footnotesize$\bm{0}$} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Fullnew_Fig5_curr.pdf}
\put(0.0,78.0){\footnotesize(b)}
\put(16.0,69.75){\footnotesize$\longleftarrow$}
\put(5.0,65.0){\footnotesize$\beta(t)$} \put(52.0,23.0){\footnotesize$t$}
\end{overpic}

Figure 5: Example of a full switching process. Panel (a) shows the (piece-
wise) trajectory in the space $\bm{u}$, as a sequence of an “expulsion”
(black), a “transfer” (red) and finally the “attractive” stage (blue). The
initial condition has been chosen as
$\bm{u}(-T_{e})=\bm{s}^{-}+\lambda(-0.1,0.05,0)$, with $T_{e}=1.0262$. The
first $\sim 133$ time units of the constructed controlling current are
reported in panel (b). The small black segment for $t\in[-T_{e},0]$
(emphasised by the arrow), represents the current injected during the
expulsion stage. We have set $\lambda=0.002$. The full set of parameters is
reported in Appendix C.

##  8 Validation tests and conclusions

The aim of this section is to provide examples of full switching processes
realised according to Thm. 4.1. To summarise, after an initial “expulsion”
stage, the evolving point is “transferred” via a CQL solution in the estimated
basin of attraction of the target equilibrium so that it can be “attracted” by
it. An example of a full switching is reported in fig. 5.
As anticipated in sec. 3, the approach proposed here is conceptually different
from the well established ballistic switching procedure, see, e.g. [dPS+15a].
From the qualitative viewpoint, these stages could be compared to our
“expulsion” and “attraction” respectively, see fig. 6, panel (a). However, as
it is evident from the latter, in order to successfully complete the
switching, the operator needs to trust a basin of attraction which is, for
instance, way larger than the one rigorously computed in this work. Clearly,
the one stated in Lem. 7.1 is nothing but a sufficient condition but it
suggests that, if not addressed with further (and highly specialised) tools
and ad-hoc arguments, the attractivity of larger sets retains no more than a
probabilistic validity. On the contrary, the proposed CQL strategy has a
genuinely deterministic character, provided that $\lambda$ is sufficiently
small (according to the bounds described in the proof) and the initial set of
initial conditions is chosen in $\mathfrak{B}_{r^{-}}(\bm{s}^{-})$.
In addition, as anticipated in sec. 3, our method only requires a minimal
amount of initial injected current during the “expulsion” stage, as the
solution is subsequently “guided” towards $\mathcal{B}_{\mathcal{E}}$ via the
constructed control, along a CQL trajectory. As a result, this remarkably
reduces the “energy” carried by $u_{3}$ and then the amplitude of the
oscillations around $\bm{s}^{+}$ when the current is switched off. See fig. 6
for a depiction of this phenomenon in a comparison with the ballistic
approach.
As a further check in relation with the robustness of this method and its
deterministic character, the outcomes of a stress test are proposed in fig. 7,
in which some errors are simulated with respect to the proposed control. More
precisely, either a time dilatation-contraction coefficient for the expulsion
time $T_{e}\leftarrow jT_{e}$ with $j=1\pm 0.02$ or for the transfer control,
i.e. ${\beta}_{tr}(t)\leftarrow\tilde{\beta}_{tr}(jt)$ are considered.

In conclusion, the present work proposes an analytical formulation apt to
realise a fully deterministic switching mechanism. The argument relies on the
concept of CQL solution, a highly non-local object at the heart of the
procedure, constructed via a perturbative approach which exploits intrinsic
technologically relevant features of the system. The mentioned CQL solutions
are realised by the explicit determination of a time dependent control,
determined by using well established tools of Hamiltonian perturbation theory.
Interestingly, the latter possess the potential to increase the range of
validity of the method to a even larger class of systems, should a higher
order analysis be considered. The numerical experiments proposed offer either
a validation or a visual interpretation of the main statement of Thm. 4.1 and
the related Lemmata.

\begin{overpic}[width=455.24408pt]{Fullnew_Fig6_u_panela.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Fullnew_Fig6_u_panelb.pdf}
\put(0.0,78.0){\footnotesize(b)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Fullnew_Fig6_u1_panelc.pdf}
\put(0.0,78.0){\footnotesize(c)} \put(4.0,66.0){\footnotesize$u_{1}$}
\put(54.0,24.0){\footnotesize$t$} \end{overpic}

\begin{overpic}[width=455.24408pt]{Fullnew_Fig6_u1_paneld.pdf}
\put(0.0,78.0){\footnotesize(d)} \put(4.0,66.0){\footnotesize$u_{1}$}
\put(54.0,24.0){\footnotesize$t$} \end{overpic}

Figure 6: Comparison between the ballistic switching (a) and the one obtained
by using CQL solutions (b). The initial condition has been set as
$\bm{u}(-T_{e})=\bm{s}^{-}+(\lambda,0,0)$ in both cases. Below, ((c) and (d),
respectively), the corresponding behaviour of the first coordinate $u_{1}(t)$.
A considerably small residual oscillation of the mentioned coordinate can be
observed in the CQL case. The remaining parameters used are specified in
Appendix C.

\begin{overpic}[width=455.24408pt]{Full_stress_minus2percent_texp.pdf}
\put(0.0,78.0){\footnotesize(a)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Full_stress_plus2percent_texp.pdf}
\put(0.0,78.0){\footnotesize(b)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Full_stress_minus2percent_coefft.pdf}
\put(0.0,78.0){\footnotesize(c)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

\begin{overpic}[width=455.24408pt]{Full_stress_plus2percent_coefft.pdf}
\put(0.0,78.0){\footnotesize(d)} \put(2.0,58.0){\footnotesize$u_{3}$}
\put(17.0,28.0){\footnotesize$u_{1}$} \put(57.0,28.0){\footnotesize$u_{2}$}
\end{overpic}

Figure 7: Switching trajectories via the CQL approach with a simulated error
in the control: panels (a) and (b) show the case of a time
contraction/dilatation of the computed “espulsion” time, more precisely
$T_{e}\leftarrow 0.98T_{e}$ and $T_{e}\leftarrow 1.02T_{e}$, respectively.
Panels (c) and (d) report an analogous test during the “transfer” stage
instead. Hence, the time variable in the control $\tilde{\beta}_{tr}(t)$ is
replaced with $t\leftarrow 0.98t$ and $t\leftarrow 1.02t$, respectively.
Taking into account the “considerable” duration of the transfer process (it
has an order of $10^{2}$ time units) the latter should be regarded as a pretty
“hard” test and evident trajectory modifications such as those shown in (c)
and (d), are reasonably expected. Full set of parameters in Appendix C.

## Appendix A

Explicit expressions of $\bm{\mathcal{F}}(\bm{u};\bm{v})$ and
$\bm{\mathcal{R}}(\bm{u};\bm{v})$:

$\displaystyle\mathcal{F}_{1}$
$\displaystyle=-u_{1}u_{3}(-v_{2}^{2}\tilde{\alpha}D_{3,2}K-v_{1}^{2}\tilde{\alpha}D_{3,1}K+v_{1}\tilde{h}_{2}-v_{1}v_{2}\tilde{D}_{2,1})\rho+u_{3}\tilde{h}_{2}+u_{1}u_{3}^{2}\tilde{\alpha}D_{3,1}$
$\displaystyle\mathcal{F}_{2}$
$\displaystyle=u_{2}u_{3}^{2}\tilde{\alpha}D_{3,2}-u_{2}u_{3}(-v_{2}^{2}\tilde{\alpha}D_{3,2}K-v_{1}^{2}\tilde{\alpha}D_{3,1}K+v_{1}\tilde{h}_{2}-v_{1}v_{2}\tilde{D}_{2,1})\rho$
$\displaystyle\mathcal{F}_{3}$
$\displaystyle=(u_{2}^{2}+u_{1}^{2})(-v_{2}^{2}\tilde{\alpha}D_{3,2}K-v_{1}^{2}\tilde{\alpha}D_{3,1}K+v_{1}\tilde{h}_{2}-v_{1}v_{2}\tilde{D}_{2,1})\rho-
u_{1}\tilde{h}_{2}$
$\displaystyle+\tilde{\alpha}(-u_{2}^{2}u_{3}D_{3,2}-u_{1}^{2}u_{3}D_{3,1})+u_{1}u_{2}\tilde{D}_{2,1}$
$\displaystyle\mathcal{R}_{1}$
$\displaystyle=u_{1}v_{2}u_{3}\tilde{\alpha}\tilde{h}_{2}K\rho+\tilde{\alpha}(u_{1}u_{2}^{2}\tilde{D}_{2,1}-u_{1}u_{2}\tilde{h}_{2})$
$\displaystyle\mathcal{R}_{2}$
$\displaystyle=u_{2}v_{2}u_{3}\tilde{\alpha}\tilde{h}_{2}K\rho+\tilde{\alpha}(u_{3}^{2}\tilde{h}_{2}+u_{1}^{2}\tilde{h}_{2}-u_{1}^{2}u_{2}\tilde{D}_{2,1})$
$\displaystyle\mathcal{R}_{3}$
$\displaystyle=-(u_{2}^{2}+u_{1}^{2})v_{2}\tilde{\alpha}\tilde{h}_{2}K\rho-
u_{2}u_{3}\tilde{\alpha}\tilde{h}_{2}$

## Appendix B

Explicit expression of the entries of
$\\{{\bm{D}}\bm{V}(\bm{\xi},\tilde{\beta})\\}$. Let us denote them with
$\lambda b_{i,j}$. We have

$\displaystyle b_{1,1}$
$\displaystyle=\xi_{3}\tilde{\beta}\lambda+\tilde{\alpha}(\tilde{D}_{2,1}\xi_{2}\Omega+D_{3,1}\xi_{3}^{2})\lambda^{2}+\tilde{\alpha}\tilde{D}_{2,1}\xi_{2}^{2}\lambda^{3}$
$\displaystyle b_{1,2}$
$\displaystyle=D_{3,2}\xi_{3}-\tilde{\alpha}\tilde{D}_{2,1}\Omega\tilde{h}_{2}\lambda+\tilde{\alpha}\tilde{D}_{2,1}(\xi_{1}\Omega-2\xi_{2}\tilde{h}_{2})\lambda^{2}+2\tilde{\alpha}\tilde{D}_{2,1}\xi_{1}\xi_{2}\lambda^{3}$
$\displaystyle b_{1,3}$
$\displaystyle=(-\tilde{\beta}\tilde{h}_{2}+\tilde{D}_{2,1}\Omega+D_{3,2}\xi_{2})+(\xi_{1}\tilde{\beta}-2\tilde{\alpha}D_{3,1}\xi_{3}\tilde{h}_{2})\lambda+2\tilde{\alpha}D_{3,1}\xi_{1}\xi_{3}\lambda^{2}$
$\displaystyle b_{2,1}$
$\displaystyle=-D_{3,1}\xi_{3}+2\tilde{\alpha}\tilde{D}_{2,1}\xi_{2}\tilde{h}_{2}\lambda^{2}-2\tilde{\alpha}\tilde{D}_{2,1}\xi_{1}\xi_{2}\lambda^{3}$
$\displaystyle b_{2,2}$
$\displaystyle=(\xi_{3}\tilde{\beta}-\tilde{\alpha}\tilde{D}_{2,1}\gamma^{2})\lambda+\tilde{\alpha}(2\tilde{D}_{2,1}\xi_{1}\tilde{h}_{2}+D_{3,2}\xi_{3}^{2})\lambda^{2}-\tilde{\alpha}\tilde{D}_{2,1}\xi_{1}^{2}\lambda^{3}$
$\displaystyle b_{2,3}$ $\displaystyle=\tilde{\beta}\Omega-
D_{3,1}\xi_{1}+(2\tilde{\alpha}D_{3,2}\xi_{3}\Omega+\xi_{2}\tilde{\beta})\lambda+2\tilde{\alpha}\xi_{3}(\tilde{D}_{2,1}\Omega+D_{3,2}\xi_{2})\lambda^{2}$
$\displaystyle b_{3,1}$
$\displaystyle=2\tilde{\beta}\tilde{h}_{2}(2\tilde{\alpha}D_{3,1}\xi_{3}\tilde{h}_{2}-2\xi_{1}\tilde{\beta}+\tilde{D}_{2,1}\xi_{2})\lambda-2\tilde{\alpha}D_{3,1}\xi_{1}\xi_{3}\lambda^{2}$
$\displaystyle b_{3,2}$
$\displaystyle=-(\tilde{D}_{2,1}\tilde{h}_{2}+2\tilde{\beta}\Omega)+(-2\tilde{\alpha}D_{3,2}\xi_{3}\Omega-2\xi_{2}\tilde{\beta}+\tilde{D}_{2,1}\xi_{1})\lambda-\tilde{\alpha}\xi_{3}(\tilde{D}_{2,1}\Omega+2D_{3,2}\xi_{2})\lambda^{2}$
$\displaystyle b_{3,3}$
$\displaystyle=-\tilde{\alpha}(D_{3,1}\tilde{h}_{2}^{2}+D_{3,2}\Omega^{2})+\tilde{\alpha}(2D_{3,1}\xi_{1}\tilde{h}_{2}-\tilde{D}_{2,1}\Omega^{2}-2D_{3,2}\xi_{2}\Omega)\lambda$
$\displaystyle-\tilde{\alpha}(\tilde{D}_{2,1}\xi_{2}\Omega+D_{3,2}\xi_{2}^{2}+D_{3,1}\xi_{1}^{2})\lambda^{2}$

## Appendix C

In this section the numerical values used in the experiments are reported. All
the experiments have the following common features

  * •

$D_{1}=0.0411$, $D_{3}=0.8527$ (see, for instance, [dPS+15a]),

  * •

$D_{2}:=D_{1}+6.51\lambda$ (perturbative setting),

  * •

$\gamma=\sqrt{1-\Omega^{2}}$.

The remaining parameters have been chosen as follows

Fig. | $D_{2}$ | $\tilde{\alpha}$ | $\tilde{\beta}_{e}$ | $\Omega$ | $\tilde{h}_{2}$ | $K$
---|---|---|---|---|---|---
2 | 0.1127 | 2 | 3 | 0.0676 | -0.4400 | 0.0700
3 | 0.1127 | 2 | 3 | 0.0676 | -0.4400 | 0.0800
4 | 0.0802 | 4 | $\times$ | 0.1799 | -1.1708 | 0.0524
5 | 0.0802 | 4 | 3 | 0.1036 | -0.6741 | 0.0308
6 | 0.0802 | 2 | 3 | 0.1036 | -0.6741 | 0.0308
7 | 0.0802 | 2 | 3 | 0.1036 | -0.6741 | 0.0308

### Acknowledgements

This work has been supported by the Italian Ministry of University and
Research, PRIN2020 funding program, grant number 2020PY8KTC.
The numerical simulations and the corresponding plots have been performed with
GNU Octave [EBHW20], whilst [Max23] has been used for the algebraic
manipulations.

## References

  * [ADS16] Dmytro Apalkov, Bernard Dieny, and Jon M Slaughter.  Magnetoresistive random access memory.  Proceedings of the IEEE, 104(10):1796–1830, 2016.
  * [Ber09] V. Berdichevsky.  Variational Principles of Continuum Mechanics: I. Fundamentals.  Interaction of Mechanics and Mathematics. Springer Berlin Heidelberg, 2009\.
  * [BLS+10] D Bedau, H Liu, JZ Sun, JA Katine, EE Fullerton, S Mangin, and AD Kent.  Spin-transfer pulse switching: From the dynamic to the thermally activated regime.  Applied Physics Letters, 97(26), 2010.
  * [BMSd03] Giorgio Bertotti, Isaak D Mayergoyz, Claudio Serpico, and Massimiliano d’Aquino.  Geometrical analysis of precessional switching and relaxation in uniformly magnetized bodies.  IEEE transactions on magnetics, 39(5):2501–2503, 2003.
  * [BSH+17] Sabpreet Bhatti, Rachid Sbiaa, Atsufumi Hirohata, Hideo Ohno, Shunsuke Fukami, and SN Piramanayagam.  Spintronics based random access memory: a review.  Materials Today, 20(9):530–548, 2017.
  * [BSM13] G Bertotti, Claudio Serpico, and ID Mayergoyz.  Probabilistic aspects of magnetization relaxation in single-domain nanomagnets.  Physical review letters, 110(14):147205, 2013.
  * [DBB+20] Thibaut Devolder, Olivier Bultynck, Paul Bouquin, Van Dai Nguyen, Siddharth Rao, D Wan, B Sorée, IP Radu, GS Kar, and S Couet.  Back hopping in spin transfer torque switching of perpendicularly magnetized tunnel junctions.  Physical Review B, 102(18):184406, 2020.
  * [DPG+20] Bernard Dieny, Ioan Lucian Prejbeanu, Kevin Garello, Pietro Gambardella, Paulo Freitas, Ronald Lehndorff, Wolfgang Raberg, Ursula Ebels, Sergej O Demokritov, Johan Akerman, et al.  Opportunities and challenges for spintronics in the microelectronics industry.  Nature Electronics, 3(8):446–459, 2020.
  * [dPS+15a] M. d’Aquino, S. Perna, C. Serpico, G. Bertotti, and I. D. Mayergoyz.  Analysis of reliable sub-ns spin-torque switching under transverse bias magnetic fields.  Journal of Applied Physics, 117(17):17B716, 03 2015.
  * [dPS+15b] Massimiliano d’Aquino, S Perna, C Serpico, G Bertotti, and ID Mayergoyz.  Analysis of reliable sub-ns spin-torque switching under transverse bias magnetic fields.  Journal of Applied Physics, 117(17), 2015.
  * [dSS+04] M d’Aquino, Werner Scholz, Thomas Schrefl, Claudio Serpico, and Josef Fidler.  Numerical and analytical study of fast precessional switching.  Journal of applied physics, 95(11):7055–7057, 2004.
  * [EBHW20] J.W. Eaton, D. Bateman, S. Hauberg, and R. Wehbring.  GNU Octave version 5.2.0 manual: a high-level interactive language for numerical computations, 2020.
  * [FW16] A. Fortunati and S. Wiggins.  Integrability and strong normal forms for non-autonomous systems in a neighbourhood of an equilibrium.  Journal of Mathematical Physics, 57(9):092703, 2016.
  * [FW19] A. Fortunati and S. Wiggins.  A Lie transform approach to the construction of Lyapunov functions in autonomous and non-autonomous systems.  Journal of Mathematical Physics, 60(8):082704, 08 2019.
  * [Gio03] A. Giorgilli.  Exponential stability of Hamiltonian systems.  In Dynamical systems. Part I, Pubbl. Cent. Ric. Mat. Ennio Giorgi, pages 87–198. Scuola Norm. Sup., Pisa, 2003.
  * [HYN+20] Atsufumi Hirohata, Keisuke Yamada, Yoshinobu Nakatani, Ioan-Lucian Prejbeanu, Bernard Diény, Philipp Pirro, and Burkard Hillebrands.  Review on spintronics: Principles and device applications.  Journal of Magnetism and Magnetic Materials, 509:166711, 2020.
  * [Kha02] H.K. Khalil.  Nonlinear Systems.  Pearson Education. Prentice Hall, 2002.
  * [KR02] Shehzaad Kaka and Stephen E Russek.  Precessional switching of submicrometer spin valves.  Applied Physics Letters, 80(16):2958–2960, 2002.
  * [Max23] Maxima.  Maxima, a computer algebra system. version 5.47.0, 2023.
  * [MBS09] I.D. Mayergoyz, G. Bertotti, and C. Serpico.  Nonlinear Magnetization Dynamics in Nanosystems.  ISSN. Elsevier Science, 2009.
  * [MSB+09] T Min, JZ Sun, R Beach, D Tang, and P Wang.  Back-hopping after spin torque transfer induced magnetization switching in magnetic tunneling junction cells.  Journal of Applied Physics, 105(7), 2009.
  * [SdBM09] Claudio Serpico, Massimiliano d’Aquino, Giorgio Bertotti, and Isaak D Mayergoyz.  Analytical description of quasi-random magnetization relaxation to equilibrium.  IEEE transactions on magnetics, 45(11):5224–5227, 2009.
  * [Slo96] John C Slonczewski.  Current-driven excitation of magnetic multilayers.  Journal of Magnetism and Magnetic Materials, 159(1-2):L1–L7, 1996\.
  * [VKL+22] Sergei M Vlasov, Grzegorz J Kwiatkowski, Igor S Lobanov, Valery M Uzdin, and Pavel F Bessarab.  Optimal protocol for spin-orbit torque switching of a perpendicular nanomagnet.  Physical Review B, 105(13):134404, 2022.