# Poster Abstract: Hierarchical Subchannel Allocation for Mode-3 Vehicle-to-
Vehicle Sidelink Communications

Luis F. Abanto-Leon  Eindhoven University of TechnologyEindhovenNetherlands
<EMAIL_ADDRESS>,  Arie Koppelaar  NXP SemiconductorsEindhovenNetherlands
<EMAIL_ADDRESS>and  Sonia Heemstra de Groot  Eindhoven University of
TechnologyEindhovenNetherlands<EMAIL_ADDRESS>
(2017)

###### Abstract.

In this poster we present a graph-based hierarchical subchannel allocation
scheme for V2V sidelink communications in Mode-3. Under this scheme, the
eNodeB allocates subchannels for in-coverage vehicles. Then, vehicles will
broadcast directly without the eNodeB intervening in the process. Therefore,
in each communications cluster, it will become crucial to prevent allocation
conflicts in time domain since vehicles will not be able to transmit and
receive simultaneously. We present a solution where the time-domain
requirement can be enforced through vertex aggregation. Additionally,
allocation of subchannels is performed sequentially from the most to the least
allocation-constrained cluster. We show through simulations that the proposed
approach attains near-optimality.

resource allocation; mode-3 V2V; sidelink

††copyright: acmlicensed††ccs: Networks Network resources
allocation††journalyear: 2017††copyright: acmlicensed††conference: SenSys ’17;
November 6–8, 2017; Delft, Netherlands††price: 15.00††doi:
10.1145/3131672.3136987††isbn: 978-1-4503-5459-2/17/11

##  1\. Motivation and Contributions

In V2V Mode-3, eNodeBs assign subchannels to vehicles in order for them to
periodically broadcast CAM messages (b2, ). A crucial aspect is to ensure that
vehicles in the same cluster will broadcast in orthogonal time subchannels111A
subchannel is a time-frequency resource chunk capable of sufficiently
conveying a CAM message. to avoid conflicts. In general, resource/subchannel
allocation problems can be represented as weighted bipartite graphs. However,
in this scenario there is an additional time orthogonality constraint which
cannot be straightforwardly handled by conventional graph matching methods
(b3, ). Thus, in our approach the mentioned constraint has been taken into
account. We also perform the allocation task in a sequential manner based on
the constrainedness of each cluster. To illustrate the gist of the problem, in
Fig. 1 we show two partially overlapping clusters where a conflict between
vehicles $V_{8}$ and $V_{10}$ is generated as the allotted subchannels are in
the same subframe.

$\tiny K$$\tiny L$ Figure 1. V2V broadcast communications scenario

##  2\. Proposed Approach

The objective is to find an assignment of subchannels to vehicles such that
the system capacity is maximized while respecting the constraints that prevent
conflicts. This allocation problem is formulated as

(1a) $\displaystyle\begin{aligned} &{\rm max}~{}{\bf c}^{T}{\bf x}\\\ &{\rm
subject~{}to}~{}{\left[\begin{array}[]{c}{\bf I}_{N\times N}\otimes{\bf
1}_{1\times L}\\\ \hline\cr{\bf U}_{J\times N}\otimes{\bf I}_{L\times
L}\end{array}\right]}\otimes{\bf 1}_{1\times K}~{}{\bf x}={\bf
1}\end{aligned}$

where ${\bf c}\in\mathbb{R}^{M},{\bf x}\in\mathbb{B}^{M}$ with $M=NLK$. The
amount of clusters is $J$ (each consisting of $N_{j}$ vehicles,
$j=1,\dots,J$); the total number of vehicles is $N$; $L$ is the number of time
subframes and $K$ represents the number of subchannels per subframe. ${\bf U}$
is the membership matrix which depicts the association of vehicles to the
different clusters. Instead of approaching (1) optimally through exhaustive
search, we can solve the allocation in a hierarchical and sequential manner
for each cluster $\mathcal{V}^{\text{(}j\text{)}}$ while retaining the
solutions from previous allocations. Thus, we solve several subproblems in
order of constrainedness, which will lead to a suboptimal solution. Each
subproblem can be modeled as a weighted bipartite graph
$G^{\text{(}j\text{)}}=\big{(}\mathcal{V}^{\text{(}j\text{)}},\mathcal{R},\mathcal{E}^{\text{(}j\text{)}}\big{)}$
as shown in Fig. 2, where vehicles and subchannels are modeled as vertices.
The subchannels $r_{k}$ in $\mathcal{R}$ have been assembled into $L$ disjoint
groups $\\{\mathcal{R}_{l}\\}_{l=1}^{L}$ called macro-vertices in order to
manage the time-domain constraints (vertex aggregation). The edge weight
$c^{\text{(}j\text{)}}_{ik}=B\log_{2}\big{(}1+\mathsf{SINR}^{\text{(}j\text{)}}_{ik}\big{)}$
represents the achievable capacity of vehicle
$v^{\text{(}j\text{)}}_{i}\in\mathcal{V}^{\text{(}j\text{)}}$ in subchannel
$r_{k}\in\mathcal{R}$. Thus, the following formulation maximizes the capacity
per cluster

(2a) $\displaystyle\begin{aligned} &{\rm max}~{}{\bf c}^{T}_{j}{\bf x}_{j}\\\
&{\rm subject~{}to}~{}{\left[\begin{array}[]{c}{\bf I}_{N_{j}\times
N_{j}}\otimes{\bf 1}_{1\times L}\\\ \hline\cr{\bf 1}_{1\times
N_{j}}\otimes{\bf I}_{L\times L}\end{array}\right]}\otimes{\bf 1}_{1\times
K}~{}{\bf x}_{j}={\bf 1}.\end{aligned}$

For completeness, we assume that $N_{1}=\dots=N_{j}=L$, thus
$\mathbf{x}_{j}=\text{[}x^{\text{(}j\text{)}}_{1,1},\dots,x^{\text{(}j\text{)}}_{L,KL}\text{]}^{T}$
and
$\mathbf{c}_{j}=\text{[}c^{\text{(}j\text{)}}_{1,1},\dots,c^{\text{(}j\text{)}}_{L,KL}\text{]}^{T}$.
Also, we claim without a proof due to space limitations, that (2) can be
recast as (3) without affecting optimality.

(3) $\begin{array}[]{lclcl}&&{\rm max}~{}{\bf d}^{T}_{j}{\bf y}_{j}\\\ &&{\rm
subject~{}to}~{}{\left[\begin{array}[]{c}{\bf I}_{L\times L}\otimes{\bf
1}_{1\times L}\\\ \hline\cr{\bf 1}_{1\times L}\otimes{\bf I}_{L\times
L}\end{array}\right]}{\bf y}_{j}={\bf 1}.\end{array}$

where ${\bf
d}_{j}=\lim_{\beta\to\infty}\frac{1}{\beta}\overset{\begin{subarray}{c}\circ\end{subarray}}{\log}\Big{\\{}({\bf
I}_{M\times M}\otimes{\bf 1}_{1\times K})\mathrm{e}^{\circ\beta{\bf
c}_{j}}\Big{\\}}$ with
$\overset{\begin{subarray}{c}\circ\end{subarray}}{\log}\\{\cdot\\}$ and
$\mathrm{e}^{\circ\\{\cdot\\}}$ representing the element-wise natural
logarithm and Hadamard exponential (b4, ), respectively. Note that ${\bf
y}_{j}$ is of lower dimensionality than ${\bf x}_{j}$ and therefore it can be
solved with less complexity.

$v^{\text{(}j\text{)}}_{1}$$v^{\text{(}j\text{)}}_{2}$⋮$v^{\text{(}j\text{)}}_{N_{j}}$$r_{1}$$r_{2}$$...$$r_{K}$$r_{K+1}$$r_{K+2}$$...$$r_{2K}$⋮$r_{K(L-1)+1}$$r_{K(L-1)+2}$$...$$r_{KL}$macro-
macro-macro-vertex $\mathcal{R}_{1}$vertex $\mathcal{R}_{2}$vertex
$\mathcal{R}_{L}$ $\mathcal{V}^{\text{(}j\text{)}}$ Vehicles $\mathcal{R}$
Subchannels Figure 2. Constrained weighted bipartite graph

##  3\. Simulations

We consider a 10 MHz channel which is divided into subchannels, each with
dimensions of 1ms in time and 1.26 MHz in frequency (b5, ). In our model, we
assume a CAM message rate of 10 Hz. In Fig. 3, we compare 4 different
algorithms in base of the average over 1000 simulations. We considered both
overlapping and non-overlapping clusters, each with at most $N=100$ vehicles.
Through simulations, we show that our scheme can attain near-optimality as its
performance is similar to exhaustive search. Considering the system average
rate criterion, our approach has an advantage over the greedy algorithm. Also,
when considering the worst-rate vehicle, our proposal excels as it is capable
of providing a higher level of fairness. In all cases, the random allocation
algorithm is outperformed by the other approaches.

Highest-Rate Vehicle System Average Rate Worst-Rate Vehicle Rate Standard
Deviation
$0$$5$$10$$8.97$$8.24$$7.14$$0.48$$8.97$$8.11$$6.94$$0.55$$8.97$$7.88$$5.35$$0.65$$8.17$$4.55$$1.43$$1.21$Rate
[Mbits / s / subchannel] Exhaustive Search Proposed Algorithm Greedy Algoritm
Random Algorithm Figure 3. Vehicles data rate
$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$$1$$2$$3$$4$$5$$6$$7$$8$Number of
vehicles ($N$)Rate [Mbits / s / subchannel] Exhaustive Search Proposed
Algorithm Greedy Algoritm Random Algorithm Figure 4. Worst-rate vehicle

Fig. 4 shows the achievable rate for the worst-rate vehicle. We observe that
when the vehicle density per cluster is low, the greedy approach attains near
optimal solutions as there are more subchannels than vehicles to serve.
However, as the density increases, especially near the overload state, its
performance dramatically drops. Nevertheless, our approach can provide a fair
allocation even in high vehicle density scenarios.

##  4\. Conclusion

We have presented a novel hierarchical subchannel allocation scheme for V2V
sidelink communications with intra-cluster conflict constraints. The approach
is capable of attaining near-optimality at less complexity than exhaustive
search.

## References

  * (1) Intelligent transport systems (its); stdma recommended parameters and settings for cooperative its; access layer part.  Technical Report ETSI TR 102 861, ETSI, January 2012.
  * (2) James Munkres.  Algorithms for the assignment and transportation problems.  SIAM J Appl Math, 5(1):32–38, March 1957.
  * (3) Fumio Hiai.  Monotonicity for entrywise functions of matrices.  Journal of Linear Algebra and its Applications, 431(8):1125–1146, September 2009.
  * (4) Technical specification group radio access network; evolved universal terrestrial radio access (e-utra); physical layer procedures; (release 14) v14.2.0.  Technical Report 3GPP TS 36.213, 3GPP, March 2017.


∎

11institutetext: S. Breiter 22institutetext: Astronomical Observatory
Institute, Faculty of Physics, Adam Mickiewicz University, Sloneczna 36,
60-286 Poznan, Poland
22email<EMAIL_ADDRESS>33institutetext: K. Langner 44institutetext:
Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz
University, Sloneczna 36, 60-286 Poznan, Poland
44email<EMAIL_ADDRESS>
# Kustaanheimo-Stiefel transformation with an arbitrary defining vector

S. Breiter  K. Langner

(Received: 7 November 2016 / Accepted: date)

###### Abstract

Kustaanheimo-Stiefel (KS) transformation depends on the choice of some
preferred direction in the Cartesian 3D space. This choice, seldom explicitly
mentioned, amounts typically to the direction of the first or the third
coordinate axis in celestial mechanics and atomic physics, respectively. The
present work develops a canonical KS transformation with an arbitrary
preferred direction, indicated by what we call a defining vector. Using a mix
of vector and quaternion algebra, we formulate the transformation in a
reference frame independent manner. The link between the oscillator and
Keplerian first integrals is given. As an example of the present formulation,
the Keplerian motion in a rotating frame is re-investigated.

###### Keywords:

