The Lissajous-Kustaanheimo-Stiefel transformation

TL;DR

This paper introduces the Lissajous-Kustaanheimo-Stiefel (LKS) transformation, a new method to analyze perturbed Kepler problems using action-angle variables that avoid singularities and do not depend on the orbital plane.

Contribution

The LKS transformation provides a novel, singularity-free set of action-angle variables for the perturbed Kepler problem, enabling better analysis of orbital stability and dynamics.

Findings

LKS transformation yields singularity-free action-angle variables.

Allows analysis of stability for most orbital configurations.

Enables reduction to lower degrees of freedom using bilinear invariant.

Abstract

The Kustaanheimo-Stiefel transformation of the Kepler problem with a time-dependent perturbation converts it into a perturbed isotropic oscillator of 4-and-a-half degrees of freedom with additional constraint known as bilinear invariant. Appropriate action-angle variables for the constrained oscillator are required to apply canonical perturbation techniques in the perturbed problem. The Lissajous-Kustaanheimo-Stiefel (LKS) transformation is proposed, leading to the action-angle set which is free from singularities of the LCF variables earlier proposed by Zhao. One of the actions is the bilinear invariant, which allows the reduction back to the 3-and-a-half degrees of freedom. The transformation avoids any reference to the notion of the orbital plane, which allowed to obtain the angles properly defined not only for most of the circular or equatorial orbits, but also for the degenerate, rectilinear ellipses. The Lidov-Kozai problem is analyzed in terms of the LKS variables, which allow a direct study of stability for all equilibria except the circular equatorial and the polar radial orbits.

Tables

Table 1. Table 1: Particular orbits and their relation to the LKS variables.

Orbit type	LKS variables	Undetermined angles
Generic circular		none
Circular, polar	$\Lambda =0,G=0,\lambda =\left(2k+1\right)\frac{\pi }{4}$	none
Circular, equatorial	$\Lambda =0,\left|G\right|=L$	$l,g,\lambda$
Generic radial		none
Radial, equatorial	$G=0,\Lambda =0,\lambda =k\frac{\pi }{2}$	none
Radial, polar	$G=0,\left|\Lambda \right|=L$	$l,g,\lambda$
Generic equatorial		none

Equations

$$\mathsf{v}=\left(v_0,\vec{v}\right)=\sum _{j=0}^3{v}_j{\mathsf{e}}_j,$$

$${\mathsf{e}}_0=(1,\vec{0}),{\mathsf{e}}_1=(0,{\vec{e}}_1),{\mathsf{e}}_2=(0,{\vec{e}}_2),{\mathsf{e}}_3=(0,{\vec{e}}_3),$$

$$\mathsf{u}+\mathsf{v}=\sum _{j=0}^3\left(u_j+v_j\right){\mathsf{e}}_j,\alpha \mathsf{v}=\sum _{j=0}^3\alpha v_j{\mathsf{e}}_j,$$

$$\mathsf{u}\cdot \mathsf{v}=\sum _{j=0}^3{u}_j{v}_j=u_0{v}_0+\vec{u}\cdot \vec{v},$$

$$\mathsf{u}\mathsf{v}=\left(u_0{v}_0-\vec{u}\cdot \vec{v},u_0\vec{v}+v_0\vec{u}+\vec{u}\times \vec{v}\right).$$

$$\overline{\mathsf{v}}=(v_0,-\vec{v}),$$

$$\mathsf{u}\wedge \mathsf{v}=\frac{\mathsf{v}\bar{\mathsf{u}}-\mathsf{u}\bar{\mathsf{v}}}{2}=\left(0,u_0\vec{v}-v_0\vec{u}+\vec{u}\times \vec{v}\right),$$

$$\kappa :\mathbb{H}\to {\mathbb{H}}^{\prime }:\mathsf{v}\mapsto \mathsf{x}=\frac{\mathsf{v}\mathsf{c}\overline{\mathsf{v}}}{\alpha },$$

$$\alpha \vec{x}=\left(v_0^2-\vec{v}\cdot \vec{v}\right)\vec{c}+2\left(\vec{c}\cdot \vec{v}\right)\vec{v}+2v_0\vec{v}\times \vec{c}=\left(\vec{c}\cdot \vec{v}\right)\vec{v}+{\left[\mathsf{v}\wedge (\mathsf{v}\wedge \mathsf{c})\right]}^{♮},$$

$$\Vert \vec{x}\Vert =r=\frac{\mathsf{v}\cdot \mathsf{v}}{\alpha }.$$

$$\mathsf{q}(\varphi )=(\cos \varphi ,\sin \varphi \vec{c}),$$

$$|\mathsf{q}(\varphi )|$$

$$\mathsf{q}(\varphi )\mathsf{q}(\psi )$$

$${\left[\mathsf{q}(\varphi )\right]}^{-1}$$

$$\mathsf{q}(\varphi )\mathsf{c}\overline{\mathsf{q}}(\varphi )$$

$$\mathsf{px}\overline{\mathsf{p}}=\frac{(\mathsf{pv})\mathsf{c}\overline{(\mathsf{pv})}}{\alpha },$$

$$\mathsf{p}(\varphi )=(\cos \varphi ,\sin \varphi \hat{\vec{x}}),$$

$$\mathsf{p}(\varphi )\mathsf{v}=\mathsf{vq}(\varphi ),\text{mbox}and\mathsf{p}(\varphi )\mathsf{v}\overline{\mathsf{q}}(\varphi )=\mathsf{v}.$$

$$\hat{\vec{x}}\cdot \vec{v}=\vec{c}\cdot \vec{v},$$

$${\vec{v}}_s=\sqrt{\alpha r}\frac{\vec{c}+\hat{\vec{x}}}{||\vec{c}+\hat{\vec{x}}||},$$

$$\mathsf{v}={\mathsf{v}}_{\mathrm{s}}\mathsf{q}(-\varphi ),$$

$$v_0$$

$$\vec{v}$$

$$\tan \beta =\sqrt{\frac{1-\vec{c}\cdot \hat{\vec{x}}}{2}}\tan \varphi .$$

$$\mathsf{v}=\sqrt{\alpha r}\left(\sin \varphi ,\cos \varphi \vec{c}\right),$$

$$\mathcal{J}(\mathsf{v},\mathsf{w})=\left(\bar{\mathsf{v}}\wedge \bar{\mathsf{w}}\right)\cdot \mathsf{c}=-v_0\vec{w}\cdot \vec{c}+w_0\vec{v}\cdot \vec{c}+\left(\vec{v}\times \vec{w}\right)\cdot \vec{c}.$$

$$\mathsf{w}=\mathsf{u}(0,\vec{f})=(-\vec{u}\cdot \vec{f},u_0\vec{f}+\vec{u}\times \vec{f}),$$

$$\mathsf{v}=\xi \mathsf{u}+\eta \mathsf{w}=\mathsf{u}(\xi ,\eta \vec{f}),$$

$$\mathsf{x}=\kappa (\mathsf{v})=\frac{\mathsf{u}(\xi ,\eta \vec{f})\mathsf{c}(\xi ,-\eta \vec{f})\overline{\mathsf{u}}}{\alpha }.$$

$$(\xi ,\eta \vec{f})\mathsf{c}(\xi ,-\eta \vec{f})=\left({\xi }^2-{\eta }^2\right)(0,\vec{c})+2\xi \eta (0,\vec{f}\times \vec{c}),$$

$$\vec{x}=\frac{{\xi }^2-{\eta }^2}{\alpha }{\hat{\vec{x}}}_1+\frac{2\xi \eta }{\alpha }{\hat{\vec{x}}}_2,$$

$${\hat{\vec{x}}}_1$$

Full text

††journal: CELE

1 Introduction

The Kustaanheimo-Stiefel (KS) transformation is probably the most renowned regularization technique for the three-dimensional Kepler problem. In the planar case, the conversion of the Kepler problem into a harmonic oscillator has been known since Goursat (1889) and Levi-Civita (1906), but its extension to the three-dimensional problem took many decades of futile efforts. Finally, Kustaanheimo (1964) discovered that the way to the third dimension is not direct, but requires a detour through a constrained problem with four degrees of freedom. The KS transformation gained popularity in the matrix-vector formulation of Kustaanheimo and Stiefel (1965), but it is much easier to interpret and generalize in the language of quaternion algebra, very closely related with the original spinor formulation of Kustaanheimo (1964).

The most common use of the KS transformation is the numerical integration of perturbed elliptic motion, where many intricacies introduced by the additional degree of freedom can be ignored, although – as recently demonstrated by Roa et al (2016) – they can be quite useful in the assessment of a global integration error. Analytical perturbation methods for KS-transformed problems often follow the way indicated by Kustaanheimo and Stiefel (1965) and developed by Stiefel and Scheifele (1971): variation of arbitrary constants is applied to constant vector amplitudes of the KS coordinates and velocities. But those who want to benefit from the wealth of canonical formalism, require a set of action-angle variables of the regularized Kepler problem.

