KS variables in rotating reference frame. Application to cometary dynamics

The problem of Keplerian motion in the uniformly rotating reference frame has been solved in terms of the Kustaanheimo-Stiefel (KS) variables using canonical formalism. No recourse to the Cartesian variables or orbital elements has been required. The form of solution is well suited for the application as a part of symplectic integrator. The results show that the motion is actually the composition of four independent harmonic oscillations and of the rotation in two specific coordinate planes and their conjugate momenta planes. As an example of application, we use the KS symplectic integrator to study the motion of comet C/1997 J2 (Meunier-Dupuoy) under the action of Galactic tides. The comet is found to follow an orbit in commensurability with the Sun motion around the Galactic centre, but the perturbations are not qualified as a resonance.


Introduction
transformation, transforms the three-dimensional Kepler problem into a four-dimensional isotropic harmonic oscillator. As a consequence, it attains both linearization and regularization of motion (Kustaanheimo and Stiefel 1965;Stiefel and Scheifele 1971 formation involves a number of subtleties related with dimension raising (Stiefel and Scheifele 1971;Kurcheeva 1977;Deprit et al. 1994), its very existence makes KS variables attractive for the application in numerical integration using symplectic methods. Thanks to the regularization, the fundamental problem of using constant integration step even for highly eccentric orbits is considerably reduced in KS variables (Mikkola 1997;Breiter 1998).
The standard form of the KS transformation is given in the reference frame with fixed direction of axes. But a reference frame with the origin in the primary body and uniformly rotating around one of its axes is a common feature of various classical problems in celestial mechanics. Most often it helps to remove explicit time dependence from potential. In artificial satellite theory, if the tesseral part of potential is included, the axes may co-rotate with the planet (e.g. Segerman and Coffey 2000). In the circular restricted three-body problem (and its variants, like the Hill problem) one of the axes passes through two primaries (e.g. Szebehely 1967). Actually, in the latter case, it is more common to set the origin at the centre of mass of the primaries, but there are notable exceptions, like Kurcheeva (1977). In the dynamics of the Oort Cloud comets, it is also advantageous to handle the Galactic tide perturbations in the heliocentric frame co-rotating with the centre of the Galaxy (e.g. Breiter et al. 2008).
The last example is noteworthy. In this particular problem, Breiter et al. (2007) worked out a symplectic integrator in KS variables using a partition of Hamiltonian into Keplerian part and perturbation. But the former part requires the knowledge of explicit solution to construct a propagator of coordinates, momenta, and optionally of their variations. Lacking a satisfactory solution in rotating reference frame, the authors stayed in the fixed coordinate system, accepting the time dependence of perturbation. The present work aims at providing the formulae necessary to construct a Keplerian motion propagator for the KS variables and their variations in rotating frame.
Actually, the problem was already discussed under the name of 'generating solution' for the circular, restricted problem of three bodies, meaning the limiting case when the ratio of masses for the two primaries tends to zero. Ahmad and Huda (1986) provided such solution using directly the KS variables, but they made an inappropriate and too particular choice of integration constants values, criticized by Hassan (1994). Hassan (1994Hassan ( , 2001 published his own solution, but it was expressed in terms of Keplerian orbit elements, requiring an additional there-and-back trip to Cartesian coordinates and momenta. In this paper we present a direct and explicit solution, with no resort to Cartesian variables or orbital elements. It can be easily incorporated into a canonical splitting method of numerical integration. As an illustrative example we show the application of such integrator to the heliocentric motion of a comet perturbed by the tidal potential of the Galaxy. We find that the test body-comet C/1997 J2 (Meunier-Dupouy), moves on an interesting, but non-resonant orbit.

Initial problem statement
Let us consider perturbed Keplerian motion of a material point with mass m p around a central body with mass m c . The origin of reference frame is placed in the centre of mass of the primary and the third axis,ẑ has a fixed direction. The remaining two axes,x andŷ rotate aroundẑ with a constant angular rate Ω. The motion can be described using Cartesian coordinates r = [x, y, z] T and conjugate momenta R = [X, Y, Z] T , which obey canonical equationṡ where the second argument of canonical Poisson brackets is the Hamiltonian function Keplerian motion results from where the third component G 3 of angular momentum vector G is the gravitational parameter of relative two-body problem, involving the Gaussian constant k, is defined as and r = r = x 2 + y 2 + z 2 .
H 1 is a perturbing potential-typically a function of coordinates r and time t (the independent variable of the differential system (1)).
Let us recall that although we use momenta and energy divided by the mass m p , the momentum is not equal to velocity, save for its third component. Indeed, Eq. (1) and the form of H implẏ In other words, the momentum vector R is the velocity vector measured in the fixed frame, whose components are resolved along instantaneous axes of the rotating frame at given epoch.

KS transformation
Within canonical formalism, the Kustaanheimo-Stiefel transformation is performed in extended phase space as a combination of point transformation of coordinates, its canonical extension for the momenta, and Sundman transformation of the independent variable. We briefly recall its basic elements, with additional observation that whenever r or R are mentioned, they stand for coordinates and momenta in the rotating frame. Using the convention of Deprit et al. (1994), who introduced an arbitrary length parameter α for the consistency of units, we define coordinates v j , with j = 0, . . . , 4, through the inverse transformation ⎡ For the record, we should warn readers that a more widespread convention of Stiefel and Scheifele (1971) uses similar u 1 = v 1 , u 2 = v 2 , u 3 = v 3 , but their u 4 = −v 0 , and α = 1. Canonical extension of new coordinates v results in their conjugate momenta V implicitly given by ⎡ The direct transformation will not be used in present studyit can be found in Deprit et al. (1994). A byproduct of the extension is the important bilinear invariant Its Poisson bracket with any function of r, R, and t (expressed in terms of KS variables) is zero, which allows to drop the occurrence of J in all formulae, provided we discuss only the motion generated by the KS image of a Hamiltonian H(r, R, t).
Finally, we introduce a new independent variable-the Sundman time τ , a close relative of the eccentric anomaly, but having the units of time, This step requires the introduction of a new pair of variables.
A new coordinate type variable v * = t + const. (let us call it formal time) will replace any possible presence of time in the Hamiltonian. The value of its canonically conjugate momentum V * , should obey at any epoch (Szebehely 1967, Sect. 6.6). In further discussion, we will not append the fifth component to formal vectors enumerating v * and V * separately. Thus the Poisson bracket, for arbitrary functions f and g, can be specified in new variables as KS transformation has two important properties: the distance r is a quadratic function of v j and where a J 2 term inside the bracket has been dropped. Both identities facilitate the derivation of the new Hamiltonian After elementary substitutions, we find K = K 0 + K 1 , with two basic components. The term K 0 , responsible for the Keplerian motion in rotating frame is and the perturbation K 1 is Observing formal similarity of K 0 to the harmonic oscillator, we should emphasize, that the frequency ω (real positive in elliptic motion) is a function of KS variables and of V * , namely where G 3 , resulting from the substitution of (8) and (9) into (4), followed by rejecting zJ /(2r) = 0, is The constant term −4μα −1 in Eq. (18) may seem redundant, because it has no influence on equations of motion. But since the initial value of V * is determined from condition K = 0, rejecting this constant would lead to a false value of the oscillator frequency ω, defined by Eq. (20). In general, the frequency ω may vary during the motion. Constant value of ω requires the original Hamiltonian H 1 which is both explicitly time-independent (to conserve V * ), and invariant with respect to rotations around z (to conserve G 3 ). Further discussion will focus on the case of Keplerian motion, where the conditions become trivially satisfied. Yet, before closing the section, we have to add an important remark.
Considering a pure, stand-alone Keplerian problem, we obviously have which can be used as a strict identity to simplify expressions. But if the Keplerian propagator serves as a part of a splitting method, like a symplectic integrator with alternating application of motions generated by K 0 and K 1 , then, instead of (22), one should substitute where E is constant or variable, depending on the properties of K 1 . But even if E varies in the general problem, its value during Keplerian propagation in a splitting algorithm remains frozen, hence we will not differentiate it with respect to τ in the next section (see Breiter 1998;Breiter et al. 2007).

Equations of motion
Using the Poisson bracket (14), we find the equations of motion generated by K 0 . Let us begin with a trivial meaning that V * = const., and an elementary which reconstructs the definition (11). In next equations we will immediately substitute r, using Eq. (15). Thus, for the KS coordinates we find and for the momenta Freezing rotation (Ω = 0) we recover the usual set of harmonic oscillator equations. But in general case, Eqs. (26) and (27) form a coupled system, because with a skew-symmetric matrix The resulting system is not linear, because r is a quadratic function of v. Fortunately, the equations of motion can be reduced to a linear system with periodic coefficients and then neatly separated to obtain an explicit solution.

Explicit solution
The differential system of Eqs. (25), (30), (31) will be solved as the initial value problem. For convenience, we assume the initial epoch of the Sundman time τ = 0. At this epoch, the time-like variable v * (0) = u * , KS coordinates are v(0) = u, and their conjugate momenta are V(0) = U. Knowing that V * is constant, we will not introduce a special symbol for its initial value.

Radius
As the first step towards solving the system, we should find r as the explicit function of τ . The success is almost guaranteed, because the solution for r(τ ) in the fixed frame is well known (up to a factor and phase, τ is the eccentric anomaly in the Kepler problem), and the distance is invariant with respect to rotation. Let us begin with the first derivative of r = α −1 v 2 , where we recover the same formula as in the fixed frame: The second differentiation involves a subtle point, signaled in the last paragraph of Sect. 2.2. Initially, we obtain Then, from Eq. (18), we can substitute with frequency 2ω, and E which either zero (pure Kepler), or follows from Eq. (23).
With the initial conditions we find the solution where is simply the half increment of eccentric anomaly. Apart from the different definition of ω, Eq. (36) agrees with the usual expressions known from the fixed frame.

Formal time v *
Once the function r(τ ) has been determined, we can find the formal time v * , which is equal to physical time t up to an arbitrary phase. Equation (11) can be solved by quadratures, with Substituting (36) we find the equivalent of Kepler's equation.
Differentiating radial solution (36) and making use of (32) and (35), we can relegate periodic terms in , obtaining although this form is useful only if the computation of the frame rotation invariant product v(τ ) · V(τ ) is done independently on v(τ ) and V(τ ).

Coordinates and momenta
Once we start considering r as a known function of time, we can observe that Eqs. (30) and (31) can be split into two separate, and formally identical subsystems. Selecting we find their equations of motion The matrix M is amenable to a standard reduction to diagonal form. Solving its characteristic equation we find four distinct imaginary eigenvalues Eigenvectors s i , associated with these values, are linearly independent: the matrix S = [s 1 |s 2 |s 3 |s 4 ] built from eigenvectors has the form so det S = ω 2 = 0. Thus, if Λ = diag(λ 1 , λ 2 , λ 3 , λ 4 ) is a diagonal matrix built from the eigenvalues, we have where Remarkably, all eigenvectors are constant in the Keplerian problem. Thus, thanks to we can perform a linear transformation of variables bringing the system (43) into diagonalized form. For g we find so, with a similar treatment of h, In these circumstances, each equation can be solved by quadratures, leading to and similarly for h j . According to Eqs. (46) and (25), the quadrature results in with the choice of signs for depending on j -similar as in the definition of λ j (46). Thus, introducing and matrix we can write the solution for principal modes With g(0) = S −1 p(0), and h(0) = S −1 q(0), where we can transform back to These can be rearranged, leading to the explicitly real matrix form Two matrices present in the final expression are: matrix W, performing the simultaneous retrograde rotations of two parts of a vector by some angle β and matrix A combining permutation and scaling of vector elements with the inverse For a specific pair of coordinate and momentum in Eqs. (62,63), we can write taking plus for j = 0, 1 and minus for j = 2, 3. Suppressing the rotation of reference frame (φ = 0), we obtain the well known expressions without the coupling of different degrees of freedom (Breiter et al. 2007).

Comments
The main result of Sect. 3 is the set of equations (40) and (67) which, given the initial values of KS coordinates, momenta and physical time, allow to compute their values at any epoch of the Sundman time τ . In other words, they constitute a Keplerian propagator for the KS variables in a rotating reference frame. The propagator can then be used, in particular, to construct a symplectic integrator, provided some appropriate perturbation part is added. Although Eq. (67) are best for practical computations, their matrix form (62) and (63) is more pleasing from the theoretical point of view. First, let us observe the way the symmetry is broken. In the fixed frame, the Keplerian motion in KS variables is totally isotropic and may be considered as four independent oscillations, i.e. clockwise rotations on the phase plane (v j , ω −1 V j ), with the momenta scaled by the frequency. One may also pick up six possible 'elliptic oscillators' with two degrees of freedom (v j , v l , V j , V l ), and each of them will admit an angular momentum integral G jl = v j V l − v l V j (Deprit 1991). Remarkably, each component of the Cartesian angular momentum G = r × R, when expressed in terms of the KS variables, is a linear combination of two different G jl .
After adding the term (−4r/α)ΩG 3 to the Keplerian Hamiltonian, we observe that G 3 is the only component of G that remains constant. Recalling Eq. (21), we can observe that G 3 = (G 03 + G 12 )/2. And, since {G 03 , G 12 } = 0, both G 03 and G 12 remain constant in the rotating frame. This fact explains the partition of (v, V) into p and q, which can be handled separately.
Subdividing Eqs. (62) and (63) we find the evolution of variables as the composition of two rotations. First, p and q are obtained according to the motion of four independent oscillators on the phase planes (v j , V j ). This requires first permutation and scaling by matrix A, to obtain independent (v j , ω −1 V j ) pairs. Then, each pair follows the clockwise (i.e. negative) rotation on the phase plane, according to the well known properties of harmonic oscillator. After the rotation, we recover the original arrangement (v j , v l , V j , V l ) by the application of A −1 . The resulting p and q are subjected to the final rotation that acts separately in the coordinate plane (v j , v l ) and in the momentum plane (V j , V l ). Since this is a passive rotation, caused by the motion of the reference frame, it is again clockwise (negative). The result being a composition of the fixed frame Keplerian motion and of the reference frame rotation, might have been guessed from the properties of the motion in Cartesian coordinates, where the two parts of the Hamiltonian (3) are in involution, i.e.
But, unfortunately, it is not the case in the KS variables, where Thus we do not exclude the possibility of a short proof that the solution must have the form (62) and (63), but-for the noncommuting Hamiltonian parts-it will not be straightforward.
Most of the comments about real and imaginary values, scattered through the text, refer to the elliptic motion (including degenerate radial orbits). But if the motion is hyperbolic, the expressions can be adapted without any difficulty. Enough to replace each ω by i|ω| and use elementary properties of sin iy = −i sinh y, and cos iy = cosh y, as did Breiter et al. (2007). In that case, the motion is a composition of hyperbolic rotation in phase plane and circular rotation in coordinate space and momenta space.

Galactic tides integrator
In order to validate the utility of presented solution, we apply it to the problem of heliocentric motion of a distant comet perturbed by Galactic tides. A symplectic integrator designed for this problem by Breiter et al. (2007) worked in a heliocentric reference frame with fixed direction of axes. But the natural framework for this problem is the system with axis Ox pointing towards the Galactic centre, Oz aligned with the angular momentum of the Sun in its motion around the Galactic centre (we assume the planar motion in the Galactic plane) and Oy completing the right-handed triad. Moreover, it is better to shift the origin to the centre of mass of the Solar system, which practically amounts to using the gravitational parameter where M is the total mass of the Sun and planets (JPL DE405 ephemeris value). The system is then described by Hamiltonian (2) with H 0 given by Eq. (3), and perturbation Similarly to Breiter et al. (2007), we use the constants G 2 = 7.0706 × 10 −16 yr −2 , G 3 = 5.6530 × 10 −15 yr −2 , and the reference frame rotation rate Ω = − √ G 2 . In terms of KS variables, the problem is defined by the Hamiltonian K, being the sum of K 0 , as defined by Eq. (18), and K 1 , which according to Eq. (19) is Technically, retaining x, y, z, and r leads to more compact expressions, although these symbols are to be understood as the functions v, defined by Eqs. (8) and (15). A symplectic integrator based on the splitting method is effectively an alternating application of two maps Φ 0,h and Φ 1,h , each of them representing the motion generated by the Hamiltonians K 0 and K 1 , respectively, with appropriately chosen time intervals h and h . In this work we use a simple, second order 'leapfrog' composition, where a single integrator step over the interval h of Sundman time amounts to where ε is, roughly speaking, the ratio of perturbation H 1 to the leading part H 0 . The map Φ 0 has been described in details in Sect. 3. The second map is elementary, because the absence of V in K 1 implies so, with constant coordinates v, equations for momenta The gradient of K 1 with respect to v can be compactly written as where Both formal time v * and its conjugate momentum V * are constant under the action of Φ 1 , because K 1 is independent on each of the variables.

The motion of C/1997 J2 (Meunier-Dupouy)
The symplectic integrator described in previous section has been applied to the comet C/1997 J2 (Meunier-Dupouy). The comet has well determined orbital elements and was discussed by Królikowska and Dybczyński (2010), who noticed it belongs to a rare class of objects having previous perihelion distance lower than the current one. The barycentric osculating elements given in Table 1 were taken from the WikiComet web page of P. Dybczyński 1 and transformed to the Galactic frame. The excessive number of digits serves only to allow independent repetition of the present calculations; actual uncertainties can be found at the WikiComet link. It is worth noting that the heliocentric orbit is hyperbolic, but in the barycentric frame we obtain a highly elliptic osculating orbit-a great challenge for a fixed step numerical integration without regularization. According to the results of Dybczyński and Królikowska (2011), neglecting the contribution of the Galactic centre (i.e. setting G 2 = 0 in Eq. (71)) does not seriously influence the computed value of the previous perihelion passage distance of C/1997 J2, leading to the difference of 10 %much less than for the remaining three comets from their Table 11. But this fact should not be taken for evidence that the Galactic centre has negligible effect on the overall motion of C/1997 J2-it is quite the opposite (P. Dybczyński, private communication).
Starting from the initial conditions computed from the elements given in Table 1, we have integrated the past motion of C/1997 J2 perturbed by the Galactic tides alone. The orbit was traced back over 3.78 Gy which amounts to 1128 of nominal Keplerian periods P ≈ 3.35 My, with a fixed time step of Sundman time τ equivalent to 0.04 of P . Actually, the integration span was chosen to cover 16 periods of the reference frame rotation P Ω = 2πΩ −1 ≈ 236.3 My.
Thanks to the use of rotating frame, we were able to check the error in the conservation of time-independent Hamiltonian K = K 0 + K 1 , evaluated at each step from the current values KS variables v and V. Its relative value (divided by V * , because K itself should be equal to 0) oscillated with an amplitude not exceeding 2 × 10 −8 , without any noticeable systematic trend.
Skipping the question of minor, short period perturbations in semi-axis, modulated by T Ω , and less interesting evolution of the mean anomaly, let us proceed to the remaining elements. Figure 1 compares the results obtained with the complete model (solid line) with those of the Galactic disc alone, i.e. obtained with G 2 = 0 (dashed line). Concerning the eccentricity e, Galactic disc perturbations generate periodic perturbations resembling a sine wave with the dominating period very close to 4T Ω (4 cycles over the time span of 16T Ω ). This commensurability of periods has no dynamical consequences in the axially symmetric disc problem, but in the presence of the Galactic centre it may trigger a resonance. In the complete model we see the presence of equally strong periodic terms, with the period of 1 2 T Ω , superimposed on the disc contribution. Inspecting the inclination I we find a similar, although more significant pattern: the disc contribution with period 4T Ω serves as an envelope for the Galactic centre terms with the period of 1 2 T Ω , but this time the latter allow transitions between prograde to retrograde motion. The argument of pericentre g and the longitude of the ascending node h (transformed to the fixed axes) behave quite similarly, remaining almost constant for most of the time, yet simultaneously jumping by 180 • every 4T Ω . The effect of Galactic centre is barely visible in both the variables.
Another source of information is the evolution of the Laplace vector e directed to the pericentre and having the length equal to the eccentricity. In Fig. 2 we show its plot in the rotating frame. Only the results for the complete model are displayed, but the picture resulting from the Galactic disc alone looks practically similar. The tip of e traces a helix, oscillating up and down. The two complete cycles with period 8T ω almost coincide; the beginning and the end of the track  Do the results of simulation support the conjecture that the motion of C/1997 J2 is significantly affected the resonance with Galactic centre? In our opinion, the answer is negative. Lacking an analytical theory that might show the absence of a relevant periodic term in the normalized Hamiltonian, we have resorted to a simple test, which is possible thanks to the property, that Galactic disk perturbations period is proportional to the mean motion. Thus, we have integrated the motion of a hypothetic comet with all osculating elements of C/1997 J2, save for the semi-axis, increased to a = 25000 au. Now, in order to cover 16T Ω , we followed past 957 of nominal Keplerian periods, in agreement with the ratio (25/22.4) 3 2 ≈ 1.18. A resonant effect should not survive such a detuning and the amplitudes should change considerably. But, as we show in Fig. 3, nothing serious has happened. In spite of different 'envelope' period, the Galactic center contribution with period 1 2 T Ω has actually grown by the factor comparable to the ratio of semi-axes. Fig. 3 Osculating eccentricity and inclination for C/1997 J2 (grey) and its twin with a = 2500 au (black) This argument alone suffices to rule out the resonant nature of the Galactic centre perturbations on the Meunier-Dupouy comet. Their relative strength is mainly the effect of rather weak perturbations due to Galactic disc. This, in turn, is related to the low latitudes of the line of apsides (see Fig. 2); together with high eccentricity it restricts the deviation of the comet from the Galactic plane, whereas the action of the disc is proportional to z.

Conclusions
Providing an explicit solution to the problem of Keplerian motion in rotating frame in terms of the Kuustanheimo-Stiefel variables was the main objective of the presented work. This is more than just an exercise in elementary Hamiltonian mechanics: the solution is an important building block for various practical applications. As an example, we have used to create a symplectic integrator for the action of Galaxy on cometary orbits. Using it, we were able to simulate the motion of comet C/1997 J2 (Meunier-Dupuoy) over almost 4 Gy in a fraction of second. The integration time was comparable to the time spent on converting the results to the osculating elements form. The results show the 1:4 commensurability between the perturbations due to the Galactic disc and the rotation of the Galactic centre with respect to the Sun. But this commensurability does not results in a significant resonance.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.