KS variables Kepler problem quaternions regularization

††journal: CELE

##  1 Introduction

There are many ways to convert the Kepler problem into the isotropic harmonic
oscillator. A comprehensive summary can be found in the monograph by Cordani
(2003) and in the work of Deprit et al (1994), the two sources overlapping
only partially. Among all the methods, the Kustaanheimo-Stiefel (KS)
transformation stands exceptional by its simplicity, popularity, and history
(traced back to Heinz Hopf, or even Carl Friedrich Gauss – see Volk (1976)).
The literature concerning theoretical and practical aspects of the KS
variables is already vast, yet there is a feature which, to our knowledge, has
not brought enough attention by either being too obvious, or by not being
realized. We mean the existence of some preferred direction in the definition
of the KS variables. Kustaanheimo, Stiefel, and most of their followers in the
realm of celestial mechanics, use the set of variables so designed, that only
the first Cartesian coordinate $x_{1}$ involves the squares of the KS
coordinates. In the present paper, we will refer to it as KS1. However, in
atomic physics a different set, to be named KS3, is considered standard, at
least since late 1970s (e.g. Duru and Kleinert, 1979). There, only the $x_{3}$
involves the squares. The choice of the preferred direction is of marginal
significance for the unperturbed Keplerian problem, but it may either
simplify, or complicate expressions resulting from added perturbations. Thus
we have found it worthwhile to establish a general KS transformation with the
preferred direction left unspecified.

In section 2 we give a brief outline of how the KS1 transformation has settled
down in the realm of celestial mechanics. Not attempting a detailed
bibliographic survey, we mark five turning points: i) the first paper of
Kustaanheimo (1964) where the transformation is born in spinor form, ii) the
metamorphosis into the L-matrix setup done by Kustaanheimo and Stiefel (1965),
iii) early quaternion formulation by Vivarelli (1983), iv) the refinement of
the quaternionic form due to Deprit et al (1994), and v) interpretation of the
KS variables as rotation parameters (Saha, 2009). Save for the last point, the
reader may observe how the special role of the $Ox_{1}$ axis direction is
transferred from one setup to another, and how the presence of a preferred
direction becomes more and more evident.

Once the stage has been set, we discuss the main theme in Section 3.
Introducing the notion of a defining vector, we build the KS transformation
(and its canonical extension) with an arbitrary preferred direction. Working
within the general quaternion and vector formalism, we refrain from using
explicit expressions in terms of coordinates. Hopefully, it should allow
capturing the intrinsic, coordinate independent features of the KS
transformation.

Section 4 considers the general properties of the Hamiltonian in KS variables.
Special attention is paid to linking the invariants of the unperturbed problem
in its two incarnations: the Kepler problem in Cartesian coordinates, and the
isotropic oscillator in KS variables.

Finally, in order to give an example that the choice of the defining vector
does matter in a perturbed problem, we return to the Kepler problem in the
rotating reference frame. Compared to earlier works which used the KS1 set,
the task can be considerably facilitated by the appropriate selection of the
preferred direction, as demonstrated in section 5.

##  2 KS1 transformation

###  2.1 The roots

What Paul Kustaanheimo announced at the Oberwolfach conference on Mathematical
Methods of Celestial Mechanics and published the same year (Kustaanheimo,
1964) is worth a brief recall, because neither Annales Universitatis
Turkuensis, nor the Publications of the Astronomical Observatory of Helsinki
are widespread enough. Moreover, an awkward notation has masked some features
that emerge immediately when the more common conventions are applied.

Given a Cartesian position vector
$\mbox{\boldmath$x$}=(x_{1},x_{2},x_{3})^{\mathrm{T}}$, Kustannheimo extended
it formally to a null 4-vector of the Minkowski space, using the length
$r=\sqrt{\mbox{\boldmath$x$}\cdot\mbox{\boldmath$x$}}=x_{0}$ as an extra
coordinate. These served as the coefficients for a linear combination of a
unit matrix ${\mathbf{\sigma}}_{0}$ and the Pauli matrices (Cartan, 1966)

$\sigma_{1}=\left(\begin{array}[]{rr}0&1\\\
1&0\end{array}\right),\quad\sigma_{2}=\left(\begin{array}[]{rr}0&-\mathrm{i}\\\
\mathrm{i}&0\end{array}\right),\quad\sigma_{3}=\left(\begin{array}[]{rr}1&0\\\
0&-1\end{array}\right).$ (1)

The result is a complex matrix

$\mathbf{S}=r\sigma_{0}+x_{1}\sigma_{3}+x_{2}\sigma_{1}-x_{3}\sigma_{2}=\left(\begin{array}[]{cc}r+x_{1}&x_{2}+\mathrm{i}\,x_{3}\\\
x_{2}-\mathrm{i}\,x_{3}&r-x_{1}\end{array}\right),$ (2)

assigned to the vector $(x_{0},x_{1},x_{2},x_{3})^{\mathrm{T}}$. The
difference between the treatment of $x_{1}$ in the diagonal and the complex
pair $x_{2}\pm\mathrm{i}x_{3}$ can be spotted already at this stage.111The
disparity of indices associated with $x_{i}$ and $\sigma_{j}$, was originally
not visible, since Kustaanheimo used a different set of the Pauli matrices
(namely: $\mathbf{i}_{x}=\sigma_{3}$, $\mathbf{i}_{y}=\sigma_{1}$, and
$\mathbf{i}_{z}=-\sigma_{2}$). It took some time until the physicists swapped
to $\mathbf{S}=\sum_{j=0}^{3}x_{j}\sigma_{j}$, the KS3 convention.

Any Hermitian matrix $\mathbf{S}$ can be expressed in terms of two complex
numbers $S_{1}$ and $S_{2}$. Using a complex 2-vector
$\mbox{\boldmath$s$}=(S_{1},S_{2})^{\mathrm{T}}$, Kustaanheimo (1964) defined
$\mathbf{S}$ in terms of its Hermitian outer product

$\mathbf{S}=2\left(\begin{array}[]{c}S_{1}\\\ S_{2}\\\
\end{array}\right)\left(\begin{array}[]{cc}\overline{S}_{1}&\overline{S}_{2}\\\
\end{array}\right)=2\left(\begin{array}[]{cc}S_{1}\overline{S}_{1}&S_{1}\overline{S}_{2}\\\
\overline{S}_{1}S_{2}&S_{2}\overline{S}_{2}\end{array}\right).$ (3)

By the equivalence of (2) and (3), the vector $s$ becomes a rank 1 spinor,222
Intriguingly, although the initial name of the KS transformation was ‘the
spinor regularization’, Kustaanheimo (1964) uses the word ‘spinor’ only twice:
once in the title, and once in the abstract. Kustaanheimo and Stiefel (1965)
use it only once – in the title. and the matrix $\mathbf{S}$ – a rank 2 spinor
associated with $x$ (c.f. Bellandi Filho and Menon, 1987; Steane, 2013).
Equating respective elements of (2) and (3), one readily finds

$\displaystyle x_{1}$ $\displaystyle=$ $\displaystyle
S_{1}\overline{S}_{1}-S_{2}\overline{S}_{2},$ $\displaystyle x_{2}$
$\displaystyle=$ $\displaystyle\overline{S}_{1}S_{2}+S_{1}\overline{S}_{2},$
(4) $\displaystyle x_{3}$ $\displaystyle=$
$\displaystyle\mathrm{i}\,\left(\overline{S}_{1}S_{2}-S_{1}\overline{S}_{2}\right),$
$\displaystyle r$ $\displaystyle=$ $\displaystyle
S_{1}\overline{S}_{1}+S_{2}\overline{S}_{2}.$

As noted by Kustaanheimo, the transformation is not unique; indeed it involves
only the products with conjugates, so using any spinor
$\mbox{\boldmath$q$}=(Q_{1},Q_{2})^{\mathrm{T}}=\mbox{\boldmath$s$}\exp{\mathrm{i}\phi}$,
leads to the same result. Geometrically, it means that any pair of complex
numbers resulting from rotations of $S_{1}$ and $S_{2}$ on a complex plane by
the same angle, generates the same position vector $x$ in (4)

The regularization of the Keplerian motion does not appear until the
coordinate transformation (4) is augmented by the time transformation of the
Sundman type. There, Kustaanheimo proposed a general formula relating the
pseudo-time $\tau$ to its physical counterpart $t$

$\frac{\mathrm{d}\tau}{\mathrm{d}t}=\frac{\beta}{r}\,\exp{\int
K\,\mathrm{d}t},$ (5)

where $K$ could be an arbitrary function of position, velocity and time. This
flexibility has never been seriously explored, and the simplest choice of
$K=0$ has became standard. The regularization converts the Kepler problem with
energy constant $h$ into a spinor oscillator problem

$\mbox{\boldmath$s$}^{\prime\prime}=\frac{h}{2}\mbox{\boldmath$s$},$ (6)

yet, this simple form does not show up, until the bilinear constraint

$\overline{S}_{1}S^{\prime}_{1}-\overline{S}^{\prime}_{1}S_{1}+\overline{S}_{2}S^{\prime}_{2}-\overline{S}^{\prime}_{2}S_{2}=0,$
(7)

is imposed,333Another wording of condition (7), provided by Kustaanheimo
(1964), is: $\overline{S}_{1}S^{\prime}_{1}+\overline{S}_{2}S^{\prime}_{2}$,
being a half of $r^{\prime}$, is a real quantity. where the prime stands for
the derivative with respect to the Sundman time $\tau$. Kustaanheimo justified
this choice by observing the invariance of the left hand side of (7) in the
perturbed Kepler problem with a particular form of perturbation (linear in
coordinates, velocities and angular momentum).

One can only speculate what would be the fate of the Kustaanheimo’s
transformation, has it not attracted the attention of Eduard Stiefel who
coauthored the paper published next year (Kustaanheimo and Stiefel, 1965). In
the new mise-en-scéne, the complex variables and unnecessary generalization
were dropped, and the discussion focused on a real 4-vector
$\mbox{\boldmath$u$}=(u_{1},u_{2},u_{3},u_{4})^{\mathrm{T}}$, containing the
parameters of the substitution

$S_{1}=u_{1}+\mathrm{i}\,u_{4},\qquad S_{2}=u_{2}-\mathrm{i}\,u_{3}.$ (8)

Then, eq. (4) takes the form

$\displaystyle x_{1}$ $\displaystyle=$ $\displaystyle
u_{1}^{2}-u_{2}^{2}-u_{3}^{2}+u_{4}^{2},$ $\displaystyle x_{2}$
$\displaystyle=$ $\displaystyle 2(u_{1}u_{2}-u_{3}u_{4}),$ (9) $\displaystyle
x_{3}$ $\displaystyle=$ $\displaystyle 2(u_{1}u_{3}+u_{2}u_{4}),$
$\displaystyle r$ $\displaystyle=$ $\displaystyle
u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+u_{4}^{2},$

and the constraint (7) is turned into

$u_{4}u_{1}^{\prime}-u_{3}u_{2}^{\prime}+u_{2}u_{3}^{\prime}-u_{1}u_{4}^{\prime}=0.$
(10)

This derivation, however, cannot be found in the paper; Kustaanheimo and
Stiefel (1965) have burnt the bridge leading back to the 1964 work and started
the presentation from the matrix equation that related $u$ with a 4-vector
$\mbox{\boldmath$x$}=(x_{1},x_{2},x_{3},0)^{\mathrm{T}}$ through the matrix
product

$\mbox{\boldmath$x$}=\mathbf{L}(\mbox{\boldmath$u$})\,\mbox{\boldmath$u$}.$
(11)

The L-matrix definition, given by Kustaanheimo and Stiefel (1965) with an
intriguing clause ‘for example’, was

$\mathbf{L}(\mbox{\boldmath$u$})=\left(\begin{array}[]{rrrr}u_{1}&-u_{2}&-u_{3}&u_{4}\\\
u_{2}&u_{1}&-u_{4}&-u_{3}\\\ u_{3}&u_{4}&u_{1}&u_{2}\\\
u_{4}&-u_{3}&u_{2}&-u_{1}\end{array}\right).$ (12)