The first step in this direction can be found in the monograph by Stiefel and Scheifele (1971), where the symplectic polar coordinates are introduced for each separate degree of freedom. However, this approach does not account for degeneracy of the problem and thus is unfit for the averaging-based perturbation techniques. Moreover, no attempt was made to relate this set with the constraint known as the ‘bilinear invariant’, effectively reducing the system to 3 degrees of freedom. Both problems have been resolved by Zhao (2015), who proposed the ‘LCF’ variables (presumably named after Levi-Civita (1906) and Féjoz (2001)). In his approach, the motion in the KS variables is considered in an osculating ‘Levi-Civita plane’ (Deprit et al, 1994) as a two degrees of freedom problem. The third degree of freedom is added by the pair of action-angle variables orienting the plane. The redundant fourth degree is hidden in the definition of the Levi-Civita plane. The transformed Keplerian Hamiltonian depends on a single action variable, the other two actions being closely related to the angular momentum and its projection on the polar axis. Interestingly, the result is identical to the ‘isoenergetic variables’ found by Levi-Civita (1913) without regularization.

The LCF variables respect the degeneracy and bring the oscillations back to three degrees of freedom. Yet they possess a significant weakness: they are founded on the orientation of a plane determined by the angular momentum. Whenever the angular momentum vanishes (even temporarily), the angles become undetermined and equations of motion are singular. It turns out that seeking the proximity to the Delaunay variables, Zhao (2015) reintroduced the singularities of unregularized Kepler problem. Of course, some singularities are inevitable when the problem having spherical topology is mapped onto a torus of action-angle variables. But there is always some freedom in the choice of the singularities. Recalling that the main purpose of regularization is to allow the study of highly elliptic and rectilinear orbits, we find it worth an effort to construct the action-angle set that – unlike the LCF variables – is regular for this class of motions.

The main goal of the present work is to derive an alternative set of the action-angle variables which is not based upon the notion of an orbital plane (thus avoiding singularities when the orbit degenerates into a straight segment), and to test it on some well known astronomical problem. Section 2 introduces some preliminary concepts related to the KS coordinate transformation in the language of quaternions. We use its generalized form with an arbitrary ‘defining vector’ (Breiter and Langner, 2017), which helps to realize how the choice of the KS1 or KS3 convention allows or inhibits the use of the Levi-Civita plane in the construction of the action-angle sets. We have also benefited from the opportunity to polish and extend the geometrical interpretation given to the KS transformation by Saha (2009). In Section 3 we complement the KS coordinates with their conjugate momenta and provide the Hamiltonian function in the extended phase space as the departure point for further transformations. Section 4 builds the new action-angle set – the Lissajous-Kustaanheimo-Stiefel (LKS) variables. Two independent Lissajous transformations are followed by a linear Mathieu transformation. In Section 5 we show how to interpret the new variables not only in terms of the Lissajous ellipses, but also by the reference to the angular momentum and Laplace vectors of the Kepler problem. As an application, we discuss the classical Lidov-Kozai problem (Section 6), showing that stability of rectilinear orbits can be discussed directly in terms of the LKS variables, which has not been possible using the Delaunay or the LCF framework. Conclusions and future prospects are presented in the closing Section 7.

2 KS transformation in quaternion form

2.1 Quaternion algebra

Adhering to the convention used by Deprit et al (1994), we treat a quaternion $\mathsf{v}\in\mathbb{H}$ as union of a scalar $v_{0}$ and a vector $\vec{v}$,

[TABLE]

where the standard basis quaternions

[TABLE]

have been defined by referring to the standard vector basis $\vec{e}_{j}$. Downgrading a ‘pure quaternion’ $\mathsf{u}=(0,\vec{u})\in\mathbb{H}^{\prime}$ to a vector $\vec{u}\in\mathbb{R}^{3}$ requires application of the projection operator $\natural$, whose action on any quaternion is $\mathsf{v}^{\natural}=(v_{0},\vec{v})^{\natural}=\vec{v}$.

As members of the Euclidean linear space $\mathbb{R}^{4}$, quaternions admit the sum and product-by-scalar rules

[TABLE]

as well as the scalar product

[TABLE]

implying the norm $|\mathsf{v}|=\sqrt{\mathsf{v}\mathbf{\cdot}\mathsf{v}}=\sqrt{v_{0}^{2}+\|\vec{v}\|^{2}}$, where $\|\vec{v}\|=\sqrt{\vec{v}\cdot\vec{v}},$ to distinguish the norms in $\mathbb{R}^{3}$ and $\mathbb{R}^{4}$.

What makes four-vectors $\mathsf{u}$ and $\mathsf{v}$ the members of the quaternion algebra $\mathbb{H}$ over $\mathbb{R}$, is the noncommutative quaternion product definition

[TABLE]

Note that $\mathbb{H}^{\prime}$ is only a linear subspace, but not a subalgebra of $\mathbb{H}$, because the quaternion product of two pure quaternions may have a nonzero scalar part.

Two other useful operations to be defined are the quaternion conjugate

[TABLE]

allowing to write $|\mathsf{v}|^{2}=\mathsf{v}\overline{\mathsf{v}}$, and the quaternion cross product

[TABLE]

always resulting in a pure quaternion, and reducing to a standard vector cross product if $u_{0}=v_{0}=0$.

2.2 KS coordinates transformation

2.2.1 Generalized definition

In a recent paper (Breiter and Langner, 2017) we have proposed a generalized form of the standard KS transformation $\kappa$ that uses an arbitrary ‘defining vector’ $\vec{c}$ with a unit norm and its respective pure quaternion $\mathsf{c}=(0,\vec{c})$, so that

[TABLE]

or, equivalently,

[TABLE]

links the KS variables quaternion $\mathsf{v}$ with the original Cartesian coordinates $\vec{x}\in\mathbb{R}^{3}$, the latter being the vector part of a pure quaternion $\mathsf{x}=(0,\vec{x})$. A real, positive parameter $\alpha$ was introduced by Deprit et al (1994). They gave it the units of length, in order to allow the KS coordinates $v_{j}$ carry the same units as $x_{j}$. We adhere to this convention for a while, although other options will be presented in Section 3. With $|\mathsf{c}|=\|\vec{c}\|=1$, the KS transformation $\kappa$ admits the well known property

[TABLE]

2.2.2 Fibres

A non-injective nature of the KS map had been known since its origins, although only recently it has been considered more an advantage than a nuisance (Roa et al, 2016).

Let us introduce a quaternion-valued function of angle $\phi$

[TABLE]

with a number of useful properties, like

[TABLE]

and special values $\mathsf{q}(0)=\mathsf{e}_{0}$, $\mathsf{q}(\pi/2)=\mathsf{c}$. The property (15) clearly implies that the KS transformation (8) is only homomorphic: given some representative KS quaternion $\mathsf{v}$, all quaternions $\mathsf{v}\mathsf{q}(\phi)$ belonging to the fibre parameterized by $0\leqslant\phi<2\pi$, render the same vector $\vec{x}$, i.e. $\kappa(\mathsf{v})=\kappa(\mathsf{v}\mathsf{q}(\phi))$. Indeed, since (15) describes the rotation of vector $\vec{c}$ around the axis $\vec{c}$, the left hand side of the equality can be substituted for $\mathsf{c}$ in eq. (8), and then $\mathsf{v}\,\mathsf{c}\,\overline{\mathsf{v}}=(\mathsf{v}\mathsf{q})\,\mathsf{c}\,\overline{(\mathsf{v}\mathsf{q})}$, leading to the same $\mathsf{x}$.

On the other hand, one might ask about the possibility of generating the fibre through the left multiplication by some quaternion function. Multiplying both sides of equality in (8) by a quaternion $\mathsf{p}$ from the left and its conjugate from the right we find the condition

[TABLE]

where the left-hand side remains equal to $\mathsf{x}=(0,\vec{x})$ only if $\mathsf{p}$ is a function

[TABLE]

that rotates vector $\vec{x}$ around itself. Thus, given some representative KS quaternion $\mathsf{v}$, we can create the fibre $\mathsf{p}(\phi)\mathsf{v}$ parameterized by $0\leqslant\phi<2\pi$, such that $\kappa(\mathsf{v})=\kappa(\mathsf{p}(\phi)\mathsf{v})=\mathsf{x}$.

The action of the fibre generators (11) and (16) with the same argument $\phi$ is equivalent; direct computation demonstrates that

[TABLE]

A kind of symmetry between the defining vector $\vec{c}$ and the normalized Cartesian position vector $\hat{\vec{x}}$ implied by the form of $\mathsf{q}(\phi)$ and $\mathsf{p}(\phi)$ manifests also in the geometrical construction of the next section.

2.3 KS quaternions more geometrico

Given the transformation (9), let us polish the geometrical interpretation of the KS variables proposed by Saha (2009). Scalar multiplication of both sides of (9) by $\mathsf{v}$ leads to the basic relation

[TABLE]

with two unit vectors $\vec{c}$ and $\hat{\vec{x}}=\vec{x}/r$. This property, valid for any scalar part $v_{0}$, means that all quaternions $\mathsf{v}$ belonging to the fibre of given Cartesian vector $\vec{x}$, have vector parts $\vec{v}=\mathsf{v}^{\natural}$ forming the same angle with $\vec{c}$ and $\vec{x}$, hence lying in the symmetry plane of this pair of vectors. The plane, marked grey in Fig. 1, contains $\vec{c}+\hat{\vec{x}}$, and is perpendicular to $\vec{c}-\hat{\vec{x}}$. The norm $|\mathsf{v}|=\sqrt{r\alpha}$ is the upper bound on the length of $\vec{v}$, so the dashed circle in Fig. 1 has the radius $\sqrt{\alpha r}$. Setting $v_{0}=0$ in equation (9), we see that $\vec{x}$ is a linear combination of $\vec{c}$ and $\vec{v}$, so the three vectors must be coplanar. Accordingly, there are exactly two pure quaternions related with $\vec{x}$: $\mathsf{v}_{\mathrm{s}}=(0,\vec{v}_{\mathrm{s}})$, and $-\mathsf{v}_{\mathrm{s}}$, where

[TABLE]

is the ‘Saha-Kustaanheimo-Stiefel (SKS)vector’ of Breiter and Langner (2017).

The entire fibre $\mathsf{v}$ can be generated from $\mathsf{v}_{\mathrm{s}}$ by the application of the generator (11),

[TABLE]

leading to

[TABLE]

The latter of the formulae is a parametric equation of an ellipse with the major semi-axis $\sqrt{\alpha r}$ and the eccentricity $\sqrt{(1+\vec{c}\mathbf{\cdot}\hat{\vec{x}})/2}$. The ellipse is drawn with a solid line in Fig. 1. The position angle $\beta$ in the figure should not be confused with the parametric longitude $\phi$; the angles are related by the formula

[TABLE]

The line segment with arrowheads at both ends in Fig. 1, complements the length of $\vec{v}$ to the full value $\sqrt{\alpha r}$, so its length can be interpreted as the absolute value of the scalar part of $\mathsf{v}$.

Of course, the generic picture shown in Fig. 1 does not include the special case of the parallel $\vec{x}$ and $\vec{c}$. If $\hat{\vec{x}}=\vec{c}$, the fibre degenerates to the set of quaternions having vector part aligned with $\vec{c}$, i.e.

[TABLE]

with $\mathsf{v}_{\mathrm{s}}=(0,\sqrt{\alpha r}\,\vec{c})$. The eccentricity of the ellipse from Fig. 1 attains the value 1, so the ellipse degenerates into a straight segment. The shaded plane from the figure is no longer defined.

But if $\hat{\vec{x}}=-\vec{c}$, the situation is different. Observing that then the ellipse from Fig. 1 turns into a circle, we conclude that the fibre consists exclusively of the pure quaternions $\mathsf{v}=(0,\sqrt{\alpha r}\,\hat{\vec{f}})$, where $\hat{\vec{f}}$ is any vector orthogonal to $\vec{c}$.

2.4 Bilinear form $\mathcal{J}$ and LC planes

2.4.1 Definitions

The skew-symmetric bilinear form $\mathcal{J}:\mathbb{H}\times\mathbb{H}\rightarrow\mathbb{R}$, introduced by Kustaanheimo (1964) and discussed in later works, can be generalized to an arbitrary defining vector $\vec{c}$ as

[TABLE]

The form plays a central role in the KS formulation of motion. If the motion can be restricted to the linear subspace of $\mathbb{H}$ spun by two basis quaternions $\mathsf{u}$ and $\mathsf{w}$, such that $\mathcal{J}(\mathsf{u},\mathsf{w})=0$, the KS transformation reduces to the Levi-Civita transformation (Levi-Civita, 1906). For this reason, a two-dimensional subspace $P$ of quaternions being the linear combinations of $\mathsf{u}$ and $\mathsf{w}$, hence such that the form $\mathcal{J}$ on any two of them equals 0, was dubbed the ‘Levi-Civita plane’ by Stiefel and Scheifele (1971). We will use the name ‘LC plane’, although, strictly speaking, a (hyper-)plane in a space of dimension 4 should be spun by 3 basis quaternions.

Repeating the proof of Theorem 3 from Deprit et al (1994) in our generalized framework, we conclude that for any unitary quaternion $\mathsf{u}$ selected for the orthonormal basis of LC plane $P$, the second basis quaternion should be

[TABLE]

where $\vec{f}$ is any unitary vector orthogonal to the defining vector, i.e. $\vec{f}\cdot\vec{c}=0$, and $||\vec{f}||=1$. This meaning of the symbol $\vec{f}$ will be held throughout the text. The basis is indeed orthonormal, since $\mathsf{u}\cdot\mathsf{w}=0$, and $|\mathsf{u}|=|\mathsf{w}|=1$, by the definition of $\mathsf{u}$ and $\vec{f}$.

2.4.2 KS map of an LC plane

Once the LC plane has been defined, a question arises about the possibility of restricting the motion in KS variables to this subspace. But such restriction implies that the motion in ‘physical’ configuration space $\mathbb{R}^{3}$ is planar.

Let us prove that KS transformation maps any quaternion in the LC plane $P$ onto a plane $\Pi$ in ‘physical’ $\mathbb{R}^{3}$ space. Using the basis of two orthonormal quaternions $\mathsf{u}$ and $\mathsf{w}=\mathsf{u}(0,\vec{f})$, we consider their linear combination

[TABLE]

with real parameters $\xi,\eta$ having the dimension of length. The KS transform of these $\mathsf{v}$, belonging to $P$, is,by the definition (8),

[TABLE]

Thanks to the orthogonality of $\vec{c}$ and $\vec{f}$, the product in the middle evaluates to

[TABLE]

so the vector part of $\kappa{(\mathsf{v})}$ is a linear combination of two fixed, orthonormal vectors

[TABLE]

where

[TABLE]

Since eq. (30) is actually a parametric equation of a plane in $\mathbb{R}^{3}$, we have demonstrated that the KS transformation of any LC plane $P\subset\mathbb{H}$ is a plane $\Pi$ in $\mathbb{R}^{3}$ (or in $\mathbb{H}^{\prime}$, depending on the context). The parameters $\xi$, $\eta$ become parabolic coordinates on $\Pi$, i.e. the usual Levi-Civita variables.

The vector normal to the plane $\hat{\vec{x}}_{3}=\hat{\vec{x}}_{1}\times\hat{\vec{x}}_{2}$, can be most easily derived in terms of the quaternion cross product (7), with the first lines of (31) and (32) substituted. Thus, letting $\mathsf{x}_{j}=(0,\hat{\vec{x}}_{j})$, and $\mathsf{b}=(0,\vec{f}\times\vec{c})$, we find

[TABLE]

Thanks to the above equation, we can relate the choice of the LC plane basis to the orientation of $\Pi$. The cosine of the angle between the defining vector $\vec{c}$ and the normal to the plane of motion $\hat{\vec{x}}_{3}$ is given by the scalar product

[TABLE]

2.4.3 KS1 and KS3 setup

Some particular choices of the first basis quaternion $\mathsf{u}$ deserve a special comment. Inspecting eq. (34), we notice three obvious cases leading to $\vec{c}$ positioned in the plane of motion: a pure scalar $\mathsf{u}=\pm(1,\vec{0})$, or pure quaternions: $\mathsf{u}=(0,\pm\vec{f})$, and $\mathsf{u}=(0,\pm\vec{c})$. The basis vectors $\hat{\vec{x}}_{1}$, resulting from eq. (31), are $\vec{c}$, $-\vec{c}$, and $\vec{c}$, respectively. The last case, i.e. $\mathsf{u}=(0,\vec{c})$, has been the most common choice in celestial mechanics since the first paper of Kustaanheimo (1964). It allows the most direct identification of the LC plane with the plane of motion, both spanned by the same vectors (or pure quaternions) $\mathsf{u}^{\natural}=\hat{\vec{x}}_{1}=\vec{c}$, and $\mathsf{w}^{\natural}=\hat{\vec{x}}_{2}=\vec{c}\times\vec{f}$. The freedom of choice for $\vec{f}$ (any vector perpendicular to $\vec{c}$) permits to identify $\vec{c}$ and $-\vec{f}$ with the basis vector $\vec{e}_{1}$ and $\vec{e}_{3}$ of the particular reference frame used to describe the planar ($x_{3}=0$) motion. For this reason, let us call the KS transformation based upon the paradigmatic choice $\vec{c}=\vec{e}_{1}$, the KS1 transformation.

Remaining in the domain of pure quaternions, let us consider $\mathsf{u}=(0,\vec{u})$. Without loss of generality, we can assume $\vec{u}=\cos{\psi}\vec{c}+\sin{\psi}\vec{f}$, with $0\leqslant\psi\leqslant\pi$. Then, according to eq. (34), we have $\vec{c}\cdot\hat{\vec{x}}_{3}=\sin{2\psi}$, so an appropriate choice of the parameter $\psi$ may lead to any orientation of the orbital plane with respect to $\vec{c}$. In particular, the defining vector will coincide with $\hat{\vec{x}}_{3}$ when $\psi=\pi/2$. The LC plane spanned by the basis quaternions