This definition leads to the first three equations (9) for $x_{j}$, whereas
the last of equations (9) had been postulated as a required property of
$\mathbf{L}(\mbox{\boldmath$u$})$. Noteworthy, the special role of $x_{1}$ has
been conserved, as visible in the first of equations (9): other coordinates
are defined by products $u_{i}u_{j}$, whereas $x_{1}$ is composed of the pure
squares $u_{i}^{2}$.

The work, published in a more widespread journal and written in a manner
friendly to the celestial mechanics audience, considerably helped to
promulgate what is now known as the Kustaanheimo-Stiefel transformation. By
the influence of the Stiefel and Scheifele (1971) monograph, the matrix
approach became paradigmatic in the celestial mechanics community, and the
‘for example’ choice (12) has been taken for granted and obvious, save for
occasional renumbering of indices and the change of sign in $u_{4}$.

###  2.2 Enter quaternions

The close relation between spinors and quaternions was known already to
Cartan. But in the framework of the L-matrix formulation, the relation of the
KS variables to the quaternion algebra is merely an additional aspect,
mentioned by Kustaanheimo and Stiefel (1965) or Stiefel and Scheifele (1971)
as an interesting, but probably unimportant curio. Setting the KS
transformation in the quaternion formalism, originated by Vivarelli (1983) and
applauded by Deprit et al (1994), offered new paths to understanding the known
properties of the transformation, as well as to its generalization to higher
dimensions.

Following Deprit et al (1994) we treat a quaternion
$\mathsf{v}=(v_{0},v_{1},v_{2},v_{3})$, or
$\mathsf{v}=(v_{0},\mbox{\boldmath$v$})$, as a union of a scalar $v_{0}$ and
of a vector $\mbox{\boldmath$v$}=(v_{1},v_{2},v_{3})^{\mathrm{T}}$.
Conjugating a quaternion, we change the signs of its vector part, i.e.

$\bar{\mathsf{v}}=(v_{0},-v_{1},-v_{2},-v_{3})=(v_{0},-\mbox{\boldmath$v$}).$
(13)

Extracting the vector part is performed by means of the operator $\natural$,
so that for $\mathsf{v}=(v_{0},\mbox{\boldmath$v$})$

$\mbox{\boldmath$v$}=\mathsf{v}^{\natural}.$ (14)

The usual scalar product, marked with a dot,

$\mathsf{v}\cdot\mathsf{w}=v_{0}w_{0}+\mbox{\boldmath$v$}\cdot\mbox{\boldmath$w$},$
(15)

is commutative, but the quaternion product

$\mathsf{v}\,\mathsf{w}=\left(v_{0}w_{0}-\mbox{\boldmath$v$}\cdot\mbox{\boldmath$w$},v_{0}\mbox{\boldmath$w$}+w_{0}\mbox{\boldmath$v$}+\mbox{\boldmath$v$}\times\mbox{\boldmath$w$}\right),$
(16)

is not. Using a standard basis

$\mathsf{e}_{0}=(1,\mbox{\boldmath$0$}),\quad\mathsf{e}_{1}=(0,\mbox{\boldmath$e$}_{1}),\quad\mathsf{e}_{2}=(0,\mbox{\boldmath$e$}_{2}),\quad\mathsf{e}_{3}=(0,\mbox{\boldmath$e$}_{3}),$
(17)

we recover the classical ‘$1ijk$’ multiplication rules of Hamilton for the
basis quaternions: $\mathsf{e}_{0}\mathsf{e}_{1}=\mathsf{e}_{1}$ for $1i=i$,
$\mathsf{e}_{1}\mathsf{e}_{2}=\mathsf{e}_{3}$ for $ij=k$, etc. The inverse of
a quaternion is, in full analogy with complex numbers,
$\mathsf{v}^{-1}=\bar{\mathsf{v}}/|\mathsf{v}|^{2}$, where the norm is
$|\mathsf{v}|=\sqrt{\mathsf{v}\cdot\mathsf{v}}$. The conjugate of a product
is, typically for noncommutative operations,

$\overline{\mathsf{u}\mathsf{v}}=\bar{\mathsf{v}}\bar{\mathsf{u}}.$ (18)

Let us observe a useful property of the mixed dot product,

$\bar{\mathsf{u}}\cdot(\mathsf{v}\mathsf{w})=\bar{\mathsf{w}}\cdot(\mathsf{u}\mathsf{v})=\bar{\mathsf{v}}\cdot(\mathsf{w}\mathsf{u}),$
(19)

echoing the mixed product rule of the standard vector algebra.

Another useful operation is called a quaternion outer product (Morais et al,
2014) or a quaternion cross product (Stiefel and Scheifele, 1971; Vivarelli,
1988; Deprit et al, 1994). Conventions vary among the authors; we adopt the
one of Deprit et al (1994)

$\mathsf{u}\wedge\mathsf{v}=\frac{\mathsf{v}\bar{\mathsf{u}}-\mathsf{u}\bar{\mathsf{v}}}{2}=\left(0,\,u_{0}\mbox{\boldmath$v$}-v_{0}\mbox{\boldmath$u$}+\mbox{\boldmath$u$}\times\mbox{\boldmath$v$}\right).$
(20)

A remarkable property of the cross product, not mentioned by Deprit et al
(1994), is a factor exchange rule

$(\mathsf{u}\mathsf{v})\wedge\mathsf{w}=\mathsf{u}\wedge(\mathsf{w}\bar{\mathsf{v}}),$
(21)

following directly from (20) and (18).

Note that the cross product always results in a quaternion with null scalar
part, called a pure vector (Deprit et al, 1994) or, more often, a pure
quaternion (e.g. Morais et al, 2014). We adopt the former convention.

Casting the KS transformation in a quaternion form is not a novelty. It can be
found already in Stiefel and Scheifele (1971). In spite of the discouraging
comments attached by the authors, Vivarelli (1983) returned to this formalism
and issued a different quaternion form of the transformation. Feeling obliged
to adhere to the ‘for example’ convention of Kustaanheimo and Stiefel (1965),
she reconstructed the transformation (9) as a quaternion product

$(x_{1},x_{2},x_{3},0)=\mathsf{u}\,\left[\mathsf{e}_{3}\,\bar{\mathsf{u}}\,\bar{\mathsf{e}}_{3}\right]=\mathsf{u}\,\mathsf{u}_{\ast},$
(22)

where $\mathsf{u}=(u_{1},u_{2},u_{3},u_{4})$, and the product in square
brackets is an ‘anti-involute’ $\mathsf{u}_{\ast}$ of $\mathsf{u}$ (an
operation, that actually amounts to a trivial change of sign in $u_{4}$)444
The ‘anti-involute’ was later reinvented as a ‘star conjugate’ by Waldvogel
(2006, 2008). His formulation of KS1, similar to (22), offers an interesting
interpretation of the invariant (43). . Of course, the presence of
$\mathsf{e}_{3}$ in (22) does not mean that the direction of $x_{3}$ gains
some special position; the transformation of Vivarelli remains the pure KS1.

Deprit et al (1994), decided to link a more natural, vector type quaternion
$\mathsf{x}=(0,x_{1},x_{2},x_{3})$ with a KS quaternion
$\mathsf{v}=(v_{0},v_{1},v_{2},v_{3})$ and then, not needing an anti-involute,
found that555setting their additional parameter $\alpha=1$.

$\mathsf{x}=\mathsf{v}\,\mathsf{e}_{1}\,\bar{\mathsf{v}},$ (23)

which is not far from the original quaternion formulation of Stiefel and
Scheifele (1971). Converting a Stiefel-Scheifele-Vivarelli quaternion
$\mathsf{u}$ to a Deprit-Elipe-Ferrer quaternion $\mathsf{v}$ can be effected
by the rule

$(v_{0},v_{1},v_{2},v_{3})\rightarrow(-u_{4},u_{1},u_{2},u_{3}).$ (24)

With this rule, the outcome of (23) is equivalent to the transformation of
Kustaanheimo and Stiefel (1965). In particular, the distinguished role of
$x_{1}$ remains unaffected, and clearly marked by the presence of
$\mathsf{e}_{1}$ in equation (23), as it will follow from the interpretation
given below.

###  2.3 Kustaanheimo-Stiefel meet Euler-Rodrigues

For practitioners, the most enjoyed property of quaternions is their
straightforward relation to rotation. Rotation of a vector $x$, formally
treated as a quaternion with null scalar part $(0,\mbox{\boldmath$x$})$, is
specified by a unit quaternion $\mathsf{q}$, including the complete
information about the rotation angle $0\leq\theta\leq\pi$, and rotation axis
given by the unit vector $n$. Then, with

$\mathsf{q}=\left(\cos\frac{\theta}{2},\sin\frac{\theta}{2}\,\mbox{\boldmath$n$}\right),$
(25)

the rotated vector $y$ is obtained through the quaternion product

$(0,\mbox{\boldmath$y$})=\mathsf{q}\,(0,\mbox{\boldmath$x$})\,\bar{\mathsf{q}}.$
(26)

If we write the unit quaternion appearing in (26) simply as
$\mathsf{q}=(q_{0},q_{1},q_{2},q_{3})$, its components are the Euler-Rodrigues
parameters (sometimes called the Cayley parameters) of the rotation matrix.
Indeed, skipping the scalar part, equation (26) may be rewritten in the matrix
form as

$\mbox{\boldmath$y$}=\mathbf{R}(\mathsf{q})\mbox{\boldmath$x$},$ (27)

with the rotation matrix

$\mathbf{R}(\mathsf{q})=\left(\begin{array}[]{rrr}q_{0}^{2}+q_{1}^{2}-q_{2}^{2}-q_{3}^{2}&2\,(q_{1}q_{2}-q_{0}q_{3})&2\,(q_{0}q_{2}+q_{1}q_{3})\\\
2\,(q_{1}q_{2}+q_{0}q_{3})&\quad
q_{0}^{2}-q_{1}^{2}+q_{2}^{2}-q_{3}^{2}&-2\,(q_{0}q_{1}-q_{2}q_{3})\\\
-2\,(q_{0}q_{2}-q_{1}q_{3})&2\,(q_{0}q_{1}+q_{2}q_{3})&\quad
q_{0}^{2}-q_{1}^{2}-q_{2}^{2}+q_{3}^{2}\\\ \end{array}\right).$ (28)

This property reveals the meaning of the quaternion KS transformation (23):

up to the reassignment (24), the normalized Kustaanheimo-Stiefel variables
$u_{i}/\sqrt{r}$ are the Euler-Rodrigues parameters of the rotation turning
the unit vector of the first axis $\mbox{\boldmath$e$}_{1}$ into the unit
radius vector $\mbox{\boldmath$x$}/r$.

Thus we find another argument in favor of the claim that the KS1
transformation attaches a special role to the axis $Ox_{1}$.

The existence of some relation between KS variables and rotation was mentioned
‘for the record’ by Stiefel and Scheifele (1971), who declared the lack of
interest in studying it closer. Then Vivarelli (1983) returned to this issue,
but her description is based on a statement that since a unit quaternion
$\mathsf{q}$ ‘represents a rotation’, so a product
$\mathsf{q}\mathsf{q}_{\ast}$ also ‘represents’ some rotation with the axis
and angle expressions provided. But, unlike (23), the assignment
$\mathsf{x}=\mathsf{q}\mathsf{q}_{\ast}$ is not a formula for rotation of some
vector, leaving the whole argument in suspense. It took years until Saha
(2009), issued an explicit reference to the rotation. His variant of the KS
transformation

$(0,\mbox{\boldmath$x$})=\bar{\mathsf{v}}\,\mathsf{e}_{3}\,\mathsf{v},$ (29)

differs from (23) in two aspects: the conjugation sequence is different, and
the basis quaternion $\mathsf{e}_{3}$ is used instead of $\mathsf{e}_{1}$ (a
rare example of the KS3 convention in celestial mechanics). According to Saha
(2009), equation (29) implies the rotation of the third axis to the $x$
direction, although actually it describes the inverse rotation: $x$ to
$\mbox{\boldmath$e$}_{3}$. Then – up to the signs mismatch – the right-hand
side of (29) coincides with the third column of the matrix
$\mathbf{R}(\mbox{\boldmath$v$})$.