[TABLE]

is mapped onto the plane of motion $\Pi$ with basis vectors $\hat{\vec{x}}_{1}=\vec{f}$, and $\hat{\vec{x}}_{2}=\vec{c}\times\vec{f}$ – both orthogonal to $\vec{c}$. Thus the choice of $\vec{c}=\vec{e}_{3}$, and $\vec{f}=\vec{e}_{1}$ leads to the KS3 transformation, which may look less attractive than KS1, with its LC plane no longer consisting of pure quaternions. Indeed, it is not practiced in celestial mechanics, save for two exceptions known to the authors (Saha, 2009; Breiter and Langner, 2017). In physics, however, the KS3 transformation is common at least since 1970’s (e.g. Duru and Kleinert, 1979; Cordani, 2003; Díaz et al, 2010; Egea et al, 2011; van der Meer et al, 2016); there are good reasons for this, but they come out only in the context of dynamics and symmetries of a perturbed Kepler (or Coulomb) problem.

3 Canonical KS variables in the extended phase space

In contrast to earlier works, let us consider from the onset a canonical problem in the extended phase space $(x^{\ast},\vec{x},X^{\ast},\vec{X})$, with a Hamiltonian

[TABLE]

where the Keplerian term

[TABLE]

depends on the Cartesian coordinates $\vec{x}$, their conjugate momenta $\vec{X}$ and the gravitational parameter $\mu$. The time-dependent perturbation $\mathcal{R}(t,\vec{x},\vec{X})$ is converted into a conservative term by substituting a formal, time-like coordinate $x^{\ast}$ for physical time $t$. The fact that $x^{\ast}(t)=t$ is a direct consequence of the way its conjugate momentum $X^{\ast}$ appears in equation (36), because

[TABLE]

and an appropriate choice of the arbitrary constant leads to the identity map of $t$ on $x^{\ast}$. The momentum $X^{\ast}$ itself evolves according to

[TABLE]

counterbalancing the variations of energy in nonconservative problems, or staying constant in the conservative case.

If the same problem is to be handled canonically in terms of the KS coordinates, their conjugate momenta $\mathsf{V}$ are implicitly defined through

[TABLE]

In this transformation we postulate

[TABLE]

to secure

[TABLE]

so that $\mathsf{X}=(0,\vec{X})$ remains a pure quaternion. The transformation $\mathbb{R}^{2}\times\mathbb{H}^{2}\rightarrow\mathbb{R}^{2}\times\mathbb{H}^{\prime}\times\mathbb{H}^{\prime}$, which maps $(v^{\ast},\mathsf{v},V^{\ast},\mathsf{V})\mapsto(x^{\ast},\vec{x},X^{\ast},\vec{X})$ according to the definitions (8), (40) and the identities $x^{\ast}=v^{\ast}$, $X^{\ast}=V^{\ast}$, is known to be weakly canonical (i.e. canonical only on a specific manifold (41)).

Let us now generalize the transformation by allowing that $\alpha$, instead of being a fixed parameter, is an arbitrary differentiable function of the energy-like momentum $X^{\ast}$ or $V^{\ast}$. A similar assumption was recently made for the Levi-Civita transformation (Breiter and Langner, 2018). The necessity, or at least usefulness of such a generalization will not be clear until the action-angle variables are introduced, but it has to be introduced already at this stage. If the generalized transformation is to be kept weakly canonical, while maintaining the direct relation $V^{\ast}=X^{\ast}$, the new formal time-like variable $v^{\ast}$ should differ from $x^{\ast}$.111Giving credit to previous applications of this idea in Breiter and Langner (2018), we have overlooked Stiefel and Scheifele (1971), much earlier than Zhao (2016). Then, the transformation

[TABLE]

conserves the Pfaffian one-form up to the total differential of a primitive function $Q$ (Arnold et al, 1997)

[TABLE]

provided

[TABLE]

and with a necessary condition of $\mathcal{J}(\mathsf{v},\mathsf{V})=0$.

It is worth noting, that with an elementary choice of $\alpha=k_{1}(X^{\ast})^{k_{2}}$, the expression in the square bracket evaluates to a single number $k_{2}$, and the multiplier $k_{1}$ has no influence on canonicity, hence it can be selected at will – for example to conserve (or to modify) the units of time and length.

In order to convert the Hamiltonian (36) into a perturbed harmonic oscillator, the independent variable has to be changed from the physical time $t$ to the Sundmann time $\tau$, related by

[TABLE]

involving $\alpha$ as a function of $V^{\ast}$ or $X^{\ast}$. Transforming the Hamiltonian (36) by the composition of $\lambda$ and $t\mapsto\tau$, we obtain

[TABLE]

where $\mathcal{R}^{\star}$ is the perturbation Hamiltonian $\mathcal{R}(x^{\ast},\vec{x},\vec{X})$ expressed in terms of the extended KS coordinates and momenta, and

[TABLE]

will have a constant value only if the original Hamiltonian $\mathcal{H}$ does not depend on time. Let us emphasize, that now every function of $\vec{x}$, when expressed in terms of $\mathsf{v}$, will generally depend on the energy-like momentum $V^{\ast}$ as well, due to its presence in $\alpha$. Noteworthy, the simplification to $\omega=1$ can be achieved by assuming $\alpha=\sqrt{8V^{\ast}}$, which makes the Sundmann time dimensionless. Choosing $\alpha=\mu/V^{\ast}$, is roughly equivalent to $\alpha=2a$, in terms of the Keplerian orbit semi-axis $a$.

4 Action-angle variables

4.1 LLC and LCF variables

When the motion is planar, with $x_{3}=0$, an appropriate action-angle set $l,g,L,G$ can be created using a combination of the Levi-Civita (Levi-Civita, 1906) and Lissajous transformations (Deprit and Williams, 1991). This approach has been recently revisited and discussed by Breiter and Langner (2018). Viewed as the special case of the KS framework, the Lissajous-Levi-Civita (LLC) variables are inherently attached to the KS1 setup, requiring the identification of the LC plane $P\subset\mathbb{H}^{\prime}$ of pure quaternions and the plane of motion $\Pi\subset\mathbb{R}^{3}$. A generalization of this approach was proposed by Zhao (2015). Roughly speaking, he attached the LC plane to an osculating plane of motion $\Pi$, and added the third action-angle pair $h,H$ orienting $\Pi$ in $\mathbb{R}^{3}$ by direct analogy with the third Delaunay pair: longitude of the ascending node, and projection of the angular momentum on the axis $\hat{\vec{x}}_{3}$. As noted by the author, this approach has the same drawbacks as in the Dalunay set – in particular, the singularity when the orbit in physical space is rectilinear, thus having no unique orbital plane.

4.2 Lissajous-Kustanheimo-Stiefel (LKS) variables

4.2.1 Intermediate set

The starting point for the new set of variables we would like to propose is completely different than in Zhao (2015). First, we choose the KS3 framework, assuming the defining vector $\vec{c}=\vec{e}_{3}$. Then, we select two subspaces of $\mathbb{H}$: $P_{03}$ with the basis $\mathsf{e}_{0},\mathsf{e_{3}}$, and $P_{12}$ spanned by $\mathsf{e}_{1}$ and $\mathsf{e}_{2}$. None of them is a Levi-Civita plane, because in the KS3 framework

[TABLE]

Thus, even for the planar case, we do not restrict motion to an invariant plane $P$, but merely project $\mathsf{v}$ on two orthogonal subspaces. The orthogonality is readily checked by

[TABLE]

On each plane, with $(i,j)=(0,3)$, or $(i,j)=(1,2)$, we perform the Lissajous transformation of Deprit (1991)

[TABLE]

Similarly to Breiter and Langner (2018), we allow $\omega>0$ to be a function of $V^{\ast}$, as given by equation (51) – both directly, and through $\alpha$. This requires a new time-like variable $s$ to be different from $v^{\ast}$, while retaining its conjugate $S=V^{\ast}$. Only then, the 1-forms are conserved up to the total differential

[TABLE]

with

[TABLE]

and

[TABLE]

where

[TABLE]

The Hamiltonian function (48) is converted into the sum of

[TABLE]

and of the perturbation $\mathcal{P}$ expressed in terms of the Lissajous variables.

This transformation is merely an intermediate step, but before the final move let us inspect the meaning and properties of the variables in the Kepler problem defined by $\mathcal{K}^{\prime}_{0}=0$. As a generic example we take a heliocentric orbit in physical phase space with the following Keplerian elements: major semi-axis $a=10\mathrm{au}$, eccentricity $e=0.5$, inclination $I=10^{\circ}$, argument of perihelion $\omega_{\mathrm{o}}=60^{\circ}$, longitude of the ascending node $\Omega=10^{\circ}$, and the initial true anomaly $f=60^{\circ}$. From these elements we compute first the position $\vec{x}(0)$ and momentum $\vec{X}(0)$, and then the representative KS3 quaternions $\mathsf{v}(0)$ and $\mathsf{V}(0)$ – an SKS vector given by equation (19), and its conjugate momentum defined by equation (40), both with $\mathsf{c}=\mathsf{e}_{3}$. These initial conditions are labeled with black dots in Fig. 2a. The ellipses described in the $(v_{0},v_{3})$ and $(v_{1},v_{2})$ planes have different semi-axes and different eccentricities; however, both are traversed in the same direction – retrograde (clockwise) in the discussed example. The retrograde motion follows from the fact that $G_{03}=G_{12}<0$ (the momenta are equal due to the postulate (41), where $\left(\bar{\mathsf{v}}\wedge\bar{\mathsf{V}}\right)\cdot\mathsf{e}_{3}=(G_{03}-G_{12})/2$). The constant angles $g_{03}$ and $g_{12}$, measured counterclockwise, position the ellipses in the coordinate planes. The initial angles $l_{03}$ and $l_{12}$ are marked according to the geometrical construction similar to that of the eccentric anomaly. Comparing our Fig. 2 with Fig. 1 of Deprit (1991), the readers may note the reverse direction of the $l_{ij}$ angle. The difference comes from the fact that Deprit (1991) assumed $G>0$, i.e. the prograde (counterclockwise) motion along the Lissajous ellipse. Yet, regardless of the sign of $G_{ij}$, equations of motion imply $\mathrm{d}l_{03}/\mathrm{d}\tau=\mathrm{d}l_{12}/\mathrm{d}\tau=\omega>0$.

Each Lissajous ellipse has the major semi-axis $a_{ij}$ and the minor semi-axis $b_{ij}$ defined by the two momenta and frequency

[TABLE]

The absolute value operator is necessary for $G_{ij}<0$, unless one adopts a convention of negative minor semi-axis for an ellipse traversed clockwise.

Another point worth observing is the ambiguity in the choice of the $l_{ij}$ and $g_{ij}$ pair. Their values are determined from two possible sets of four equations linking quadratic forms of $v_{i},v_{j},V_{i},V_{j}$ with sine and cosine functions of the angles. Regardless of whether we use

[TABLE]

or

[TABLE]

the solution will always result in two pairs: $(l_{ij},g_{ij})$, and $(l_{ij}+\pi,g_{ij}+\pi)$ – both giving the same values of the sine and cosine.222The statements about ‘the Lissajous variables […] determined unambiguously from the Cartesian variables’ made by Deprit (1991) should not be taken too literally. In other words, one of the two minor semi-axes in each of the ellipses in Fig. 2 can be chosen at will as the reference one.

Recalling the fibration property of the KS variables, we have plotted the ellipses obtained from the same Cartesian $\vec{x}(0)$ and $\vec{X}(0)$, but with the KS3 initial conditions $\mathsf{v}(0)$ and $\mathsf{V}(0)$ right-multiplied by $\mathsf{q}(\pi/2)=\mathsf{c}=\mathsf{e}_{3}$, according to equation (11) in the KS3 case. The results are displayed in Fig. 2b. Not only the initial conditions, abut the entire ellipses are rotated by $90^{\circ}$ in the $(v_{0},v_{3})$ plane, and by $-90^{\circ}$ in the $(v_{1},v_{2})$ plane. The momenta $L_{ij},G_{ij}$, and the angles $l_{ij}$ remain intact, compared to Fig. 2a. The new angles positioning the ellipses are $g^{\prime}_{03}=g_{03}+\pi/2$, and $g^{\prime}_{12}=g_{12}-\pi/2$, but their sum has not changed: $g^{\prime}_{03}+g^{\prime}_{12}=g_{03}+g_{12}.$

4.2.2 Final transformation

Bearing in mind the example shown in Fig. 2, we can establish the final set of the LKS variables by defining four action-angle pairs

[TABLE]

with $s$ and $S$ retained unaffected. One may easily verify that (64) amounts to an elementary Mathieu transformation, thus the complete composition

[TABLE]

is a weakly canonical, dimension raising transformation. The Hamiltonian $\mathcal{H}$ from equation (36) is transformed into

[TABLE]

where

[TABLE]

and $\mathcal{Q}$ is the pullback of $\frac{4r}{\alpha(S)}\mathcal{R}(x^{\ast},\vec{x},\vec{X})$ by $\zeta$.

Expressing the Cartesian variables from the initial extended phase space in terms of the LKS variables, we first introduce six actions-dependent coefficients

[TABLE]

allowing a compact formulation of the expressions for coordinates

[TABLE]

and momenta

[TABLE]

where

[TABLE]

Finally, the ‘time deputy’ variable $x^{\ast}=t$ is linked with the formal time-like variable $s$ through

[TABLE]

We have skipped the explicit expression of the KS variables, because it can be immediately obtained from the substitution of (64) into (54-57).

Two features of the above expressions for $\mathsf{x}$, $\mathsf{X}$, and $x^{\ast}$ deserve special attention. First, none of them depends on $\gamma$, which means that any dynamical system primarily defined in terms of $\vec{x}$, $\vec{X}$, and time, conserves the value of $\Gamma$. Secondly, the expressions for the Cartesian coordinates and momenta in the extended phase space do not depend on the particular choice of $\alpha(S)$ and $\omega(S)$; the choice affects only the form of the Hamiltonian $\mathcal{M}$.

5 LKS variables and orbital elements

Let us interpret the variables forming the LKS set – first the momenta, and then their conjugate angles – by showing their relation to the Keplerian elements or the Delaunay variables.

5.1 LKS momenta

Comparing equations (42) and (72) one immediately finds that

[TABLE]

when $\vec{c}=\vec{e}_{3}$, so observing that $\mathcal{J}(\mathsf{v},\mathsf{V})=0$ is the fundamental assumption of the KS transformation since the time of Kustaanheimo (1964), there is no other choice than $\Gamma=0$. Recalling the absence of its conjugate angle $\gamma$ in the Hamiltonian, $\Gamma=0$ is the integral of motion.

The meaning of $G$ becomes clear once we find the pull-back of the orbital angular momentum $\vec{G}_{\mathrm{o}}$ by $\zeta$, obtaining

[TABLE]

Thus, setting $\Gamma=0$, we find the momentum $G$ to be twice the projection of the orbital angular momentum on the third axis (i.e. twice the Delaunay action $H_{\mathrm{o}}$). Whenever the Hamiltonian admits the rotational symmetry around $\vec{e}_{3}$, the momentum $G$ will be the first integral of the system.

Proceeding to the momentum $L$, we have to distinguish the pure Kepler problem and the perturbed one. In the former case, we can set $\mathcal{M}_{0}=0$ in equation (66), finding

[TABLE]

at $\Gamma=0$. Moreover, in the pure Kepler problem, the momentum $S$ can be expressed in terms of the major semi-axis $a$ as $S=\mu/(2a)$, which justifies the direct link between the values of $L$ and of the Delaunay action $L_{\mathrm{o}}$

[TABLE]

The two restrictive clauses of the previous sentence (‘values’ and ‘pure Kepler’) deserve comments. Equation (78) does not imply differential relations, because, for example, $\partial\vec{x}/\partial L_{\mathrm{o}}\neq 2\partial\vec{x}/\partial L$ (c.f. Deprit and Williams, 1991). Moreover, the values of $L$ and $2L_{\mathrm{o}}$ generally differ in a perturbed problem due to the fact that $L_{\mathrm{o}}$ is always defined by $\mathcal{H}_{0}$ alone, whereas the definition of LKS momentum $L$ depends on the complete Hamiltonian $\mathcal{H}_{0}+\mathcal{R}$ through the value of $S=X^{\ast}$ (the latter fixed by the restriction to the manifold $\mathcal{H}=0$).

Similar intricacies are met for the momentum $\Lambda$, which turns out to be related with the Laplace (eccentricity) vector $\vec{e}$, or rather the Laplace-Runge-Lenz vector $\vec{J}=L_{\mathrm{o}}\vec{e}$, having the dimension of angular momentum. In the pure Kepler problem, substituting $\Gamma=0$, we find

[TABLE]

Thus the momentum $\Lambda$ has been identified as twice the projection of the Laplace-Runge-Lenz vector on the third axis, yet this equality, using the property $2SL_{\mathrm{o}}^{2}=\mu^{2}$, holds only in the pure Kepler problem. In the perturbed case, one should refer to the general definition of $\vec{e}$ in terms of the KS variables (Breiter and Langner, 2017).

Closing the discussion of the momenta, let us collect the bounds on their values:

[TABLE]

By the construction, the value of $L=L_{12}+L_{03}$ must be nonnegative; but $L=0$ implies the permanent location at the origin ($\vec{x}=\vec{X}=\vec{0}$), so we exclude it. The momenta $\Lambda$ and $G$ may be either positive or negative, but the above inequality guarantees that all coefficients in equations (67) are real.

5.2 LKS angles

As already mentioned, the angle $\gamma$ is a cyclic variable, absent in the pullback of any Hamiltonian $\mathcal{H}$ by $\zeta$. Actually, $\gamma$ is the ‘KS angle’ parameterizing the the fibre of KS variables $(\mathsf{v},\mathsf{V})$ mapped into the same point in the $(\vec{x},\vec{X})$ phase space. Thus, unless we are interested in some topological stability issues (Roa et al, 2016), the angle can be ignored.