##  3 KS transformation with arbitrary defining vector

###  3.1 Point transformation

The most straightforward generalization of the quaternion formulation proposed
by Deprit et al (1994) is to consider an arbitrary unit quaternion

$\mathsf{c}=\sum_{j=0}^{3}c_{j}\mathsf{e}_{j},$ (30)

and the transformation

$\alpha\mathsf{x}=\mathsf{v}\,\mathsf{c}\,\bar{\mathsf{v}},$ (31)

where a positive real parameter $\alpha$, having the dimension of length, is
introduced as in Deprit et al (1994) to allow the components of $\mathsf{v}$
have the same dimension as $\mathsf{x}$, as well as to facilitate a
comfortable units choice later on. Since the scalar component of (31) is
$x_{0}=c_{0}\,\mathsf{v}\cdot\mathsf{v}$, and $c_{0}$ does not appear in the
vector part of $\mathsf{x}$, we may simply set $c_{0}=0$. Thus
$\mathsf{c}=(0,\mbox{\boldmath$c$})$, and its vector part will be called a
defining vector. By definition, we assume the unit length of the defining
vector $||\mbox{\boldmath$c$}||=1$ throughout the text.

The remaining subsystem of (31) may be set in the matrix-vector form

$\alpha\mbox{\boldmath$x$}=\mathbf{R}(\mathsf{v})\,\mbox{\boldmath$c$},$ (32)

where $\mathbf{R}$ is defined in eq. (28). Alternatively, we can represent
(32) as

$\alpha\mbox{\boldmath$x$}=\left(v_{0}^{2}-\mbox{\boldmath$v$}\cdot\mbox{\boldmath$v$}\right)\,\mbox{\boldmath$c$}+2\left(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$v$}\right)\mbox{\boldmath$v$}+2v_{0}\mbox{\boldmath$v$}\times\mbox{\boldmath$c$},$
(33)

or

$\alpha\mbox{\boldmath$x$}=\left(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$v$}\right)\,\mbox{\boldmath$v$}+\left[\mathsf{v}\wedge(\mathsf{v}\wedge\mathsf{c})\right]^{\natural}.$
(34)

The fibration property, known since Kustaanheimo (1964), may be stated in the
general case as follows: quaternions $\mathsf{v}$ and

$\displaystyle\mathsf{w}$
$\displaystyle=\mathsf{v}\,(\cos{\phi},\sin{\phi}\,\mbox{\boldmath$c$})=$
$\displaystyle=\left(v_{0}\cos{\phi}-(\mbox{\boldmath$v$}\cdot\mbox{\boldmath$c$})\sin{\phi},\,\cos{\phi}\,\mbox{\boldmath$v$}+(v_{0}\mbox{\boldmath$c$}+\mbox{\boldmath$v$}\times\mbox{\boldmath$c$})\sin{\phi}\right),$
(35)

generate the same vector $x$ for all values of angle $\phi$. In other words,
the point $x$ maps onto a fiber consisting of all quaternions $\mathsf{w}$
generated from a given representative $\mathsf{v}$. The proof is elementary,
once we recall equation (18). Then

$\mathsf{w}(0,\mbox{\boldmath$c$})\bar{\mathsf{w}}=\mathsf{v}\left[(\cos{\phi},\sin{\phi}\,\mbox{\boldmath$c$})\,(0,\mbox{\boldmath$c$})\,(\cos{\phi},-\sin{\phi}\,\mbox{\boldmath$c$})\right]\bar{\mathsf{v}}=\mathsf{v}(0,\mbox{\boldmath$c$})\bar{\mathsf{v}},$
(36)

because the part in square brackets describes the rotation of $c$ around
itself. Thus, to a given rotation/scaling matrix $\mathbf{R}(\mathsf{v})$,
exactly two quaternions can be assigned ($\mathsf{v}$ and $-\mathsf{v}$), but
the product of $\mathbf{R}$ and a specified vector $c$ allows more freedom.
This means also that the ‘geometrical interpretation’, stated in Section 2.3,
refers to only one representative of the fiber.

According to the fibration property, inverting the transformation (32) amounts
to picking up some particular $\mathsf{v}$ that serves as the generator of the
fiber. Since $\mathbf{R}$ is homogenous of degree 2, introducing a unit
quaternion $\mathsf{q}=\mathsf{v}/|\mathsf{v}|$ we obtain

$\frac{\alpha\mbox{\boldmath$x$}}{|\mathsf{v}|^{2}}=\mathbf{R}(\mathsf{q})\,\mbox{\boldmath$c$},$
(37)

where $\mathbf{R}(\mathsf{q})\in\mathrm{SO}(3,\mathbb{R})$. Accordingly

$||\mbox{\boldmath$x$}||=r=\frac{\mathsf{v}\cdot\mathsf{v}}{\alpha},$ (38)

since $||\mbox{\boldmath$c$}||=1$ by the assumptions. Thus, the transformation

$\hat{\mbox{\boldmath$x$}}=\mathbf{R}(\mathsf{q})\,\mbox{\boldmath$c$},\qquad\mbox{where}\quad\hat{\mbox{\boldmath$x$}}=\frac{\mbox{\boldmath$x$}}{r},$
(39)

has the meaning of rotation from $c$ to $\hat{\mbox{\boldmath$x$}}$.

Recalling the axis-angle decomposition (25) we can aim at some ‘natural’
choice of $\mathsf{q}$ based upon rotation axis $n$ and angle $\theta$
resulting from elementary vector identities. Thus, the vector part of
$\mathsf{q}$ is

$\mbox{\boldmath$q$}=\sin{\frac{\theta}{2}}\,\mbox{\boldmath$n$}=\frac{\mbox{\boldmath$c$}\times\hat{\mbox{\boldmath$x$}}}{2\,q_{0}},$
(40)

whereas,

$q_{0}=\cos{\frac{\theta}{2}}=\sqrt{\frac{1+\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}}{2}}.$
(41)

Note the singular case $\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}=-1$,
when the actual choice should be $\mathsf{q}=(0,\mbox{\boldmath$n$})$, with an
arbitrary unit vector $n$ orthogonal to $c$ (and thus to
$\hat{\mbox{\boldmath$x$}}$). For a (perturbed) Kepler problem, it may only
happen on a polar orbit, with $c$ placed in the (osculating) orbit plane.
Another problematic situation concerns the collision
$\mbox{\boldmath$x$}=\mbox{\boldmath$0$}$, when the KS quaternion
$\mathsf{v}=\mathsf{0}$ simply cannot be normalized to $q$, and the notion of
rotation is inappropriate.

Thus we first propose an inversion rule

$\mathsf{v}=\sqrt{\frac{\alpha}{2}}\,\left(\sqrt{r+\mbox{\boldmath$c$}\cdot\mbox{\boldmath$x$}},\,\frac{\mbox{\boldmath$c$}\times\mbox{\boldmath$x$}}{\sqrt{r+\mbox{\boldmath$c$}\cdot\mbox{\boldmath$x$}}}\right),\qquad\mbox{for}\quad\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}>-1,$
(42)

and

$\mathsf{v}=\sqrt{\alpha
r}\,\left(0,\mbox{\boldmath$n$}\right),\quad\mbox{\boldmath$n$}\cdot\mbox{\boldmath$x$}=0,\qquad\mbox{for}\quad\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}=-1.$
(43)

It differs from the rules of Stiefel and Scheifele (1971), effectively based
upon the sign of $\mbox{\boldmath$c$}\cdot\mbox{\boldmath$x$}$ (with
$\mbox{\boldmath$c$}=\mbox{\boldmath$e$}_{1}$).

An interesting alternative was proposed by Saha (2009), who used the inversion
rule implying a pure vector $(0,\mbox{\boldmath$v$})$ form. We can obtain it
from (42) through a quaternion product of $\mathsf{v}$ with
$(0,\mbox{\boldmath$c$})$ or its conjugate (both being particular cases of
(35)). The result is indeed considerably simpler:

$\mathsf{v}=\pm\sqrt{\frac{\alpha}{2r\left(1+\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}\right)}}\,\left(0,\,\mbox{\boldmath$x$}+r\mbox{\boldmath$c$}\right),\qquad\mbox{for}\quad\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}>-1,$
(44)

and (43) otherwise. The sign plus or minus can be chosen at will (both
$\mathsf{v}$ and $-\mathsf{v}$ belong to the same fiber). In practice, while
converting a sequence of positions forming an orbit $\mbox{\boldmath$x$}(t)$,
we can swap the signs at the instants, where the motion in $\mathsf{v}$
appears discontinuous. Note that the choice (44) allows plotting the evolution
of KS variables in $\mathbb{R}^{3}$, furnishing spatial trajectory
$\mbox{\boldmath$v$}(t)$. The vector $v$ defined by (44) will be called an SKS
vector (after Saha, Kustaanheimo, and Stiefel). In order to distinguish
general formulae from the ones referring to the SKS vector, we will use the
subscript ‘s’ for the latter.

Any KS quaternion $\mathsf{v}$ with $v_{0}\neq 0$ can be reduced to the SKS
vector by the product

$(0,\mbox{\boldmath$v$}_{\mathrm{s}})=\pm\,\mathsf{v}\,\mathsf{q}_{\mathrm{s}},$
(45)

where the gauge function is a unit quaternion

$\mathsf{q}_{\mathrm{s}}=\frac{\left(\mbox{\boldmath$v$}\cdot\mbox{\boldmath$c$},\,v_{0}\mbox{\boldmath$c$}\right)}{\sqrt{v_{0}^{2}+(\mbox{\boldmath$v$}\cdot\mbox{\boldmath$c$})^{2}}}.$
(46)

This reduction rule can be used if KS coordinates are followed without
reference to the Cartesian position $x$.

###  3.2 Canonical extension

####  3.2.1 KS momenta

Being interested in Hamiltonian formulation of the Kepler problem, we need to
match the KS variables $\mathsf{v}$ with their conjugate momenta $\mathsf{V}$.
This goal can be achieved by a dimension raising Mathieu transformation
$\mbox{\boldmath$X$}\cdot\mathrm{d}\mbox{\boldmath$x$}=\mathsf{V}\cdot\mathrm{d}\mathsf{v}$,
where $X$ are the momenta conjugate to the Cartesian $x$ coordinates. We will
also use a formal quaternion $\mathsf{X}$, postulating $X_{0}=0$. Following
the standard procedure (Kurcheeva, 1977; Deprit et al, 1994), we generalize it
from $\mathsf{e}_{1}$ to a unit quaternion
$\mathsf{c}=(0,\mbox{\boldmath$c$})$ obtaining

$\mathsf{X}=\frac{\mathsf{V}\mathsf{c}\bar{\mathsf{v}}}{2r},$ (47)

wherefrom the vector part is

$X$
$\displaystyle=\frac{(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$v$})\,\mbox{\boldmath$V$}+\left[\,\mathsf{V}\wedge(\mathsf{v}\wedge\mathsf{c})\right]^{\natural}}{2r}=$
$\displaystyle=\frac{1}{2r}\left[(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$v$})\,\mbox{\boldmath$V$}+(v_{0}V_{0}-\mbox{\boldmath$v$}\cdot\mbox{\boldmath$V$})\,\mbox{\boldmath$c$}+(V_{0}\mbox{\boldmath$v$}+v_{0}\mbox{\boldmath$V$})\times\mbox{\boldmath$c$}+(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$V$})\,\mbox{\boldmath$v$}\right].$
(48)

Since the scalar component of (47) should be null, we obtain the constraint

$X_{0}=\frac{\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}}{2r}=0,\qquad\mbox{where}\qquad\mbox{\boldmath$J$}=-v_{0}\mbox{\boldmath$V$}+V_{0}\mbox{\boldmath$v$}+\mbox{\boldmath$v$}\times\mbox{\boldmath$V$}.$
(49)

Remarkably the vector $J$, orthogonal to $c$, is directly related to the
quaternion cross product (20)

$\mathsf{J}=(0,\mbox{\boldmath$J$})=\bar{\mathsf{v}}\wedge\bar{\mathsf{V}},$
(50)

so we can rewrite the condition (49) as

$X_{0}=\frac{\left(\bar{\mathsf{v}}\wedge\bar{\mathsf{V}}\right)\cdot\mathsf{c}}{2r}=0,$
(51)

valid regardless of $c_{0}$. This is the general equivalent of the KS1 phase
space constraint (7) or (10) for the transformation (31). Using
$\mathsf{c}=\mathsf{e}_{1}$, we can recover the formula of Deprit et al
(1994),

$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$e$}_{1}=\mathsf{v}\cdot(\mathsf{V}\mathsf{e}_{1})=0.$
(52)

The presence of constraint (51) allows a unique determination of $V_{0}$ in
terms of the remaining variables. This, combined with a possibility of
reduction to $\mbox{\boldmath$v$}_{\mathrm{s}}$ and $\mathsf{V}_{\mathrm{s}}$
implies that in spite of using eight variables, we follow the dynamics of a
system with effectively three degrees of freedom.

The inverse of the transformation (53) is given as a quaternion product

$\mathsf{V}=\frac{2\,\mathsf{X}\,\mathsf{v}\,\bar{\mathsf{c}}}{\alpha},$ (53)

or, explicitly (setting $X_{0}=0$)

$\displaystyle V_{0}$
$\displaystyle=\frac{2\,\left(\bar{\mathsf{v}}\wedge\mathsf{X}\right)\cdot\mathsf{c}}{\alpha}=\frac{2\,(v_{0}\mbox{\boldmath$c$}+\mbox{\boldmath$v$}\times\mbox{\boldmath$c$})\cdot\mbox{\boldmath$X$}}{\alpha},$
$V$
$\displaystyle=\frac{2}{\alpha}\left[(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$v$})\,\mbox{\boldmath$X$}+(\mbox{\boldmath$v$}\cdot\mbox{\boldmath$X$})\,\mbox{\boldmath$c$}-(\mbox{\boldmath$c$}\cdot\mbox{\boldmath$X$})\,\mbox{\boldmath$v$}+v_{0}(\mbox{\boldmath$c$}\times\mbox{\boldmath$X$})\right].$
(54)

Notably, the definition of the new momenta is given in a mixed form, involving
old momenta $X$ and new coordinates $v$. It means, that for a given set of
values $\mbox{\boldmath$x$},\mbox{\boldmath$X$}$, the phase space fiber
contains not only the family of coordinates $\mathsf{v}$ implied by (35), but
also the family of momenta $\mathsf{V}$ – one quaternion for each member of
(35). Thus, if we try to provide the explicit form of
$\mathsf{V}(\mbox{\boldmath$x$},\mbox{\boldmath$X$})$, we may choose some
particular representative of the fiber. Let us comfortably choose the SKS
vector (44), because then we can simplify expressions, remaining on the ground
of usual vector calculus.

In the absence of $v_{0}$, the scalar part $V_{0}$ simplifies to
$V_{0\mathrm{s}}=-2(\mbox{\boldmath$v$}_{\mathrm{s}}\times\mbox{\boldmath$X$})\cdot\mbox{\boldmath$c$}/\alpha$.
Substituting $\mbox{\boldmath$v$}_{\mathrm{s}}$ from (44), with the plus sign
selected, into the first of equations (54), we find

$V_{0\mathrm{s}}=-\sqrt{\frac{2}{\alpha\,r\left(1+\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}\right)}}\,(\mbox{\boldmath$x$}\times\mbox{\boldmath$X$})\cdot\mbox{\boldmath$c$},$
(55)

whereas for $\hat{\mbox{\boldmath$x$}}=-\mbox{\boldmath$c$}$, the scalar part
of $\mathsf{V}$ is $V_{0\mathrm{s}}=0$. In both cases the conclusion is the
same: whenever the SKS vector is taken for coordinates, the scalar part of the
KS momenta quaternion is a product of a coordinates dependent factor and the
projection of angular momentum on the defining vector $c$.

The vector part $V$ is also linked with familiar quantities when the same SKS
vector is used, leading to

$\mbox{\boldmath$V$}_{\mathrm{s}}=\sqrt{\frac{2}{\alpha\,r\left(1+\mbox{\boldmath$c$}\cdot\hat{\mbox{\boldmath$x$}}\right)}}\,\left(r\mbox{\boldmath$X$}+(\mbox{\boldmath$x$}\cdot\mbox{\boldmath$X$})\,\mbox{\boldmath$c$}+(\mbox{\boldmath$x$}\times\mbox{\boldmath$X$})\times\mbox{\boldmath$c$}\right),$
(56)

where the Cartesian momentum, radial velocity and angular momentum appear.

Finding $\mathsf{V}_{\mathrm{s}}$ is possible without the knowledge of $x$ and
$X$. Let us multiply both sides of eq. (53) by a quaternion product
$\mathsf{c}\mathsf{q}_{\mathrm{s}}\bar{\mathsf{c}}$. Then, the equation
becomes

$\mathsf{V}\mathsf{c}\mathsf{q}_{\mathrm{s}}\bar{\mathsf{c}}=\frac{2\,\mathsf{X}\,(\mathsf{v}\,\mathsf{q}_{\mathrm{s}})\,\bar{\mathsf{c}}}{\alpha},$
(57)

and we see that its left hand side should represent $\mathsf{V}_{\mathrm{s}}$.
But one may easily verify that
$\mathsf{c}\mathsf{q}_{\mathrm{s}}\bar{\mathsf{c}}=\mathsf{q}_{\mathrm{s}}$,
so if momenta $\mathsf{V}$ are determined by $X$ with an arbitrary quaternion
$\mathsf{v}$, then the momenta $\mathsf{V}_{\mathrm{s}}$ determined by the
same $X$ and the equivalent SKS vector $\mbox{\boldmath$v$}_{\mathrm{s}}$ are

$\mathsf{V}_{\mathrm{s}}=\pm\mathsf{V}\,\mathsf{q}_{\mathrm{s}},$ (58)

where the sign choice should be the same as in (45).

##  4 KS Hamiltonian and its invariants

###  4.1 Perturbed Hamiltonian and equations of motion

As long as the defining vector is constant, the canonical KS transformation is
time independent, so the Hamiltonian function transforms without a remainder.
Thus, the first step is to obtain
$\mathcal{H}^{\star}(\mathsf{v},\mathsf{V},t)=\mathcal{H}(\mbox{\boldmath$x$},\mbox{\boldmath$X$},t)$,
where $\mathcal{H}$ is a perturbed two body problem Hamiltonian

$\mathcal{H}(\mbox{\boldmath$x$},\mbox{\boldmath$X$},t)=\mathcal{H}_{0}(\mbox{\boldmath$x$},\mbox{\boldmath$X$})+\mathcal{R}(\mbox{\boldmath$x$},\mbox{\boldmath$X$},t),$
(59)

with the Keplerian part

$\mathcal{H}_{0}=\frac{\mbox{\boldmath$X$}\cdot\mbox{\boldmath$X$}}{2}-\frac{\mu}{r},$
(60)

depending on the gravitational parameter $\mu$, and we make no assumptions
about the order of magnitude for the perturbation $\mathcal{R}$.

The distance $r$ may be considered a known function of $v$ thanks to (38), so
we retain this symbol in $\mathcal{H}_{0}$. The square of
$|\mbox{\boldmath$X$}|$ is easily found from
$\mathsf{X}\bar{\mathsf{X}}=\mbox{\boldmath$X$}\cdot\mbox{\boldmath$X$}+X_{0}^{2}$.
Substituting eq. (47), using the rule (18) and defining $X_{0}$ by (49), we
find666The equivalent equation (22) of Deprit et al (1994) is incomplete by
the omission of $X_{0}^{2}$.

$\mbox{\boldmath$X$}\cdot\mbox{\boldmath$X$}=\frac{\alpha}{4r}\mathsf{V}\cdot\mathsf{V}-\frac{(\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$})^{2}}{4r^{2}}.$
(61)

Thus the transformed Keplerian Hamiltonian is

$\mathcal{H}^{\star}_{0}(\mathsf{v},\mathsf{V})=\frac{\alpha}{8r}\mathsf{V}\cdot\mathsf{V}-\frac{\mu}{r}-\frac{(\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$})^{2}}{8r^{2}}.$
(62)

Dropping the last, zero valued term in (62) is allowed, but it should not be
done without reflection. It is to be remembered that while $x_{0}=0$ is the
property of the point transformation (31) itself, $X_{0}=0$ is only
postulated. Two separate questions should be addressed. First: is the value of
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$ conserved during the motion?
Here the answer is conditionally positive: it is not so for an arbitrarily
invented Hamiltonian function of $\mathsf{v}$ and $\mathsf{V}$. But if the
Hamiltonian is a transformed function
$\mathcal{H}(\mbox{\boldmath$x$},\mbox{\boldmath$X$},t)$, then its KS image
$\mathcal{H}^{\star}(\mathsf{v},\mathsf{V},t)$ conserves
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$, because all the Poisson’s
brackets
$\left\\{\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$},x_{j}\right\\}=\left\\{\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$},X_{j}\right\\}=0$,
for $0\leq j\leq 3$, similarly to the argument of Deprit et al (1994). The
second question, less often considered, is: does the presence of
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$ influence the solution? Worth
asking, because a zero valued function may still have nonzero derivatives. For
a while the answer is obvious, because $\mathcal{H}_{0}^{\star}$ contains the
square $(\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$})^{2}$, so its gradient
will vanish. But the problem may reappear in the context of variational
equations or when the rotating reference frame will be considered.

The Hamiltonian (60) remains singular at $r=0$, so we proceed with the Sundman
transformation (5)

$\frac{\mathrm{d}\tau}{\mathrm{d}t}=\frac{\beta}{r},$ (63)

with an arbitrary parameter $\beta$. In the canonical framework, switching
from physical time $t$ to the Sundman time $\tau$ requires the transition to
the extended phase space, appending to the KS coordinates and momenta a
conjugate pair $v^{\ast}$ and $V^{\ast}$. The former is an imitator of the
physical time (up to an additive constant); the latter serves to fix a zero
energy manifold for the motion, being a doppelganger of the Hamiltonian
$\mathcal{H}^{\ast}$. Thus, on the manifold $\mathcal{H}^{\star}+V^{\ast}=0$,
we can divide the extended Hamiltonian by the right hand side of (63),
obtaining the Hamiltonian function

$\mathcal{K}(\mathsf{v},\mathsf{V},v^{\ast},V^{\ast})=\frac{\alpha}{8\beta}\mathsf{V}\cdot\mathsf{V}-\frac{\mu}{\beta}+\frac{r}{\beta}\mathcal{R}^{\star}(\mathsf{v},\mathsf{V},v^{\ast})+\frac{V^{\ast}r}{\beta}=0,$
(64)

independent on the new time variable $\tau$. What remains, is a judicious
choice of $\alpha$ and $\beta$. Assuming

$\beta=\frac{\alpha}{4},\qquad\mbox{hence}\quad\frac{\mathrm{d}\tau}{\mathrm{d}t}=\frac{\alpha}{4r},$
(65)

we secure $\mathsf{v}^{\prime}=\mathsf{V}$. If then $\alpha$ is equal to major
axis (i.e. $\alpha=2a$) of the elliptic orbit, $\tau$ will run on average at
half rate of $t$ for the Kepler problem, in accord with the angle doubling
property of the KS transformation.

Although the constant term has no influence on equations of motion, we retain
it for its role in fixing the $\mathcal{K}=0$ manifold. And so we finally set
up the KS Hamiltonian