The only fast angular variable is $l$. As expected, its values in the pure Kepler problem are equal to a half of the orbital eccentric anomaly $E$. Indeed,

[TABLE]

and both angles are equal to 0 at the pericentre. Once again, this direct relation does not survive the addition of the perturbation. Nevertheless, it also reveals the nature of equation (74) as a generalized Kepler’s equation.

The two remaining angles are more unusual. A quick look at equations (76) and (79) might suggest that the role of $g$ and $\lambda$ in $\vec{G}_{\mathrm{o}}$ and $\vec{J}$ is similar. But if the norms of the vectors are evaluated, one finds

[TABLE]

The absence of $g$ proves it to be some rotation angle; the presence of $\lambda$ means that this angle plays a different role (and is somehow related with the eccentricity $e$).

More light is shed on this problem if we introduce the vectors

[TABLE]

These are essentially the so-called Cartan or Pauli vectors (Cordani, 2003), except that we use the sign of $\vec{N}$ opposite to the usual convention. Both the vectors have the same norm $M=N=L/4=L_{\mathrm{o}}/2$, and lie either in the plane perpendicular to orbit, or along a degenerate radial orbit direction. The angle $\theta$ they form depends on the eccentricity alone, because

[TABLE]

Obviously, $\theta$ is the upper bound for the angle $\theta^{\prime}$ between the projections of the Cartan vectors on the coordinate plane $(x_{1},x_{2})$

[TABLE]

Using equations (84), (85), and (87), one finds

[TABLE]

Let us make $\theta^{\prime}$ an oriented angle by postulating that it is measured from $\vec{N}^{\prime}$ to $\vec{M}^{\prime}$, counterclockwise. Then its sine is given by

[TABLE]

Thus we have identified the angle $\lambda$ as the quarter of the angle between the projections of the Cartan vectors on the reference plane $(x_{1},x_{2})$, measured from $\vec{N}^{\prime}$ to $\vec{M}^{\prime}$. Finding $\theta^{\prime}$ from the eccentricity-dependent $\theta$ involves orbital inclination and the argument of pericentre, which means that $\lambda$ is a function of $e$, $I$, and $\omega_{\mathrm{o}}$.

Interestingly, whenever the argument of pericentre $\omega_{\mathrm{o}}$ exists, the statement $\sin{4\lambda}=0$ means $\cos{\omega_{\mathrm{o}}}=0$. Thus, any $\lambda=k\pi/4$ refers to $\omega_{\mathrm{o}}=\pi/2$ or $\omega_{\mathrm{o}}=3\pi/2$.

Once we have interpreted $\lambda$, the meaning of $g$ comes out of equations (84) and (85): let us create the sum of normalized vectors $\vec{M}^{\prime}/\|\vec{M}^{\prime}\|+\vec{N}^{\prime}/\|\vec{N}^{\prime}\|$ and let us rotate the resulting vector by $\pi/2$ counterclockwise, obtaining

[TABLE]

This formula suggests that $g$ is a half of the longitude of $\vec{M}_{\mathrm{m}}$, or of $-\vec{M}_{\mathrm{m}}$, depending on the sign of $\cos{2\lambda}$. Whichever the case, changing the value of $g$ we perform a simultaneous rotation of both $\vec{N}^{\prime}$ and $\vec{M}^{\prime}$ by the same angle. Indirectly, it means the rotation of the orbital plane (if it exists) around the third axis, which makes $g$ a relative of the ascending node longitude.

5.3 Special orbit types

Let us inspect how some specific orbit types are mapped onto the LKS variables. The discussion is restricted to the elliptic orbits ($0\leqslant e\leqslant 1$) in the pure Kepler problem.

5.3.1 Circular orbits

Circular orbits, having $e=0$, are characterized by $\vec{J}=\vec{0}$, hence all must possess $\Lambda=0$, since $2\Lambda=\vec{J}\cdot\vec{e}_{3}$. Then, the norm of the Laplace-Runge-Lenz vector (83) simplifies, thanks to $C_{1}=C_{2}=\sqrt{L^{2}-G^{2}}/2$, and equating its square to 0 we find the condition

[TABLE]

Setting $\lambda=(2k+1)\pi/4$, $k\in\mathbb{Z}$, leads to generic circular orbits with the inclination $I=\arccos{(G/L)}$, including circular polar orbits when $G=0$. However, if $|G|=L$, then the first factor is null regardless of $\lambda$. This is the case of circular orbits in the ‘equatorial plane’ $(x_{1},x_{2})$: prograde for $G=L$, or retrograde for $G=-L$.

The values of $\lambda$ mentioned above well coincide with the interpretation from section 5.2. In circular orbits, the Cartan vectors $\vec{N}$ and $\vec{M}$ are collinear and opposite, thus the angle $\theta=\pi$, and its projection $\theta^{\prime}$ remains $\pm\pi$ as long as the orbit is not equatorial. Thus $\lambda=\theta^{\prime}/4=\pm\pi/4$, plus any multiple of $(2\pi)/4$.

Another explanation of the LKS variables for $e=0$ can be given by inspecting the Lissajous ellipses in Fig. 2. The orbital distance $r$ is the sum of $\rho^{2}_{03}=v_{0}^{2}+v_{3}^{2}$ and $\rho^{2}_{12}=v_{1}^{2}+v_{2}^{2}$, both divided by $\alpha$. In order to secure a constant $r=(\rho^{2}_{12}+\rho^{2}_{03})/\alpha$, it is not necessary that both $\rho_{ij}$ are constant; enough if they oscillate with the same amplitude and a phase shift of $\pm\pi/2$. Equal amplitudes result from $L_{12}=L_{03}$ (because $G_{12}=G_{03}$ by $\Gamma=0$), hence $\Lambda=L_{12}-L_{03}=0$. The phase shift condition is given by $l_{12}-l_{03}=2\lambda=(2k+1)\pi/2$, which means the values of $\lambda$ as above.

The case of constant $\rho_{ij}$, mentioned above, should be related with some special kind of a circular orbit. Indeed, since it needs $L_{12}=|G|/2=L_{03}$, i.e. two circles of equal radii in Fig. 2, we obtain the circular equatorial orbits with $\Lambda=0$ and $|G|=L$ (prograde or retrograde, depending on the sign of $G$). Observe that due to the lack of distinct semi-axes in the two circles, the angles $l_{ij}$, and $g_{ij}$ are undefined, and so are $l$, $g$, $\gamma$, and $\lambda$. But still one can use properly defined ‘longitudes’ $l+g$ or $l-g$ in the prograde, and retrograde cases, respectively – at least until some ‘virtual singularities’ appear (Henrard, 1974). In terms of the Cartan vectors, $\vec{M}=-\vec{N}=\vec{G}_{\mathrm{o}}$, so $\vec{M}^{\prime}=\vec{N}^{\prime}=\vec{0}$, making the angles $g$ and $\lambda$ undetermined.

5.3.2 Radial orbits

Rectilinear (radial) orbits require $\vec{G}_{\mathrm{o}}=\vec{0}$, hence $G=0$ and $\vec{G}_{\mathrm{o}}\cdot\vec{G}_{\mathrm{o}}=0$. According to equations (67), $G=0$ means $C_{1}=C_{2}=\sqrt{L^{2}-\Lambda^{2}}/2$, wherefrom equation (82) implies

[TABLE]

Regardless of $\lambda$, it is satisfied by $|\Lambda|=L$, i.e. by polar radial orbits with $\vec{J}=(\Lambda/2)\,\vec{e}_{3}$. For all other directions of the Laplace-Runge-Lenz vector, radial orbits need $\lambda=k\pi/2$, where $k\in\mathbb{Z}$; this time $\Lambda$ can be arbitrary, with $\Lambda=0$ indicating an equatorial radial orbit. In terms of the Cartan vectors, $\vec{G}_{\mathrm{o}}=0$ means $\vec{N}=\vec{M}$, so their angle $\theta=0$ is projected as $\theta^{\prime}=0+2k\pi$, which (divided by 4) gives the above values of $\lambda$.

In terms of the Lissajous ellipses in $(v_{1},v_{2})$ and $(v_{0},v_{3})$ planes from Fig. 2, $G=0$ means that both degenerate into straight segments. The motion along the segments must obey $l_{12}=l_{03}+k\pi$, to guarantee that $v_{0}=v_{1}=v_{2}=v_{3}=0$ at the same epoch. The direction of $\vec{x}(\mathsf{v})$ is determined by the difference of lengths of the two segments: equatorial orbits result if the segments have the same length, whereas polar orbits require that one of the segments collapses into a point. In the latter case, $l$ and $\lambda$ are undetermined, but $l+\lambda=l_{02}$ or $l-\lambda=l_{03}$ retain a well defined meaning for an appropriate sign of $\Lambda$. Problems with the definition of $g_{ij}$ due to the vanishing minor axes are only apparent, because they can be solved by an alternative definition: instead of ‘position angle of the minor semi-axis’, one can equally well say ‘position angle of the major semi-axis minus $\pi/2$’.

5.3.3 Equatorial orbits

Since the Laplace-Runge-Lenz vector lies in the plane $(x_{1},x_{2})$ for the equatorial orbits, all they must have $\Lambda=0$, whence $C_{1}=C_{2}=\sqrt{L^{2}-\Lambda^{2}}/2$. Moreover, $G_{\mathrm{o}}^{2}=(G/2)^{2}$, which leads to the condition

[TABLE]

The case of $|G|=L$ brings us back to the circular equatorial orbits, already discussed. Other values of $G$ require $\lambda=k\pi/2$, where $k\in\mathbb{Z}$. These are the same values as in the case of radial orbits, which makes sense, because $G=0$ should bring us to the radial equatorial orbit.

For an elliptic ($e\neq 0$) equatorial orbit, the Cartan vectors $\vec{N}$ and $\vec{M}$ may form different angles $\theta$, but since they lie in a polar plane, the projection of these angle is always $\theta^{\prime}=0$, exactly as in the radial orbit case – thus the same values of $\lambda$.

The two Lissajous ellipses in Fig. 2 must have the same semi-axes, and $l_{12}=l_{03}+k\pi$. This is necessary to obtain $v_{1}^{2}+v_{2}^{2}=v_{0}^{2}+v_{3}^{2}$, which guarantees $x_{3}=0$ for all epochs, according to equation (9) in the KS3 setup.

5.3.4 Polar orbits

Polar orbits are generically indicated by the simple condition $G=0$. It is only in the special cases where the angle $\lambda$ comes into play: circular polar orbits ($\Lambda=0$) need $\lambda=(2k+1)\pi/4$, whereas radial polar orbits $(|\Lambda|=L)$ are the ones where $\lambda$ is undetermined. Since $G=0$, both Lissajous ellipses degenerate into segments, but their lengths may be different, and the phase shift arbitrary.

5.3.5 Singularities

Inspecting specific types of orbits we met the situations, where $\lambda$ and $g$ become undetermined: circular equatorial orbits with $|G|=L,\Lambda=0$ and rectilinear polar orbits with $|\Lambda|=L,G=0$. These four points are the vertices of the square on the $(G,\Lambda)$ plane defined by the constraint $|G|+|\Lambda|\leqslant L$. However, all four edges of the square leave the angles undetermined. This is related to the fact, that:

• a)

$L=G+\Lambda$ (upper right edge in Fig. 4) means $L_{03}=G_{03}$, i.e. prograde circular motion on $(v_{0},v_{3})$ plane with undetermined $l_{03}$ and $g_{03}$ (but $l_{03}+g_{03}$ is well defined),

• b)

$L=-G+\Lambda$ (upper left edge in Fig. 4) means $L_{03}=-G_{03}$, i.e. retrograde circular motion on $(v_{0},v_{3})$ plane with undetermined $l_{03}$ and $g_{03}$ (but $l_{03}-g_{03}$ is well defined),

• c)

$L=G-\Lambda$ (lower right edge in Fig. 4) means $L_{12}=G_{12}$, i.e. prograde circular motion on $(v_{1},v_{2})$ plane with undetermined $l_{12}$ and $g_{12}$ (but $l_{12}+g_{12}$ is well defined),

• d)

$L=-G-\Lambda$ (lower left edge in Fig. 4) means $L_{12}=G_{12}$, i.e. retrograde circular motion on $(v_{1},v_{2})$ plane with undetermined $l_{12}$ and $g_{12}$ (but $l_{12}-g_{12}$ is well defined).

The Keplerian orbits obtained by mapping the edges of the $(G,\Lambda)$ square onto $\vec{J}$ and $\vec{G}_{0}$ or $\vec{x}$ and $\vec{X}$, have $e=\sin{I}$, and $\sin\omega_{\mathrm{o}}=\pm 1$. Thus the vertices in Fig. 4 are $\sin{I}=e=0$ (left and right) and $\sin{I}=e=1$ (top and bottom). Along the edges, half of the coefficients (67) does vanish, and one of the vanishing coefficients is always $C_{1}$ or $C_{2}$, which implies that either $\vec{M}^{\prime}$ or $\vec{N}^{\prime}$ is a null vector, so the angles $\lambda$ and $g$ become undefined.

6 Application to the Lidov-Kozai problem

6.1 Derivation of the secular model

In order to test the LKS variables in a nontrivial astronomical problem, let us revisit the Lidov-Kozai resonance arising in the artificial satellites theory (Lidov, 1962) or asteroid dynamics (Kozai, 1962). In this already classical problem, the Keplerian motion of a small body (a satellite or an asteroid) around a central mass with the gravitational parameter $\mu$ (a planet or the Sun) is influenced by a distant perturber with the gravitational parameter $\mu^{\prime}$ (the Sun or a planet, respectively). The origin of the reference frame is attached to the central mass, the plane $(x_{1},x_{2})$ coincides with the orbital plane of the perturber, and the third axis basis vector $\vec{e}_{3}$ is directed along the angular momentum of the perturber. Further, let us assume that the perturber moves on a circular orbit with the mean motion $n_{\mathrm{p}}$, so its position vector is

[TABLE]

Compared to the small body, whose position vector is $\vec{x}$, the perturber is distant, i.e. $||\vec{x}||/||\vec{r}_{\mathrm{p}}||=r/a_{\mathrm{p}}$ is small enough to approximate the perturbing function by the second degree Legendre polynomial term. Thus we obtain the problem with the Hamiltonian $\mathcal{H}$ from equation (36) with the perturbation

[TABLE]

The perturbation is time-dependent, so – after substituting (94) – we replace $t$ by its formal twin $x^{\ast}$, obtaining

[TABLE]

Choosing $\omega=1$, and $\alpha=\sqrt{8X^{\ast}}=\sqrt{8S}$, we apply the LKS transformation, setting $\Gamma=0$, because we are not interested in the evolution of the KS angle $\gamma$. The resulting Hamiltonian (65) is

[TABLE]

with the Keplerian part

[TABLE]

For a while, the perturbation $\mathcal{Q}$ will be given in an intermediate form without the explicit substitution of the LKS variables into $\vec{x}$ and $r$, which leads to a relatively concise form

[TABLE]

where, according to equation (74),

[TABLE]

By the choice of $\alpha$, the Sundman time $\tau$ is dimensionless and the unperturbed motion gives

[TABLE]

with all the remaining variables constant. Solving (101) we find $l=\tau+l_{0}$. The value of $S$ is set to give $\mathcal{M}=0$, but ignoring the contribution of $\mathcal{Q}$ we may estimate that $s\approx\tau/n$, where $n$ is the Keplerian mean motion.

According to the standard Lie transform method (e.g. Ferraz-Mello, 2007), the mean variables can be introduced by a nearly canonical transformation that converts $\mathcal{M}$ into $\mathcal{N}=\mathcal{N}_{0}+\mathcal{Q}^{\prime}$, with $\mathcal{N}_{0}=\mathcal{M}_{0}$ and $\mathcal{Q}^{\prime}$ being constant along the phase trajectory generated by $\mathcal{N}_{0}$. Up to the first order, the new perturbation $\mathcal{Q}^{\prime}$ is simply the average of $\mathcal{Q}$ with respect to $\tau$, assuming $l=\tau+l_{0}$ and $s=\tau/n$.

Since the perturber has been assumed distant, its mean motion $n_{p}$ is small compared to $n$ and both frequencies can be treated as irrational; even if they are not, the resonance will occur in high degree harmonics with practically negligible amplitudes. In these circumstances, any product of sine or cosine of $2n_{\mathrm{p}}s=2(n_{\mathrm{p}}/n)\tau$ with a function which is either constant or $2\pi$-periodic in $\tau$, has the zero average.333The general definition of the average for a function $f(\tau)$ is $\lim_{\tau\rightarrow\infty}\tau^{-1}\int_{0}^{\tau}f(\tau^{\prime})\mathrm{d}\tau^{\prime}$, so its value for a quasi-periodic function is null. When $f(\tau)$ is $T$-periodic, this definition simplifies to the standard $T^{-1}\int_{0}^{T}f(\tau)\mathrm{d}\tau$. Then $\mathcal{Q}^{\prime}$ simplifies to

[TABLE]

Thus we obtain the first order approximation of the secular system

[TABLE]

where the mean variables should be given different symbols, but we adhere to a widespread habit of distinguishing the mean and the osculating variables by context. The following study of motion generated by $\mathcal{N}$ will refer only to the mean variables, so no confusion should occur.

Both $\mathcal{N}$ and the classical secular Hamiltonian of the Lidov-Kozai problem share the same property: they are reduced to 1 degree of freedom. In our case it is the canonically conjugate pair $(\lambda,\Lambda)$ instead of the usual Delaunay pair of the argument of pericentre and the angular momentum norm. All other momenta are constant and will be treated as parameters. However, there is a fundamental difference between our formulation and the classical approach: the equations of motion for $\lambda$ and $\Lambda$ are not singular for most of the radial orbits.

6.2 Secular motion and equilibria