$\displaystyle\mathcal{K}(\mathsf{v},\mathsf{V},v^{\ast},V^{\ast})$
$\displaystyle=\mathcal{K}_{0}(\mathsf{v},\mathsf{V},V^{\ast})+\mathcal{P}(\mathsf{v},\mathsf{V},v^{\ast}),$
(66) $\displaystyle\mathcal{K}_{0}(\mathsf{v},\mathsf{V},V^{\ast})$
$\displaystyle=\frac{1}{2}\,\mathsf{V}\cdot\mathsf{V}+\frac{4V^{\ast}}{\alpha^{2}}\mathsf{v}\cdot\mathsf{v}-\frac{4\mu}{\alpha},$
(67) $\displaystyle\mathcal{P}(\mathsf{v},\mathsf{V},v^{\ast})$
$\displaystyle=\frac{4r}{\alpha}\mathcal{R}^{\star}(\mathsf{v},\mathsf{V},v^{\ast}).$
(68)

The value of $V^{\ast}$ should secure $\mathcal{K}=0$.

Equations of motion resulting from (66) are those of a perturbed harmonic
oscillator with frequency

$\omega_{0}=\frac{2\sqrt{2V^{\ast}}}{\alpha},$ (69)

namely

$\displaystyle\mathsf{v}^{\prime}$
$\displaystyle=\left\\{\mathsf{v},\mathcal{K}\right\\}=\mathsf{V}+\sum_{j=0}^{3}\frac{\partial\mathcal{P}}{\partial
V_{j}}\mathsf{e}_{j},$ (70) $\displaystyle\mathsf{V}^{\prime}$
$\displaystyle=\left\\{\mathsf{V},\mathcal{K}\right\\}=-\omega_{0}^{2}\,\mathsf{v}-\sum_{j=0}^{3}\frac{\partial\mathcal{P}}{\partial
v_{j}}\mathsf{e}_{j},$ (71) $\displaystyle(v^{\ast})^{\prime}$
$\displaystyle=\left\\{v^{\ast},\mathcal{K}\right\\}=\frac{4r}{\alpha},$ (72)
$\displaystyle(V^{\ast})^{\prime}$
$\displaystyle=\left\\{V^{\ast},\mathcal{K}\right\\}=\frac{\partial\mathcal{P}}{\partial
v^{\ast}},$ (73)

where the prime marks the derivative with respect to Sundman time $\tau$. We
can note, that using an arbitrary defining vector $c$ has no influence on the
unperturbed problem in KS variables. Whatever change results from a particular
choice of $c$, may be revealed only by the form taken by $\mathcal{P}$ in a
specific problem.

###  4.2 Invariants

Including the perturbation $\mathcal{P}$, we can only mention two invariants:
if $\mathcal{P}$ does not depend explicitly on $v^{\ast}$ (hence on time $t$),
the momentum $V^{\ast}=-E=\mathrm{const}$, where $E$ is the total energy;
regardless of $\mathcal{P}$, the scalar product
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}=0$ is also invariant, as already
mentioned. Thus, let us focus on the first integrals of the unperturbed system
$\mathcal{K}_{0}$, i.e. the two body problem.

There are two points of view for the unperturbed system. We can see it as a
four dimensional isotropic harmonic oscillator with its own first integrals.
But we can also see it as a transformed Kepler problem with the well known
first integrals, potentially expressible in terms of the oscillator constants.

####  4.2.1 Oscillator

In the absence of perturbation, we can consider the Hamiltonian
$\mathcal{K}_{0}$ as a separable system of four independent oscillators

$\mathcal{K}_{0}=-\frac{4\mu}{\alpha}+\sum_{j=0}^{3}\mathcal{N}_{j}=0,\qquad\mbox{where}\quad\mathcal{N}_{j}=\frac{V_{j}^{2}}{2}+\frac{\omega_{0}^{2}v_{j}^{2}}{2}=E_{j}.$
(74)

Each Hamiltonian $\mathcal{N}_{j}$ is a first integral in involution with the
rest $\left\\{\mathcal{N}_{i},\mathcal{N}_{j}\right\\}=0$, but only three of
them are independent, since their sum is fixed by (74). Yet, the system is
superintegrable and more integrals can be found. First, we can introduce a
four-dimensional variant of the Fradkin tensor $\mathbf{F}$ (Fradkin, 1967)

$F_{ij}=\frac{V_{i}V_{j}}{\omega_{0}}+\omega_{0}v_{i}v_{j}.$ (75)

The symmetric matrix $\mathbf{F}$ contains 10 different first integrals (not
all independent), including four diagonal terms $F_{ii}=2E_{i}/\omega_{0}$.

Another set of integrals constitutes an antisymmetric angular momentum matrix
$\mathbf{L}$

$L_{ij}=v_{i}V_{j}-v_{j}V_{i},$ (76)

with 6 distinct elements (again, not all independent).

Considering $F_{ij}$ and $L_{ij}$ as the generators of Hamiltonian equations,
we signal their different, somewhat complementary roles. Equations

$\mathsf{v}^{\prime}=\left\\{\mathsf{v},F_{ij}\right\\}=\frac{V_{j}\mathsf{e}_{i}+V_{i}\mathsf{e}_{j}}{\omega_{0}},\qquad\mathsf{V}^{\prime}=\left\\{\mathsf{V},F_{ij}\right\\}=-\omega_{0}\left(v_{j}\mathsf{e}_{i}+v_{i}\mathsf{e}_{j}\right),$
(77)

define a phase plane rotation: either in a phase plane
$(v_{i},V_{i}/\omega_{0})$ (diagonal terms $F_{ii}$), or in two phase planes
$(v_{i},V_{i}/\omega_{0})$ and $(v_{j},V_{j}/\omega_{0})$. The angular
momentum terms $L_{ij}$ lead to equations

$\mathsf{v}^{\prime}=\left\\{\mathsf{v},L_{ij}\right\\}=-v_{j}\mathsf{e}_{i}+v_{i}\mathsf{e}_{j},\qquad\mathsf{V}^{\prime}=\left\\{\mathsf{V},L_{ij}\right\\}=-V_{j}\mathsf{e}_{i}+V_{i}\mathsf{e}_{j},$
(78)

that generate rotation on a coordinate plane $(v_{i},v_{j})$ and on a momentum
plane $(V_{i},V_{j})$.

Two important cross products that appeared in the KS transformation are
expressible in terms of $L_{ij}$ and thus are first integrals of unperturbed
motion:

$\displaystyle\left(\mathsf{v}\wedge\mathsf{V}\right)^{\natural}$
$\displaystyle=\left(L_{01}+L_{23}\right)\mbox{\boldmath$e$}_{1}+\left(L_{02}+L_{31}\right)\mbox{\boldmath$e$}_{2}+\left(L_{03}+L_{12}\right)\mbox{\boldmath$e$}_{3},$
(79)
$\displaystyle\left(\bar{\mathsf{v}}\wedge\bar{\mathsf{V}}\right)^{\natural}$
$\displaystyle=\left(L_{10}+L_{23}\right)\mbox{\boldmath$e$}_{1}+\left(L_{20}+L_{31}\right)\mbox{\boldmath$e$}_{2}+\left(L_{30}+L_{12}\right)\mbox{\boldmath$e$}_{3}.$
(80)

Accordingly, the condition (51) can be seen as a constraint on the angular
momentum of the oscillator.

####  4.2.2 Kepler problem

The energy integral of the Kepler problem has been already discussed in
Section 4, so let us pass to the two vector-valued first integrals: angular
momentum and Laplace (Runge-Lenz) vector.

Expressing the angular momentum

$\mbox{\boldmath$G$}=\mbox{\boldmath$x$}\times\mbox{\boldmath$X$},$ (81)

in terms of the KS variables is a formidable task if a brute force attack is
attempted by substitution of (33) and the second line of (48) into (81). It is
much better to start from plugging in a cross product
$\mathsf{x}\wedge\mathsf{X}$, so that

$\mbox{\boldmath$x$}\times\mbox{\boldmath$X$}=\left[\,\mathsf{x}\wedge\mathsf{X}\right]^{\natural}+X_{0}\mbox{\boldmath$x$},$
(82)

with $x_{0}=0$. Substituting the quaternion product forms (31) for
$\mathsf{x}$ and (47) for $\mathsf{X}$, one can resort to the factor exchange
rule (21)

$\mathsf{x}\wedge\mathsf{X}=\frac{\left(\mathsf{v}\mathsf{c}\bar{\mathsf{v}}\right)\wedge\left(\mathsf{V}\mathsf{c}\bar{\mathsf{v}}\right)}{2\alpha
r}=\frac{\left(\mathsf{v}\mathsf{c}\right)\wedge\left(\mathsf{V}\mathsf{c}\bar{\mathsf{v}}\mathsf{v}\right)}{2\alpha
r}=\frac{\mathsf{v}\wedge\left(\mathsf{V}\mathsf{c}\bar{\mathsf{c}}\right)}{2}=\frac{\mathsf{v}\wedge\mathsf{V}}{2},$
(83)

where we used $\bar{\mathsf{v}}\mathsf{v}=\alpha r$, and
$\mathsf{c}\bar{\mathsf{c}}=1$. Thus we obtain

$\mbox{\boldmath$G$}=\mbox{\boldmath$x$}\times\mbox{\boldmath$X$}=\frac{\left[\,\mathsf{v}\wedge\mathsf{V}\right]^{\natural}}{2}+X_{0}\mbox{\boldmath$x$}=\frac{1}{2}\left[\,\mathsf{v}\wedge\mathsf{V}+\frac{\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}}{\alpha
r}\,\mathsf{v}\mathsf{c}\bar{\mathsf{v}}\right]^{\natural}.$ (84)

An analogous expression was obtained by Deprit et al (1994) with a statement:
‘Proof. - By straightforward calculation using Symbol Processor’. We have
decided to provide the proof in full length (or rather shortness), as a good
example of the situation where quaternion formalism beats standard vector
calculus.

If the expression (84) is used for the sole purpose of computing the value of
$G$, then $X_{0}=0$ can be safely set. But if $G$ is to appear in the
perturbing Hamiltonian $\mathcal{P}$, one should not forget that derivatives
of $X_{0}$ do not vanish in general.

Linking the Keplerian $G$ with the oscillator’s angular momentum tensor
$\mathbf{L}$ is straightforward: inserting (79) into (84) results in

$\mbox{\boldmath$G$}=\frac{L_{01}+L_{23}}{2}\mbox{\boldmath$e$}_{1}+\frac{L_{02}+L_{31}}{2}\mbox{\boldmath$e$}_{2}+\frac{L_{03}+L_{12}}{2}\mbox{\boldmath$e$}_{3}+X_{0}\mbox{\boldmath$x$}.$
(85)

Notably, the defining vector $c$ has no direct effect on the direction of the
angular momentum. Its action is only indirect, through the constraint
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}=0$.

The Laplace vector $e$ is primarily given with a cross product of momentum and
angular momentum, but in order to express it in terms of KS variables, a
version resulting form the ‘BAC-CAB’ identity is more convenient

$\mu\mbox{\boldmath$e$}=\mbox{\boldmath$X$}\times\mbox{\boldmath$G$}-\mu\hat{\mbox{\boldmath$x$}}=\left(\mbox{\boldmath$X$}\cdot\mbox{\boldmath$X$}-\mu
r^{-1}\right)\,\mbox{\boldmath$x$}-(\mbox{\boldmath$x$}\cdot\mbox{\boldmath$X$})\,\mbox{\boldmath$X$}.$
(86)

Most of the building blocks are ready in equations (33), (34), (48), and (61).
The remaining product is elementarily found applying the identity (19)

$\mbox{\boldmath$x$}\cdot\mbox{\boldmath$X$}=\mathsf{x}\cdot\mathsf{X}=\frac{(\mathsf{v}\mathsf{c}\bar{\mathsf{v}})\cdot(\mathsf{V}\mathsf{c}\bar{\mathsf{v}})}{2\alpha
r}=\frac{\bar{\mathsf{V}}\cdot\left(\mathsf{c}\bar{\mathsf{v}}\overline{(\mathsf{v}\mathsf{c}\bar{\mathsf{v}})}\right)}{2\alpha
r}=\frac{\bar{\mathsf{V}}\cdot\left(\mathsf{c}\bar{\mathsf{v}}\mathsf{v}\bar{\mathsf{c}}\bar{\mathsf{v}}\right)}{2\alpha
r}=\frac{\bar{\mathsf{V}}\cdot\bar{\mathsf{v}}}{2}=\frac{\mathsf{v}\cdot\mathsf{V}}{2}.$
(87)

And so the Laplace vector can be expressed as a quaternion product

$\mu\mbox{\boldmath$e$}=\frac{1}{2r}\left[\left(\left(\frac{\mathsf{V}\cdot\mathsf{V}}{2}-\frac{2\mu}{\alpha}+\frac{(\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$})^{2}}{2\alpha
r}\right)\mathsf{v}-\frac{\mathsf{v}\cdot\mathsf{V}}{2}\mathsf{V}\right)\mathsf{c}\bar{\mathsf{v}}\right]^{\natural}.$
(88)

The formula does not look friendly, but one should expect it to be expressible
in terms of Fradkin integrals. Indeed, after rather tedious manipulations, we
have found that

$\mu\mbox{\boldmath$e$}=-\frac{\alpha\omega_{0}}{4}\,\mathbf{E}\mbox{\boldmath$c$}-X_{0}\mbox{\boldmath$G$}+\frac{\alpha\mathcal{K}_{0}}{4r}\mbox{\boldmath$x$},$
(89)

where the second term to the right vanishes due to the KS constraint
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}=0$, and the third is null on the
Keplerian manifold $\mathcal{K}_{0}=0$. It is not by chance, that the form of
matrix $\mathbf{E}$

$\mathbf{E}=\left(\begin{array}[]{ccc}E_{11}&\quad F_{12}-F_{03}&\quad
F_{13}+F_{02}\\\ F_{12}+F_{03}&E_{22}&\quad F_{23}-F_{01}\\\
F_{13}-F_{02}&\quad F_{23}+F_{01}&E_{33}\\\ \end{array}\right),$ (90)

with diagonal terms

$\displaystyle E_{11}$ $\displaystyle=$
$\displaystyle\frac{F_{00}+F_{11}-F_{22}-F_{33}}{2}=\frac{\mathcal{K}_{0}}{\omega_{0}}+\frac{4\mu}{\alpha\omega_{0}}-F_{22}-F_{33},$
$\displaystyle E_{22}$ $\displaystyle=$
$\displaystyle\frac{F_{00}-F_{11}+F_{22}-F_{33}}{2}=\frac{\mathcal{K}_{0}}{\omega_{0}}+\frac{4\mu}{\alpha\omega_{0}}-F_{11}-F_{33},$
(91) $\displaystyle E_{33}$ $\displaystyle=$
$\displaystyle\frac{F_{00}-F_{11}-F_{22}+F_{33}}{2}=\frac{\mathcal{K}_{0}}{\omega_{0}}+\frac{4\mu}{\alpha\omega_{0}}-F_{11}-F_{22},$

mimics the matrix $\mathbf{R}(\mbox{\boldmath$v$})$ defined by eq. (28),
present in the KS transformation formula (32). In contrast to the angular
momentum $G$, the defining vector $c$ is explicitly present in the definition
of $e$.

###  4.3 Dynamical role of the invariant
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$

In Sections 4.1 and 4.2.2, some warnings have been issued concerning the
presence of the invariant $\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}=0$,
which should not be dropped blindly in some expressions. Let us now inspect
its influence, by considering a Hamiltonian
$\mathcal{M}=\Psi\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$, where $\Psi$ is
an arbitrary function of KS variables. Canonical equations of motion generated
by $\mathcal{M}$, can be cast into a quaternion product form

$\mathsf{v}^{\prime}=\left\\{\mathsf{v},\mathcal{M}\right\\}=\Psi\,\mathsf{v}\bar{\mathsf{c}},\qquad\mathsf{V}^{\prime}=\left\\{\mathsf{V},\mathcal{M}\right\\}=\Psi\,\mathsf{V}\bar{\mathsf{c}}.$
(92)

The solution of this system is a quaternion product

$\mathsf{v}=\mathsf{u}\left(\cos{\phi},\sin{\phi}\,\mbox{\boldmath$c$}\right),\qquad\mathsf{V}=\mathsf{U}\left(\cos{\phi},\sin{\phi}\,\mbox{\boldmath$c$}\right),$
(93)

with arbitrary constants $\mathsf{u},\mathsf{U}$, provided $\phi$ is a
function of time (possibly implicit) satisfying $\phi^{\prime}=-\Psi$. But,
recalling the fiber definition (35), we see that the resulting evolution of KS
variables (93) happens on a fiber referring to constant values of the
Cartesian variables.

Considering any Hamiltonian $\mathcal{K}+\mathcal{M}$, where $\mathcal{K}$ is
a KS transform of some function of the Cartesian variables, we recall that
both terms commute ($\left\\{\mathcal{K},\mathcal{M}\right\\}=0$), so the the
solution will be a direct composition of flows resulting from both
Hamiltonians separately. As it follows, adding or retracting a
$\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}$ term to a Hamiltonian function
can modify a trajectory in the KS phase space, but it has no influence on the
resulting motion in the Cartesian variables
$\mbox{\boldmath$x$},\mbox{\boldmath$X$}$. We encourage an interested reader
to consult a related work of Roa et al (2016) in the framework of the KS1 set.

##  5 Kepler problem in rotating reference frame

###  5.1 General equations of motion

The Kepler problem in a uniformly rotating reference frame is a necessary
building block for a number of dynamical problems handled by analytical
perturbation techniques or symplectic integrators with partitioned
Hamiltonian. Recently, the problem has been solved by Langner and Breiter
(2015), where the account of earlier works can be found as well. Since Langner
and Breiter (2015) solved the problem in the frame rotating around
$\mbox{\boldmath$e$}_{3}$ using the KS1 set (based upon the defining vector
$\mbox{\boldmath$c$}=\mbox{\boldmath$e$}_{1}$), it may be interesting to
confront the solution with a new one, assuming an arbitrary direction of
rotation axis and benefiting from the freedom in the defining vector choice.
Intuitively, selecting $c$ directed along the rotation axis seems most
appropriate, so we assume the angular velocity vector of the reference frame
to be $\Omega\mbox{\boldmath$c$}$ from the onset.

The KS transformation, as described in Section 3, will be applied to the
coordinates $x$ and momenta $X$ of the rotating frame. If at the epoch $t=0$
the rotating frame and the fixed frame axes coincide, then the transformation
linking $\mbox{\boldmath$x$},\mbox{\boldmath$X$}$ with the fixed frame
coordinates $\mbox{\boldmath$x$}_{\mathrm{f}}$ and momenta
$\mbox{\boldmath$X$}_{\mathrm{f}}$

$\mbox{\boldmath$x$}=\mathbf{R}(\mathsf{q})\mbox{\boldmath$x$}_{\mathrm{f}},\qquad\mbox{\boldmath$X$}=\mathbf{R}(\mathsf{q})\mbox{\boldmath$X$}_{\mathrm{f}},$
(94)

involves rotation matrix from eq. (28) with the rotation quaternion

$\mathsf{q}=\left(\cos{\left(\frac{\Omega
t}{2}\right)},-\sin{\left(\frac{\Omega
t}{2}\right)}\,\mbox{\boldmath$c$}\right).$ (95)

The transformation is canonical and, being time-dependent, it creates the
remainder $-\Omega\mbox{\boldmath$G$}\cdot\mbox{\boldmath$c$}$ supplementing
the transformed Hamiltonian. Accordingly, the Hamiltonian function to be
considered is a sum of $\mathcal{K}_{0}$ from eq. (66) and of

$\mathcal{P}=-\frac{4r\Omega}{\alpha}\,\mbox{\boldmath$G$}\cdot\mbox{\boldmath$c$}.$
(96)

###  5.2 Simplification

Recalling the conclusion of Section 4.3, and benefiting from the choice of
$c$, we can modify $\mathcal{P}$ and use

$\mathcal{P}_{\mathrm{m}}=\mathcal{P}-\frac{2\Omega}{\alpha}\left(\mbox{\boldmath$x$}\cdot\mbox{\boldmath$c$}+r\right)\,\mbox{\boldmath$J$}\cdot\mbox{\boldmath$c$}.$
(97)

According to equations (84), (79) and (50), the modified term is simply

$\mathcal{P}_{\mathrm{m}}=-\frac{4r}{\alpha}\Omega H,\qquad
H=\left(\mbox{\boldmath$v$}\times\mbox{\boldmath$V$}\right)\cdot\mbox{\boldmath$c$},$
(98)

so it contains only the vector parts of the KS quaternions, save for $v_{0}$
present in $r$.

Deriving canonical equations of motion from

$\mathcal{K}=\frac{\mathsf{V}\cdot\mathsf{V}}{2}+\frac{4V^{\ast}\,\mathsf{v}\cdot\mathsf{v}}{\alpha^{2}}-\frac{4\mu}{\alpha}-\frac{4\,\Omega\,\mathsf{v}\cdot\mathsf{v}}{\alpha^{2}}\,\left(\mbox{\boldmath$v$}\times\mbox{\boldmath$V$}\right)\cdot\mbox{\boldmath$c$},$
(99)

we observe that they neatly split into a scalar part

$v_{0}^{\prime}=V_{0},\qquad V_{0}^{\prime}=-w^{2}v_{0},$ (100)

and the vector part

$\mbox{\boldmath$v$}^{\prime}=\mbox{\boldmath$V$}-\frac{4r}{\alpha}\left(\Omega\mbox{\boldmath$c$}\times\mbox{\boldmath$v$}\right),\qquad\mbox{\boldmath$V$}^{\prime}=-w^{2}\mbox{\boldmath$v$}-\frac{4r}{\alpha}\left(\Omega\mbox{\boldmath$c$}\times\mbox{\boldmath$V$}\right),$
(101)

where the frequency

$w=\frac{2\sqrt{2(V^{\ast}-\Omega H)}}{\alpha},$ (102)

is a constant of motion, because
$H^{\prime}=\left\\{H,\mathcal{K}\right\\}=0$, and the last two equations of
motion are

$(v^{\ast})^{\prime}=\frac{4r}{\alpha},\qquad(V^{\ast})^{\prime}=0.$ (103)

Initial conditions for this system at $\tau=0$ will be

$\mathsf{v}(0)=\mathsf{u},\qquad\mathsf{V}(0)=\mathsf{U},\qquad
v^{\ast}(0)=t=0.$ (104)

Thus, the situation is much more comfortable than in Langner and Breiter
(2015) and all earlier works. Equations (100), describing a simple, one-
dimensional harmonic oscillator, are easily solved, rendering

$v_{0}=\cos{(w\tau)}\,u_{0}+\sin{(w\tau)}\,\frac{U_{0}}{w},\qquad
V_{0}=-w\sin{(w\tau)}\,u_{0}+\cos{(w\tau)}\,U_{0},$ (105)

similarly to the fixed frame case.

Looking at the equations (101), we recognize two parts referring to harmonic
oscillator dynamics and to the kinematics of rotation. Introducing the cross
product matrix

$\mathbf{C}=\Omega\,\left(\begin{array}[]{ccc}0&-c_{3}&c_{2}\\\
c_{3}&0&-c_{1}\\\ -c_{2}&c_{1}&0\\\ \end{array}\right),$ (106)

we rewrite (101) in the vector-matrix form

$\frac{\mathrm{d}\mbox{\boldmath$v$}}{\mathrm{d}\tau}=\mbox{\boldmath$V$}-\frac{4r}{\alpha}\mathbf{C}\mbox{\boldmath$v$},\qquad\frac{\mathrm{d}\mbox{\boldmath$V$}}{\mathrm{d}\tau}=-w^{2}\mbox{\boldmath$v$}-\frac{4r}{\alpha}\mathbf{C}\mbox{\boldmath$V$},$
(107)

which suggest to postulate the solution

$\mbox{\boldmath$v$}=\mathbf{A}\,\left(b_{1}\mbox{\boldmath$u$}+b_{2}\mbox{\boldmath$U$}\right),\qquad\mbox{\boldmath$V$}=\mathbf{A}\,\left(b_{3}\mbox{\boldmath$u$}+b_{4}\mbox{\boldmath$U$}\right),$
(108)