Let us set

[TABLE]

The equations of motion derived from (104) are

[TABLE]

Integral curves of this system are plotted in Fig. 3 for three values of $G$: $0.9L$, $0.75L$ and [math]. The phase plane has been clipped to $-\pi\leqslant\lambda\leqslant\pi$, because the reaming range of $\lambda$ is a simple duplication of the plotted phase portrait.

Referring to Table 1, one can check that a generic radial orbit does not introduce a singularity into equations (107). Indeed, $G=0$, and $\lambda=k\pi/2$ result in a well defined

[TABLE]

It means that a radial orbit is not en equilibrium, unless $\Lambda=0$, which is exactly the case of a radial orbit in the equatorial plane. Observing that for $G=0$ the points $(\lambda=k\pi/2,\Lambda=0)$ are well defined local minima of $\mathcal{N}$, we are able to state that radial orbits in the orbital plane of the perturber are stable444The word ‘stable’ is a bit paradoxical in this context, because it means that the motion starting in such an orbit will inevitably end up in collision with the central body. equilibria.The bottom panel of Figure 3 confirms this observation: the points $(0,0)$, $(90^{\circ},0)$, and $(-90^{\circ},0)$ are surrounded by closed, oval shape contours. Intersection of any of the integral curves plotted in the bottom panel with the vertical lines $\lambda=0$, or $\lambda=\pm 90^{\circ}$ marks a temporary passage through the radial orbit degeneracy.

As far as the polar radial orbits (with $|\Lambda|=L$) are concerned, equations (107) become singular, but this singularity is purely geometrical. Such orbits should be located at the upper and lower edges of the bottom panel in Fig. 3, where $\lambda$ is undetermined. But since the integral curves approaching the edges become parallel to them, one should expect that polar radial orbits are stable equilibria (which is actually the case, if the analysis is performed in terms of vectors $\vec{G}_{\mathrm{o}}$ and $\vec{J}$, or simply observing that for $G=0$ the Hamiltonian $\mathcal{N}$ has the local maxima at $\Lambda=\pm L$ regardless of the value of $\lambda$).

Actually, The presence of $\Lambda$ as a factor of the first of equations (107) means that for any value of $|G|\neq L$, the equilibria exist at $(\lambda=j\,\pi/4,\Lambda=0)$, as seen in Fig. 3. For even $j=2k$, the equilibria refer to equatorial orbits with the eccentricity depending on $G$ through $e=\sqrt{1-(G/L)^{2}}$. It is easy to check the they are the local minima of the Hamiltonian $\mathcal{N}$, hence the equatorial orbits are stable. The circular equatorial case with $|G|=L$ is problematic, because then the upper and lower limits of $\Lambda$ merge, and in order to prove that these are actually the stable equilibria one has to resort to the analysis of $\vec{G}_{\mathrm{o}}$ and $\vec{J}$ vectors.

For odd $j=(2k+1)$, the equilibria are circular orbits with inclinations depending on $G$ (equatorial if $G=0$, prograde for $G>0$ and retrograde when $G<0$). Their stability depends on the ratio $G/L$. Unlike in the Delaunay chart, variational equations can be formulated directly in the phase plane of $(\lambda,\Lambda)$, leading to the eigenvalues that are pure imaginary for $(G/L)^{2}>3/5$. Thus circular orbits are stable for inclinations below $I_{1}=\arccos{\sqrt{3/5}}\approx 39^{\circ}\!.23$ and above $I_{2}=\arccos{-\sqrt{3/5}}\approx 140^{\circ}\!.77$. At these critical values a bifurcation occurs: when $(G/L)^{2}<3/5$ circular orbits become unstable and two stable equlibria are created at $(\lambda=(2k+1)\,\pi/4,\Lambda=\pm\Lambda_{\mathrm{c}})$ (see the middle panel of Fig. 3). Recall that, in general case of inclined, elliptic orbits, this value of $\lambda$ means the argument of pericentre equal $\pi/2$ or $3\pi/2$. The value of $\Lambda_{\mathrm{c}}$ is the root of the first of equations (107) with $\Lambda\neq 0$ and $\cos{4\lambda}=-1$, i.e.

[TABLE]

leading to

[TABLE]

These are the classical equilibria of the Lidov-Kozai problem – the only ones that can be analyzed directly in the Delaunay variables. In terms of the orbital elements, equation (111) is equivalent to the well known condition (Lidov, 1962)

[TABLE]

The equilibria can be located in the middle panel of Fig. 3 at $\lambda=\pm 45^{\circ}$, $\Lambda\approx\pm 0.115\,L$.

Figure 4 shows all the equilibria and their stability, with the dashed lines marking the unstable equilibrium. The edges of the $(G,\Lambda)$ square (upper and lower boundaries of the plots in Fig. 3) may not be attached to any of the values of $\lambda$, but we added the black dots at the corners to show the stable equilibria of the special type as the natural limits of the stable branches (solid lines).

It is not unusual that all action-angle-like variables with bounded momentum suffer from indeterminate angle at the boundary of its conjugate. The LKS variables cannot be different, even if many cases, problematic in the Delaunay chart, have been located inside the boundaries of $\Lambda$. For each value of $G\neq 0$ (and $|G|\neq L$) there exist integral curves passing through both the extremes: $\Lambda=L-|G|$ and $\Lambda=-L+|G|$. In the top or the middle panel of Fig. 3 they are seen as four disjoint fragments; for example, the two open curves approaching the edges at $\lambda\pm 22^{\circ}\!.5$ are the fragments of such an integral curve. There is no singularity in these orbits (see Sect. 5.3.5) other than the indeterminacy of longitude at the poles of a sphere (Deprit, 1994).

7 Conclusions

While commenting a transformation due to Fukushima, Deprit (1994) observed that it amounts to swapping singularities, and immediately added ‘This remark is not meant to diminish its practical merit, quite the contrary.’ The LKS variables we have presented also ‘trade in singularities’, but the rule of trade we propose is to spare the radial, rectilinear orbits (except the polar ones) at the expense of some other types. The exceptions include mostly a family of expendable, rank-and-file orbits with $e=\sin{I}$ and the lines of apsides perpendicular to the lines of nodes – the cases easily tractable without the KS regularization and unlikely to focus attention by becoming equilibria in typical problems of celestial mechanics. More we regret the problems caused by the polar radial, and equatorial circular orbits. Nevertheless, we believe that more has been gained than lost. Enough to enumerate the orbits that remain regular points in our chart: circular inclined, equatorial elliptic, and all radial (except the polar ones). Thanks to refraining from the use of orbital plane in their construction, the LKS variables are better fitted to study highly elliptic orbits than any other action-angle set known to the authors.

The analysis of the quadrupole Lidov-Kozai problem in Section 6 suggests that the LKS variables may be a handy tool in the analysis of the more problematic cases, like the eccentric, octupolar Lidov-Kozai problem. In the latter, the ‘orbital flip’ phenomenon occurs: changing the direction of motion with the passage through an equatorial rectilinear orbit phase (Lithwick and Naoz, 2011). Previous attempts to discuss this phenomenon in terms of the action-angle variables (e.g Sidorenko, 2018) faced the problems which may possibly be resolved with the newly presented parametrization.

Some of the readers might be sceptical about the unnecessary duplication of the phase space resulting from the LKS transformation $\zeta$. Indeed, Fig. 3 covers the whole phase space of in terms of the argument of pericentre $\omega_{\mathrm{o}}$, although it has been clipped to the half range of $\lambda$. This feature can be trivially removed by means of a symplectic transformation $(\lambda,\Lambda)\rightarrow(2\lambda,\Lambda/2)$, and similarly for other conjugate pairs. We have not made this move in the present work for the sake of retaining the fundamental, angle-halving property of both the Levi-Civita and the Kustaanheimo-Stiefel transformations. Avoiding factor 2 in the arguments of sines and cosines in equations (69) and (72), we would introduce the factor $\frac{1}{2}$ in the expressions for $\mathsf{v}$ and $\mathsf{V}$. Let us mention that the restriction of the LKS transformation to $(\mathsf{v},\mathsf{V})\rightarrow(l,g,h,\gamma,L,G,H,\Gamma)$ can be useful also in the studies of perturbed, four degrees of freedom oscillators, not necessarily resulting from the KS transformation (e.g. Crespo et al, 2015; van der Meer et al, 2016). In that case, unwanted spurious singularities may arise in course of the Birkhoff normalization, when the multiple of angle does not properly match the power of action.

Having based the LKS variables upon the KS3 variant of the KS transformation, we do not exclude a possibility of performing a similar construction within the KS1 framework. But then the $G$ and $\Lambda$ variables will be the projections of the angular momentum and the Laplace-Runge-Lenz vectors on the $x_{1}$ axis. With such a choice, the Lidov-Kozai Hamiltonian (104) would depend on both $g$ and $\lambda$, with the rotation symmetry hidden deeply in some complicated function of all variables, instead of the obvious $G=\mbox{const}$.

**Compliance with ethical standards

** 

Conflict of interest The authors S. Breiter and K. Langner declare that they have no conflict of interest.