involving a common matrix $\mathbf{A}$ and four scalars $b_{j}$, with the
initial conditions $\mathbf{A}=\mathbf{I}$, $b_{2}=b_{3}=0$, and
$b_{1}=b_{4}=1$ at $\tau=t=v^{\ast}=0$. Substitution into (107) leads to

$\displaystyle\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}\tau}\left(b_{1}\mbox{\boldmath$u$}+b_{2}\mbox{\boldmath$U$}\right)+\mathbf{A}\left(\frac{\mathrm{d}b_{1}}{\mathrm{d}\tau}\mbox{\boldmath$u$}+\frac{\mathrm{d}b_{2}}{\mathrm{d}\tau}\mbox{\boldmath$U$}\right)$
$\displaystyle=$
$\displaystyle\mathbf{A}\,\left(b_{3}\mbox{\boldmath$u$}+b_{4}\mbox{\boldmath$U$}\right)$
(109)
$\displaystyle-\frac{4r}{\alpha}\mathbf{C}\mathbf{A}\,\left(b_{1}\mbox{\boldmath$u$}+b_{2}\mbox{\boldmath$U$}\right),$
$\displaystyle\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}\tau}\left(b_{3}\mbox{\boldmath$u$}+b_{4}\mbox{\boldmath$U$}\right)+\mathbf{A}\left(\frac{\mathrm{d}b_{3}}{\mathrm{d}\tau}\mbox{\boldmath$u$}+\frac{\mathrm{d}b_{4}}{\mathrm{d}\tau}\mbox{\boldmath$U$}\right)$
$\displaystyle=$
$\displaystyle-w^{2}\mathbf{A}\,\left(b_{1}\mbox{\boldmath$u$}+b_{2}\mbox{\boldmath$U$}\right)$
(110)
$\displaystyle-\frac{4r}{\alpha}\mathbf{C}\mathbf{A}\,\left(b_{3}\mbox{\boldmath$u$}+b_{4}\mbox{\boldmath$U$}\right).$

Collecting the term preceded by $\mathbf{A}$, we obtain the system

$\displaystyle\left(\frac{\mathrm{d}b_{1}}{\mathrm{d}\tau}-b_{3}\right)\mbox{\boldmath$u$}+\left(\frac{\mathrm{d}b_{2}}{\mathrm{d}\tau}-b_{4}\right)\mbox{\boldmath$U$}$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$0$},$
$\displaystyle\left(\frac{\mathrm{d}b_{3}}{\mathrm{d}\tau}+w^{2}b_{1}\right)\mbox{\boldmath$u$}+\left(\frac{\mathrm{d}b_{4}}{\mathrm{d}\tau}+w^{2}b_{2}\right)\mbox{\boldmath$U$}$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$0$},$ (111)

with an obvious solution

$b_{1}=b_{4}=\cos{(w\tau)},\qquad b_{2}=\frac{\sin{(w\tau)}}{w},\qquad
b_{3}=-w\sin{(w\tau)},$ (112)

actually known from (105). In the remaining part of (109) and (110) we change
the independent variable using (103), and letting $v^{\ast}=t$ for brevity, we
find

$\displaystyle\frac{4r}{\alpha}\left(\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}+\mathbf{C}\mathbf{A}\right)\left(b_{1}\mbox{\boldmath$u$}+b_{2}\mbox{\boldmath$U$}\right)$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$0$},$ (113)
$\displaystyle\frac{4r}{\alpha}\left(\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}+\mathbf{C}\mathbf{A}\right)\left(b_{3}\mbox{\boldmath$u$}+b_{4}\mbox{\boldmath$U$}\right)$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$0$},$ (114)

solved by the orthogonal matrix $\mathbf{A}$, which represents rotation around
$c$ by an angle $(-\Omega t)$. Thus, the final solution consists of the scalar
equations (105) and the vector system

$v$ $\displaystyle=$
$\displaystyle\mathbf{R}(\mathsf{q})\left[\cos{(w\tau)}\,\mbox{\boldmath$u$}+\sin{(w\tau)}\,\frac{\mbox{\boldmath$U$}}{w}\right],$
$V$ $\displaystyle=$
$\displaystyle\mathbf{R}(\mathsf{q})\left[-w\sin{(w\tau)}\,\mbox{\boldmath$u$}+\cos{(w\tau)}\,\mbox{\boldmath$U$}\right],$
(115)

where $\mathsf{q}$ is defined as in (95). Since the solutions for physical
time $v^{\ast}(\tau)$ and distance $r(\tau)$ do not depend on rotation of the
reference frame, we omit them – the readers may find them in Stiefel and
Scheifele (1971), Langner and Breiter (2015) or any other KS-related text.

It is common to select the reference frame rotation axis as
$\mbox{\boldmath$c$}=\mbox{\boldmath$e$}_{3}$. In that case, the solution is
further simplified, because then the rotation matrix $\mathbf{R}(\mathsf{q})$
does not influence $v_{3}$ and $V_{3}$. Thus the appropriate choice of the
rotation axis and of the defining vector leads to the simplest form of the
solution, where only two degrees of freedom, namely $v_{1},V_{1}$, and
$v_{2},V_{2}$, are affected by the rotation.

##  6 Concluding remarks

We dare to hope that pinpointing the presence of the defining vector in the
Kustaanheimo-Stiefel transformation may help in both the understanding and the
efficient use of this ingenuous device. We do encourage those of the readers
who practice the use of KS variables, to choose the defining vector best
suited for the problem at hand, instead of inertially following the ‘for
example’ choice made by Kustaanheimo and Stiefel. Actually, a common habit in
physics is to align the third axis with a symmetry axis of the potential. It
means that the KS3 set ($\mbox{\boldmath$c$}=\mbox{\boldmath$e$}_{3}$) should
be widespread, which is true in physics but not yet in celestial mechanics. Of
course (as noticed by a reviewer), from a purely formal point of view, the
adjustment of, say, KS1 to a different preferred direction $c$ may be achieved
by means of a rotation matrix $\mathbf{M}\in\mathrm{SO(4)}$, applied to the
left-hand side of (11) or (23). But then, in the quaternion formalism, one
should associate to $\mathbf{M}$ a unit quaternion $\mathsf{m}$, so that

$\alpha\,\mathsf{m}\mathsf{x}\bar{\mathsf{m}}=\mathsf{v}\mathsf{e}_{1}\bar{\mathsf{v}},$
(116)

and, finally

$\alpha\,\mathsf{x}=\bar{\mathsf{m}}\mathsf{v}\mathsf{e}_{1}\bar{\mathsf{v}}\mathsf{m}.$
(117)

The next step towards the form (31) is blocked by the lack of commutativity in
the quaternion product, leaving the relation of $\mathsf{m}$ to $c$ unclear,
and simplicity is lost, unless some trivial $\mathsf{m}$ has been considered.

Another way of using an arbitrary preferred direction is implicitly present in
the theory of generalized L-matrices worked out by Poleshchikov (2003), yet
its geometrical interpretation, analogous to the one we propose, would need an
additional effort.

Performing the canonical extension of the KS coordinates (point)
transformation, we went a step further than usual, providing the direct
expression of new momenta $\mathsf{V}$ in terms of $x$ and $X$. It has
revealed the dependence of the KS momenta on the Cartesian angular momentum
vector. The explicit relation between Fradkin tensor and Laplace vector,
derived in this paper, is another point of novelty, at least to our knowledge.

While working on some parts of the present study, we have been occasionally
surprised by the power of the quaternion algebra when applied to the KS
formulation of Deprit et al (1994). It was not our intention to contradict the
part of their conclusions that praised symbolic processors and lengthy
calculations, but some proofs happened to be shorter than expected and we
could not help it.

###### Acknowledgements.

Not by chance, the text was completed and submitted in the day marking the
tenth anniversary of the death of Prof. André Deprit. We would like to
dedicate the work to the memory of this unforgettable man of science and a
benevolent spirit of our Observatory. We thank reviewers for their comments
and bibliographic suggestions.

## References

  * Bellandi Filho and Menon (1987) Bellandi Filho J, Menon M (1987) The Kustaanheimo-Stiefel transformation in a spinor repesentation. Revista Brasileira de Fisica 17(2):302–309
  * Cartan (1966) Cartan E (1966) The Theory of Spinors. M.I.T. Press, Cambridge, Massachusetts
  * Cordani (2003) Cordani B (2003) The Kepler Problem. Springer Basel AG, Basel-Boston-Berlin
  * Deprit et al (1994) Deprit A, Elipe A, Ferrer S (1994) Linearization: Laplace vs. Stiefel. Celestial Mechanics and Dynamical Astronomy 58:151–201, DOI 10.1007/BF00695790
  * Duru and Kleinert (1979) Duru I, Kleinert H (1979) Solution of the path integral for the H-atom. Phys Let 84B(2):185–188
  * Fradkin (1967) Fradkin D (1967) Existence of the dynamic symmetries $O_{4}$ and $SU_{3}$ for all classical central potential problems. Progress in Theoretical Physics 37(5):798–812
  * Kurcheeva (1977) Kurcheeva IV (1977) Kustaanheimo-Stiefel regularization and nonclassical canonical transformations. Celestial Mechanics 15:353–365, DOI 10.1007/BF01228427
  * Kustaanheimo (1964) Kustaanheimo P (1964) Spinor regularization of the Kepler motion. Annales Universitatis Turkuensis, Series A 73:1–7
  * Kustaanheimo and Stiefel (1965) Kustaanheimo P, Stiefel E (1965) Perturbation theory of Kepler motion based on spinor regularization. J Reine Angew Math 218:204–219
  * Langner and Breiter (2015) Langner K, Breiter S (2015) KS variables in rotating reference frame. Application to cometary dynamics. Astrophysics and Space Science 357:153, DOI 10.1007/s10509-015-2384-6
  * Morais et al (2014) Morais JP, Georgiev S, Sprössig W (2014) Real Quaternionic Calculus Handbook. Birkhäuser, Basel, DOI 10.1007/978-3-0348-0622-0
  * Poleshchikov (2003) Poleshchikov SM (2003) Regularization of motion equations with L-transformation and numerical integration of the regular equations. Celestial Mechanics and Dynamical Astronomy 85:341–393
  * Roa et al (2016) Roa J, Urrutxua H, Peláez J (2016) Stability and chaos in Kustaanheimo-Stiefel space induced by the Hopf fibration. Monthly Notices RAS 459:2444–2454, DOI 10.1093/mnras/stw780, 1604.06673
  * Saha (2009) Saha P (2009) Interpreting the Kustaanheimo-Stiefel transform in gravitational dynamics. Mon Not RAS 400:228–231, DOI 10.1111/j.1365-2966.2009.15437.x, 0803.4441
  * Steane (2013) Steane AM (2013) An introduction to spinors. ArXiv e-prints 1312.3824
  * Stiefel and Scheifele (1971) Stiefel E, Scheifele G (1971) Linear and Regular Celestial Mechanics. Springer-Verlag, Berlin, Heidelberg, New York
  * Vivarelli (1983) Vivarelli MD (1983) The KS-transformation in hypercomplex form. Celestial Mechanics 29:45–50, DOI 10.1007/BF01358597
  * Vivarelli (1988) Vivarelli MD (1988) Geometrical and physical outlook on the cross product of two quaternions. Celestial Mechanics 41:359–370
  * Volk (1976) Volk O (1976) Miscellanea from the history of celestial mechanics. Celestial Mechanics 14:365–382, DOI 10.1007/BF01228523
  * Waldvogel (2006) Waldvogel J (2006) Quaternions and the perturbed Kepler problem. Celestial Mechanics and Dynamical Astronomy 95:201–212, DOI 10.1007/s10569-005-5663-7
  * Waldvogel (2008) Waldvogel J (2008) Quaternions for regularizing Celestial Mechanics: the right way. Celestial Mechanics and Dynamical Astronomy 102:149–162, DOI 10.1007/s10569-008-9124-y