A&A 415, 1187-1199 (2004)
DOI: 10.1051/0004-6361:20034058

Gauge freedom in the N-body problem of celestial mechanics

M. Efroimsky¹ - P. Goldreich²

1 - US Naval Observatory, Washington, DC 20392, USA
2 - IAS, Princeton NJ 08540 & CalTech, Pasadena, CA 91125, USA

Received 8 July 2003 / Accepted 22 October 2003

Abstract
The goal of this paper is to demonstrate how the internal symmetry of the N-body celestial-mechanics problem can be exploited in orbit calculation.

We start with summarising research reported in (Efroimsky 2002, 2003; Newman & Efroimsky 2003; Efroimsky & Goldreich 2003) and develop its application to planetary equations in non-inertial frames. This class of problems is treated by the variation-of-constants method. As explained in the previous publications, whenever a standard system of six planetary equations (in the Lagrange, Delaunay, or other form) is employed for N objects, the trajectory resides on a 9(N-1)-dimensional submanifold of the 12(N-1)-dimensional space spanned by the orbital elements and their time derivatives. The freedom in choosing this submanifold reveals an internal symmetry inherent in the description of the trajectory by orbital elements. This freedom is analogous to the gauge invariance of electrodynamics. In traditional derivations of the planetary equations this freedom is removed by hand through the introduction of the Lagrange constraint, either explicitly (in the variation-of-constants method) or implicitly (in the Hamilton-Jacobi approach). This constraint imposes the condition (called "osculation condition'') that both the instantaneous position and velocity be fit by a Keplerian ellipse (or hyperbola), i.e., that the instantaneous Keplerian ellipse (or hyperbola) be tangential to the trajectory. Imposition of any supplementary constraint different from that of Lagrange (but compatible with the equations of motion) would alter the mathematical form of the planetary equations without affecting the physical trajectory.

However, for coordinate-dependent perturbations, any gauge different from that of Lagrange makes the Delaunay system non-canonical. Still, it turns out that in a more general case of disturbances dependent also upon velocities, there exists a "generalised Lagrange gauge'', i.e., a constraint under which the Delaunay system is canonical (and the orbital elements are osculating in the phase space). This gauge reduces to the regular Lagrange gauge for perturbations that are velocity-independent.

Finally, we provide a practical example illustrating how the gauge formalism considerably simplifies the calculation of satellite motion about an oblate precessing planet.

Key words: celestial mechanics - reference systems - solar system: general - methods: N-body simulations - methods: analytical - methods: numerical

   
1 Introduction

1.1 Prefatory notes

On the 6th of November 1766 young geometer Giuseppe Lodovico Lagrangia, invited from Turin at d'Alembert's recommendation by King Friedrich the Second, succeeded Euler as the Director of Mathematics at the Berlin Academy. Lagrange held the position for 20 years, and this fruitful period of his life was marked by an avalanche of excellent results, and by three honourable prizes received by him from the Académie des Sciences of Paris. All three prizes (one of which he shared with Euler) were awarded to Lagrange for his contributions to celestial mechanics. Among these contributions was a method introduced by Newton and Euler and furthered by Lagrange for his studies of planet-perturbed cometary orbits and only later applied to planetary motion (Lagrange 1778, 1783, 1808, 1809, 1810). The method was based on an elegant mathematical device, the variation of constants emerging in solutions of differential equations. This approach was pioneered by Newton in his unpublished Portsmouth papers and described very succintly in Cor. 3 ans 4 of Prop. 17 in Book I of his "Principia". The first attempts of practical implementation of this tool were presented in a paper on Jupiter's and Saturn's mutual disturbances, submitted by Euler to a competition held by the French Academy of Sciences (Euler 1748), and in the treatise on the lunar motion, published by Euler in 1753 in St. Petersburg (Euler 1753). However, it was Lagrange who revealed the full power of this approach.

Below we shall demonstrate that the equations for the instantaneous orbital elements possess a hidden symmetry similar to the gauge symmetry of electrodynamics. Derivation of

Figure 1: Two coplanar ellipses, having one common focus, are assumed to rotate about this focus, always remaining within their plane. Let a planet be located at one of the points of the ellipses' intersection, P, and assume that the rotation of the ellipses is such that the planet is always at the instantaneous point of their intersection. At some instant of time, the rotation of one ellipse will be faster than that of the other. On these grounds one may state that the planet is rapidly moving along the slower rotating ellipse. On the other hand, though, it is also true that the planet is slowly moving along the faster rotating ellipse. Both viewpoints are equally valid, because one can divide, in an infinite number of ways, the actual motion of the planet into a motion along some ellipse and a simultaneous evolution of this ellipse. The Lagrange constraint (7) singles out, of all the sequences of evolving ellipses, that unique ellipse sequence which is always tangential to the physical velocity of the planet.

Open with DEXTER

the Lagrange system involves a mathematical operation equivalent to the choice of a specific gauge. As a result, trajectories get constrained to some 9-dimensional submanifold in the 12-dimensional space constituted by the orbital elements and their time derivatives, as demonstrated in (Efroimsky 2002, 2003; Newman & Efroimsky 2003; Efroimsky & Goldreich 2003). However, the choice of this submanifold is essentially ambiguous, and this ambiguity gives birth to an internal symmetry. The symmetry is absent in the two-body case, but comes into being in the N-body setting ($N \ge 3$) where each body follows an ellipse of varying shape whose time evolution contains an inherent ambiguity.

For a simple illustration of this point imagine two coplanar ellipses sharing one focus (Fig. 1). Suppose they rotate at different rates in their common plane. Let a planet be located at one of the intersection points of these ellipses. The values of the elliptic elements needed to describe its trajectory would depend upon which ellipse was chosen to parameterise the orbit. Either set would be equally legitimate and would faithfully describe the physical trajectory. Thus we see that there exists an infinite number of ways of dividing the actual planet's movement into motion along its orbit and the simultaneous evolution of the orbit. Although the physical trajectory is unique, its description (i.e., its parametrisation in terms of Kepler's elements) is not. A map between two different (though physically equivalent) sets of orbital elements is a symmetry transformation (a gauge transformation, in physicists' jargon).

Lagrange never dwelled upon that point. However, in his treatment (based on direct application of the method of variation of constants) he passingly introduced a convenient mathematical condition which removed the said ambiguity. This condition and possible alternatives to it will be the topics of Sects. 1-3 of this paper.

In 1834-1835 Hamilton put forward his theory of canonical transformations. Several years later this approach was furthered by Jacobi who brought Hamilton's technique into astronomy and, thereby, worked out a new method of deriving the planetary equations (Jacobi 1866), a method that was soon accepted as standard. Though the mathematical content of the Hamilton-Jacobi theory is impeccably correct, its application to astronomy contains a long overlooked aspect that needs attention. This aspect is: where is the Lagrange constraint tacitly imposed, and what happens if we impose a different constraint? This issue will be addressed in Sect. 4 of our paper.

The main line of reasoning and the principal results presented in this paper are the following. In a concise introduction, to be presented in the next subsection, we derive Eq. (16) which is the most general form of the gauge-invariant perturbation equation of celestial mechanics, written in terms of a disturbing force. Then we transform it into (25), which is the most general gauge-invariant perturbation equation expressed through a disturbance of the Lagrangian. The next step is to show that, in the case of velocity-dependent perturbations, the equations for the Delaunay elements are, generally, not canonical. However, they become canonical in a special gauge (which is, generally, different from the customary Lagrange constraint).

Our other goal is to explore how the freedom of gauge choice is supplemented by the freedom of choice of a coordinate frame in which to implement the method of variation of constants. This investigation, carried out in Sect. 3, provides a physical example of a nontrivial gauge being instrumental in simplifying orbit calculations. In this way we demonstrate that gauge freedom in celestial dynamics is not only a cultural acquisition, but can be exploited to improve computations.

   
1.2 Osculating elements vs. orbital elements

We start in the spirit of Lagrange. Before addressing the N-particle case, Lagrange referred to the reduced two-body problem,

 

$${{\ddot {\vec {r}}}} + \frac{\mu}{r^2}\; \frac{{{\vec r}}}{r} = 0,$$

$$\mbox{{$\vec r$ }}\;\equiv\;\mbox{{$\vec r$ }}_{{\rm planet}} - \... ...m sun}},\;\;\;\;\;\;\;\;\;\;\;\;\;\mu\;\equiv\;G(m_{\rm planet} + m_{\rm sun}),$$ (1)

whose generic solution, a Keplerian ellipse or a hyperbola, can be expressed, in some fixed Cartesian frame, as

$$% \mbox{{$\vec r$ }}= {\mbox{{$\vec f$ }}} \left(C_1, ..., C_... ...{$\vec r$ }}}} = {\vec g}\left(C_1,\; ...,\; C_6, ~t \right),$$ (2)

where

$$% {{{\vec{g}}}} \equiv \left(\frac{\partial { {\mbox{{$\vec f$ }}}}}{ \partial t}\right)_{C={\rm {const}}}\cdot$$ (3)

Since the problem (1) constitutes by three second-order differential equations, its general solution naturally contains six adjustable constants C_i. At this point it is irrelevant which particular set of the adjustable parameters is employed. (It may be, for example, a set of Lagrange or Delaunay orbital elements or, alternatively, a set of initial values of the coordinates and velocities.)

Following Lagrange (1808-1810), we take ${{\mbox{{$\vec f$ }}}}$ as an ansatz for a solution to the N-particle problem, the disturbing force acting at a particle being denoted by  $\Delta{{\vec F}}$:

$$% {\ddot{\vec r}} + \frac{\mu}{r^2}\;\frac{{\vec r }}{r} = {\Delta {\vec F}}.$$ (4)

Now the "constants" become time dependent:

$$% {\vec r } = {\mbox{{$\vec f$ }}} \left(C_1(t), ..., C_6(t), ~t \right),$$ (5)

whence the velocity

$$% \frac{{\rm d} \mbox{{$\vec r$ }}}{{\rm d}t} = \frac{ \parti... ...x{{$\vec f$ }}}}{\partial C_i}\;\frac{{\rm d} C_i}{{\rm d} t},$$ (6)

acquires an additional term, $\sum ({\partial{\mbox{{$\vec f$ }}}}/{\partial C_i})({{\rm d}C_i}/{{\rm d}t})$.

Substitution of (5) and (6) into the perturbed equation of motion (4) leads to three independent scalar second-order differential equations which contain one independent parameter, time, and six functions C_i( t). These are to be found from the three scalar equations comprised by the vector Eq. (4), and this makes the essence of the variation-of-constants method. However, the latter task cannot be accomplished in a unique way because the number of variables exceeds, by three, the number of equations. Hence, though the physical trajectory (comprised by the locus of points in the Cartesian frame and by the values of velocity at each of these points) is unique, its parametrisation through the orbital elements is ambiguous. This circumstance was appreciated by Lagrange, who amended the system, by hand, with three independent conditions,

$$% \sum_i \;\frac{\partial {\mbox{{$\vec f$ }}}}{\partial C_i}\;\frac{{\rm d} C_i}{{\rm d} t} = 0,$$ (7)

and thus made it solvable. His choice of constraints was motivated by both physical considerations and mathematical expediency. Physically, the set of functions  $\left(C_1(t), ..., C_6(t)\right)~$ can be interpreted as parameters of an instantaneous ellipse (in a bound-orbit case) or an instantaneous hyperbola (in a flyby situation). Lagrange found it natural to make the instantaneous orbital elements C_i such that, at each moment of time, this ellipse (hyperbola) coincides with the unperturbed (two-body) orbit that the body would follow if the disturbances were to cease instantaneously. In other words, these ellipses (hyperbolae) are tangential to the physical trajectory. This is the well known condition of osculation, which is mathematically implemented by (7). The appropriate orbital elements are called osculating elements. This way, Lagrange restricted his use of the orbital elements to elements that osculate in the reference frame wherein ansatz (5) is employed. Lagrange never explored alternative options; he simply imposed (7) and used it to derive his famous system of equations for orbital elements.

Such a restriction, though physically motivated, is completely arbitrary from the mathematical point of view. While the imposition of (7) considerably simplifies the subsequent calculations it in no way influences the shape of the physical trajectory and the rate of motion along it. As emphasised in Efroimsky (2002, 2003) and Newman & Efroimsky (2003), a choice of any other supplementary constraint

$$% \sum_i \;\frac{\partial {\mbox{{$\vec f$ }}}}{\partial C_i}... ...} C_i}{{\rm d} t} = {\bf {\vec \Phi}}(C_{1},~...,~C_{6},~t),$$ (8)

$\mbox{{\boldmath$\vec \Phi$ }}$ being an arbitrary function of time and parameters C_i, would lead to the same physical orbit and the same velocities. Substitution of the Lagrange constraint (7) by its generalisation (8) would not influence the motion of the body, but would alter its mathematical description (i.e., would entail different solutions for the orbital elements). Such invariance of a physical picture under a change of parametrisation is called gauge freedom or gauge symmetry. It parallels a similar phenomenon well known in electrodynamics and, therefore, has similar consequences. On the one hand, a "good" choice of gauge often simplifies solution of the equations of motion. (In Sect. 3 we provide a specific application to motion in an accelerated coordinate system.) On the other hand, one should expect a displacement in the gauge function  $\mbox{{\boldmath$\vec \Phi$ }}$, owing to the accumulated error in the constants.

A simple case of the same motion being described by two different families of instantaneous ellipses is presented in Fig. 1. Two coplanar ellipses have a common focus and are rotating about it, in the same plane but at different rates. Assume that one ellipse is rotating rapidly and another slowly. Then it will be legitimate to state that the point of their intersection, P, is moving rapidly along a slowly rotating ellipse. At the same time, it will be right to say that it is moving slowly along a swiftly rotating ellipse.

Derivation of the conventional Lagrange and Delaunay planetary equations by the variation-of-constants method incorporates the Lagrange constraint (Brouwer & Clemence 1961). Both systems of equations get altered under a different gauge choice, as we now show. The essence of a derivation suitable for a general choice of gauge starts with (6) from which the formula for the acceleration follows:

 

$$% \frac{{\rm d}^2 \mbox{{$\vec r$ }}}{{\rm d}t^2} = \frac{\partia... ... C_i}\;\frac{{\rm d} C_i}{{\rm d} t} + \frac{{\rm d}{\bf\vec \Phi}}{{\rm d}t} = \frac{\partial^2 {{\mbox{{$\vec f$ }}}}}{\partial^2 t} + \sum_i\;... ... }\;\frac{{\rm d} C_i}{{\rm d} t} + \frac{{\rm d}{\bf\vec \Phi}}{{\rm d}t}\cdot$$ (9)

Together with the equation of motion (4), it yields:

 

$$% \frac{\partial^2 {{\mbox{{$\vec f$ }}}} }{\partial t^2} + \frac... ...d} C_i}{{\rm d} t} + \frac{{\rm d}{\bf\vec \Phi}}{{\rm d}t} = \Delta{{\vec F}},$$ (10)

$$r \equiv \vert\mbox{{$\vec r$ }}\vert = \vert{{\mbox{{$\vec f$ }}}}\vert.$$

The vector function ${{\mbox{{$\vec f$ }}}}$ was from the beginning introduced as a Keplerian solution to the two-body problem; hence it must obey the unperturbed Eq. (1). On these grounds the above formula simplifies to:

$$% \sum_i\; \frac{\partial {\vec g}}{\partial C_i}\;\frac{{\rm... ...{{\vec F}}} - \frac{{\rm d}{\bf {\vec{\Phi}}}}{{\rm d}t}\cdot$$ (11)

This equation describes the perturbed motion in terms of the orbital elements. Together with constraint (8) it constitutes a well-defined system which may be solved with respect to  ${\rm d}C_i/{\rm d}t$. An easy way to do this is to use the elegant trick suggested by Lagrange: to multiply the equation of motion by  $\partial{{\mbox{{$\vec f$ }}}}/\partial{C_n}$ and to multiply the constraint by  $-\partial{{\vec g} }/\partial{C_n}$. The former operation results in

$$% \frac{\partial {\mbox{{$\vec f$ }}}}{\partial C_n} ~\left( ... ...$ }}}}{\partial C_n}\;\frac{{\rm d} \bf \vec \Phi}{{\rm d}t},$$ (12)

while the latter gives

$$% -\frac{\partial \vec g}{\partial C_n}~ \left(\sum_j~\frac{... ...ht) = -\frac{\partial \vec g}{\partial C_n} \;{\bf\vec \Phi}.$$ (13)

Having summed these two equalities, we arrive at:

$$% \sum_j\;[C_n\;C_j]\;\frac{{\rm d}C_j}{{\rm d}t} = \frac{\... ... \frac{\partial {\vec g}}{\partial C_n} \;{\bf {\vec \Phi}},$$ (14)

$[C_n\;C_j]$ standing for the unperturbed (i.e., defined as in the two-body case) Lagrange brackets:

$$% [C_n\;C_j] \equiv \frac{\partial {{{\mbox{{$\vec f$ }}}}}}{... ...partial C_j}~ \frac{\partial {{{\vec g}}}}{\partial C_n}\cdot$$ (15)

It was selected above that $\mbox{{\boldmath$\vec \Phi$ }}$ is a function of time and of the "constants" C_i, but not of their time derivatives. Under this convention, the above equation may be shaped into a more convenient form by splitting the full derivative of  $\mbox{{\boldmath$\vec \Phi$ }}$ into a partial and a convective one, and by moving the latter to the left-hand side:

 

$$% \sum_j\;\left([C_n\;C_j] + \frac{\partial \mbox{{$\vec f$ }}}{\... ...ac{\partial \bf {\vec \Phi}}{\partial C_j}\right)~\frac{{\rm d}C_j}{{\rm d}t} = \frac{\partial {\mbox{{$\vec f$ }}}}{\partial C_n}\; {\Delta {\ve... ...}{\partial t} - \frac{\partial {\vec g}}{\partial C_n} \;{\bf {\vec \Phi}}\cdot$$ (16)

This is the most general form of the gauge-invariant perturbation equations of celestial mechanics. In the Lagrange gauge, when the  $\mbox{{\boldmath$\vec \Phi$ }}$ terms are absent, we can obtain an immediate solution for the individual ${\rm d}C_i/{\rm d}t$ by exploiting the well known expression for the Poisson-bracket matrix which is inverse to the Lagrange-bracket one and is offered in the literature for the two-body problem. (Be mindful that our brackets (15) are defined in the same manner as in the two-body case; they contain only functions  ${{\mbox{{$\vec f$ }}}}$ and $\vec g$ that are defined in the unperturbed, two-body, setting.) In an arbitrary gauge the presence of the term proportional to  $\partial{\mbox{{\boldmath$\vec \Phi$ }}}/\partial C_j$ on the left-hand side of (16) complicates the solution for  ${\rm d}C_i/{\rm d}t$, but only to the extent of requiring the resolution of a set of six simultaneous linear algebraic equations.

To draw to a close, we again emphasise that, for fixed interactions and initial conditions, all possible (i.e., compatible and sufficient) choices of gauge conditions expressed by the vector function  $\mbox{{\boldmath$\vec \Phi$ }}$ lead to a physically equivalent picture. In other words, the real trajectory is invariant under reparametrisations permitted by the ambiguity of the choice of gauge. This invariance has the following meaning. Suppose the equations of motion for  C₁, ..., C₆, with some gauge condition  ${\bf\vec \Phi}$ imposed, render the solution  C₁(t), ..., C₆(t). The same equations of motion, with another gauge  $\tilde{\mbox{{\boldmath$\vec \Phi$ }}}$enforced, furnish the solution  ${\tilde{C}}_{i}(t)$ that has a different functional form. Despite this difference, both solutions, C_i(t) and  ${\tilde{C}}_{i}(t)$, when substituted back in (5), yield the same curve  ${{\vec r}}(t)$ with the same velocities  $\dot{{\vec r}}(t)$. In mathematics this situation is called a fiber bundle, and it gives birth to a 1-to-1 map of C_i(t) onto  ${\tilde{C}}_{i}(t)$, which is merely a reparametrisation. In physics this map is called a gauge transformation. The entire set of these reparametrisations constitute a group of symmetry known as a gauge group.

Just as in electrodynamics, where the fields $\vec E$ and $\vec B$ stay invariant under gradient transformations of the four-potential $A^{\mu}$, so the invariance of the orbit implements itself through the form-invariance of expression (5) under the aforementioned map. The vector  $\mbox{{$\vec r$ }}$ and its full time derivative $\dot{\mbox{{$\vec r$ }}}$, play the role of the physical fields $\vec E$ and $\vec B$, while the Keplerian coordinates  C₁, ..., C₆play the role of the four-potential $A^{\mu}$.

No matter whether the role of constants C_i is played by the six Kepler variables  $\left(e,\;a,\;M_0,\;\omega,\;\Omega,\;{\it i}~\right)$or by some six combinations thereof (like, say, the Delaunay set (27) or the Jacobi set or the Poincare set), these constants are always called "orbital elements''. Only in the case when the condition (7) of Lagrange is imposed will these orbital elements be called "osculating''. This is called an osculating solution, because the particles' positions and velocities are instantaneously tangent to, i.e., touch, the Keplerian ellipses (or hyperbolae). It is easy to see that, by virtue of the differentiation chain rule, the osculation condition (7) remains form-invariant under a transition from the Kepler set of elements to any other set. This is most natural, because the osculation condition is purely geometric one and does not depend upon our preferences in the choice of convenient parameters.

A comprehensive discussion of all the above-raised issues can be found in Efroimsky (2003). The interconnection between the internal symmetry and multiple time scales, both in celestial mechanics and in a more general ODE context, is addressed in Newman & Efroimsky (2003).

2 Planetary equations in an arbitrary gauge

2.1 Lagrangian and Hamiltonian perturbations

We can proceed further by restricting the class of perturbations we consider to those in which $\Delta{{\vec F}}$ is derivable from a perturbed Lagrangian. Such a restricted class of disturbances is still broad enough to encompass most applications of celestial-mechanics perturbation theory. This happens largely because we shall consider Lagrangian perturbations dependent not only upon coordinates but also upon velocities. It will enable us to describe in a perturbative manner not only physical but also velocity-dependent, inertial forces emerging in non-inertial frames of reference. This, in its turn, will provide us with an opportunity to work in the coordinate system precessing with the planet of interest.

Let the unperturbed dynamics be described by the unit-mass Lagrangian ${\cal L}_0({{\vec r}},~{\dot{\vec r}},t) = {\dot{\vec r}}^{\left.~\right. 2}/2~-~U({\vec r},\;t)$, canonical momentum ${{\vec{p}}} = {\dot{\vec r}}$ and Hamiltonian ${\cal H}_0({{\vec r}},~{{\vec p}},~t) = {{\vec p}}^{\left.~\right. 2}/2~+~U({\vec r},t)$. The disturbed motion will be described by the new, perturbed, functions:

 

$$% {\cal L}({{\vec r}},~{\dot{\vec r}},~t)={\cal L}_0+ \Delta {\ca... ...}{2} - U({\vec r},~t) + \Delta {\cal L} ({\vec r}, ~{ { \dot { \vec {r}}}},~t),$$ (17)

 

$$% {{\vec p}} = {{\dot {\vec r}}} + \frac{\partial ~\Delta {\cal L}}{\partial {{\dot{\vec r}}}},$$ (18)

and

 

$${\cal H}({{\vec r}},~{{\vec p}},~t) = {\vec p}~{{\dot {\vec r}}} - {\cal L} = \frac{{\vec p}^{\left. ~\right. 2}}{2} + U + \Delta {\cal H},$$

$$\Delta {\cal H}({{\vec r}},~{{\vec p}},~t) \equiv - \Delta {\cal ... ...}~\left(\frac{\partial ~\Delta {\cal L}}{\partial {{\dot {\vec r}}}} \right)^2,$$ (19)

$\Delta {\cal H}$ being introduced as a variation of the functional form, i.e., as $\Delta {\cal H}~\equiv~{\cal H}({{\vec r}},~{{\vec p}},~t) - {\cal H}_0({{\vec r}},~{{\vec p}},~t)$. The Euler-Lagrange equations written for the perturbed Lagrangian (17) will give:

$$% {{\ddot {\vec r}}} = - \frac{\partial U}{\partial {{\vec r}}} + \Delta {\vec F},$$ (20)

where the new term

$$% \Delta {\vec F} \equiv \frac{\partial ~\Delta {\cal L}}{\pa... ...{\partial ~\Delta {\cal L}}{\partial {{\dot {\vec r}}}}\right)$$ (21)

is the disturbing force. Expression (21) reveals that whenever the Lagrangian perturbation is velocity-independent, it plays the role of the disturbing function. Generally, though, the disturbing force is not equal to the gradient of  $\Delta {\cal L}$, but has an extra term generated by the velocity dependence.

Examples in which a velocity dependent $\Delta {\cal L}$ has been used in a celestial mechanics setting include: the treatment of inertial forces in a coordinate system tied to the spin axis of a precessing planet (Goldreich 1965) and the velocity-dependent corrections to Newton's law of gravity in the relativistic two-body problem (Brumberg 1992).

Our next goal will be to employ the above formula (21) to translate the gauge-invariant perturbation Eq. (16) from the language of disturbing forces into that of Lagrangian disturbances.

2.2 Gauge-invariant planetary equations

When the expression (21) for the most general force emerging within the Lagrangian formalism is substituted into the gauge-invariant perturbation Eq. (14), it yields:

 

$$% \sum_j\;[C_n\;C_j]\;\frac{{\rm d}C_j}{{\rm d}t} = \frac{\partia... ...{\vec r}}}\right) - \frac{\partial {\vec g}}{\partial C_n} \;{\bf {\vec \Phi}},$$ (22)

$[C_n\;C_j]$ being the Lagrange brackets defined through (15). Since, for a velocity-dependent disturbance, the chain rule

 

$$% \frac{\partial ~\Delta {\cal L}}{\partial C_n} = \frac{\partial... ... ~\frac{\partial ({{\vec g}} + \mbox{{\boldmath$\vec \Phi$ }}) }{\partial C_n},$$ (23)

takes place, formula (22) may be reshaped into:

 

$$% \sum_j\;[C_n\;C_j]\;\frac{{\rm d}C_j}{{\rm d}t} = \frac{\partia... ...\Phi } } + \frac{\partial ~\Delta {\cal L}}{\partial \dot{\vec r}} \right)\cdot$$ (24)

Next we group terms so that the gauge function $\mbox{{\boldmath$\vec \Phi$ }}$ everywhere appears added to $\partial (\Delta {\cal L})/\partial {\dot{\vec r}}$, and bring the term proportional to  ${\rm d}C_j/{\rm d}t$to the left hand side of the equation:

 

$$\sum_j\;\left\{\;[C_n\;C_j] + \frac{\partial {\mbox{{$\vec f$ }}}... ...box{{\boldmath$\vec \Phi$ }}} \right)\;\right\} \frac{{\rm d}C_j }{{\rm d}t } =$$

$$\frac{\partial }{\partial C_n}~\left[\Delta {\cal L} + \frac{1}{2}~\left(\frac{\partial ~\Delta {\cal L}}{\partial \dot{\vec r}} \right)^2 \right]$$

$$-\left( \frac{\partial \vec g}{\partial C_n} + \frac{\partial \mb... ... \Phi}} + \frac{\partial ~\Delta {\cal L}}{\partial \dot{\vec r }} \right)\cdot$$ (25)

These modifications help us to recognise the special nature of the gauge $\mbox{{\boldmath$\vec \Phi$ }}= - \partial (\Delta {\cal L})/\partial {\dot{\vec r}}$, which will be the subject of discussion in the next subsection.

Contrast (25) with (16): while (16) expresses the variation-of-constants method in the most general form it can have in terms of the disturbing force $\Delta {\vec F}({\vec r},~{\dot{\vec r}},~t)$, Eq. (25) renders the most general form in the language of a Lagrangian perturbation  $\Delta {\cal L}({\vec r},~{\dot{\vec r}},~t)$.

The applicability of so generalised planetary equations in analytical calculations is complicated by the nontrivial nature of the left-hand sides of (16) and (25). Nevertheless, the structure of these left-hand sides leaves room for analytical simplification in particular situations. One such situation is when the gauge is chosen to be

 

$$% \mbox{{\boldmath$\vec \Phi$ }}= -\frac{\partial~\Delta {\cal L}}{\partial {\dot{\vec r}}} + \mbox{{\boldmath$\vec{\eta}$ }}({\vec{r}},~t),$$ (26)

$\mbox{{\boldmath$\vec{\eta}$ }}({\vec{r}},~t)$ being an arbitrary vector function linear in ${\vec r}$. (It may be, for example, proportional to ${\vec r}$ or may be equal, say, to a cross product of ${\vec r}$ by some time-dependent vector.) Under these circumstances the left-hand side in (25) reduces to the Lagrange brackets. The situation becomes especially simple when $\partial \Delta {\cal L}/\partial \dot{\vec{r}}$ happens to be linear in ${\vec r}$, in which case we may put $\mbox{{\boldmath$\vec \eta$ }}(\mbox{{$\vec r$ }},~t) = \partial \Delta {\cal L}/\partial \dot{\vec{r}}$ and, thus, employ the trivial Lagrange gauge $\mbox{{\boldmath$\vec \Phi$ }}= 0$ instead of the generalised Lagrange gauge. We shall encounter one such example in Sect. 3.4.

As already stressed above, the Lagrange brackets are gauge-invariant because functions ${{\mbox{{$\vec f$ }}}}$ and $\vec g$ are defined within the unperturbed, two-body, problem (1-3) that lacks gauge freedom. For this reason, one may exploit, to solve (25), the well-known expression for the inverse of matrix $[C_i\;C_j]$. Its elements are simplest (and are either zero or unity) when one chooses as the "constants'' the Delaunay set of orbital variables (Plummer 1918):

 

$${C}_i = \left\{L,\;G,\;H,\;M_0,\;\omega,\;\Omega\right\}$$

$$L \equiv \sqrt{\mu a},\;\; G \equiv \sqrt{\mu a \left(1 - e^2\right)},\;\; H \equiv \sqrt{\mu a \left(1 - e^2\right)}\;\cos {\it i},$$ (27)

where $\mu \equiv G(m_{\tiny {{\rm sun}}} + m_{\tiny {{\rm planet}}})$and $\left(e,\;a,\;M_0,\;\omega,\;\Omega,\;{\it i}\right)$are the Keplerian elements: e and a are the eccentricity and semimajor axis, M₀ is the mean anomaly at epoch, and the Euler angles  $\omega,\;\Omega,\;{\it i}$ are the the argument of pericentre, the longitude of the ascending node, and the inclination, respectively.

The simple forms of the Lagrange and Poisson brackets in Delaunay elements is the proof of these elements' canonicity in the unperturbed, two-body, problem: the Delaunay elements give birth to three canonical pairs  (Q_i, P_i) corresponding to a vanishing Hamiltonian: $(L,\;-~M_0),\;(G,\;-~\omega),\;(H,\;-~\Omega)$. In a perturbed setting, when only a position-dependent disturbing function  $R({\vec r},~t)$ is "turned on", it can be expressed through the Lagrangian and Hamiltonian perturbations in a simple manner, $R({\vec r},~t)~=~\Delta {\cal L}({\vec r},~t)~=~-~\Delta {\cal H} ({\vec r},~t)$, as can be seen from the formulae presented in the previous subsection. Under these circumstances, the Delaunay elements remain canonical, provided the Lagrange gauge is imposed (Brouwer & Clemence 1961). This long known fact can also be derived from our Eq. (25): if we put $\mbox{{\boldmath$\vec \Phi$ }}~=~0$ and assume $\Delta {\cal L}$ velocity-independent, we arrive to

 

$$% \sum_j\;[C_n\;C_j]\;\frac{{\rm d}C_j}{{\rm d}t} = \frac{\partial~ \Delta {\cal L} }{\partial C_n},$$ (28)

where

 

$$% \Delta {\cal L}\! = \! \Delta {\cal L}\left({\mbox{{$\vec f$ }}... ...(C,~t),\;t\right)= -\Delta {\cal H}\left({\mbox{{$\vec f$ }}}(C,~t),\;t\right).$$ (29)

This, in its turn, results in the well known Lagrange system of planetary equations, provided the parameters C_i are chosen as the Kepler elements. In case they are chosen as the Delaunay elements, then (28) leads to the standard Delaunay equations, i.e., to a symplectic system wherein the pairs $(L,\;-~M_0)$, $(G,\;-~\omega),\;(H,\;-~\Omega)$ again play the role of canonical variables, but the Hamiltonian, in distinction to the unperturbed case, no longer vanishes, instead being equal to  $\Delta {\cal H}\;=\;-\;\Delta {\cal L}$.

In our more general case, the perturbation depends also upon velocities (and, therefore, $\Delta {\cal L}$ is no longer equal to  $-\Delta {\cal H}$). Beside this, the gauge $\mbox{{\boldmath$\vec \Phi$ }}$ is set arbitrary. As demonstrated in Efroimsky (2003), under these circumstances the gauge-invariant Delaunay-type system is no longer canonical. However, it turns out that this system regains the canonical form in one special gauge, one that coincides with the Lagrange gauge when the perturbation bears no velocity dependence. The issue is explained at length in our previous paper (Efroimsky & Goldreich 2003). Here we offer a brief synopsis of this study.

2.3 Generalised Lagrange gauge wherein the Delaunay-type system is canonical

Equation (22) was cast in the form of (25) not only to demonstrate the special nature of the gauge

 

$$% \mbox{{\boldmath$\vec \Phi$ }}= -\frac{\partial~\Delta {\cal L}}{\partial {\dot{\vec r}}},$$ (30)

but also to single out the terms in the square brackets on the right-hand side of (25): together, these terms give exactly the Hamiltonian perturbation. Thus, we come to the conclusion that in the special gauge (30) our Eq. (25) simplifies to

 

$$% \sum_j\;[C_n\;C_j]\;\frac{{\rm d}C_j}{{\rm d}t} = -\frac{\partial\; \Delta {\cal H} }{\partial C_n}\cdot$$ (31)

As emphasised in the preceding subsection, the gauge invariance of definition (15) enables us to use the standard (Lagrange-gauge) expressions for the matrix inverse to [C_n C_j], to get the planetary equations from (31).

Comparing (31) with (28), we see that in the general case of an arbitrary $\Delta {\cal L}({\vec r},~{\dot{\vec r}},~t)$ one arrives from (31) to the same equations as from (28), except that now they contain $-\Delta {\cal H}$ instead of $\Delta {\cal L}$. When the orbit is parametrised by the Delaunay variables, those equations take the form:

 

$$% \frac{{\rm d}L}{{\rm d}t}= \frac{\partial ~\Delta {\cal H}}{\pa... ...frac{{\rm d}(-M_0)}{{\rm d}t} = - \frac{\partial ~\Delta {\cal H}}{\partial L},$$ (32)

 

$$% \frac{{\rm d}G}{{\rm d}t} = \frac{\partial ~\Delta {\cal H}}{\p... ...{{\rm d}(-\omega) }{{\rm d}t} = - \frac{\partial ~\Delta {\cal H}}{\partial G},$$ (33)

 

$$% \frac{{\rm d} H}{{\rm d}t} = \frac{\partial ~\Delta {\cal H}}{\... ...m d} (-\Omega)}{{\rm d}t} = - \frac{\partial ~\Delta {\cal H}}{\partial H}\cdot$$ (34)

which is a symplectic system. For this reason we name this special gauge the "generalised Lagrange gauge". In any different gauge $\mbox{{\boldmath$\vec \Phi$ }}$ the equations for the Delaunay variables would contain $\mbox{{\boldmath$\vec \Phi$ }}$-dependent terms and would not be symplectic. (Those gauge-invariant equations, for both Lagrange and Delaunay elements are presented in Efroimsky 2003 and Efroimsky & Goldreich 2003.) This analysis proves the following.

Theorem: though the gauge-invariant equations for Delaunay elements are, generally, not canonical, they become canonical in the "generalised Lagrange gauge".

That this Theorem is not merely a mathematical coincidence, but has deep reasons beneath it, will be shown in Sect. 4 where the subject is to be approached from the Hamilton-Jacobi viewpoint.

The above Theorem gives one example of the gauge formalism being of use: an appropriate choice of gauge can considerably simplify the planetary equations (in this particular case, it makes them canonical).

According to (18), the momentum can be written as:

$$% {\vec{p}} = {\dot{\vec{r}}} + \frac{\partial \Delta {\cal L... ... \frac{\partial \Delta {\cal L}}{\partial {{\dot{\vec{r}}}}},$$ (35)

which, in the generalised Lagrange gauge (30), simply reduces to:

$$% {\vec{p}} = {{\vec{g}}}.$$ (36)

Vector $\vec g$ was introduced back in (2-3) to denote the functional dependence of the unperturbed velocity upon the time and the parameters C_i. In the unperturbed, two-body, setting this velocity is equal to the momentum canonically conjugate to the position ${\vec r}$ (this is obvious from (18), for zero $\Delta {\cal L}$). This way, in the unperturbed case equality (36) is fulfilled trivially. The fact that it remains valid also under perturbation means that, in the said gauge, the canonical momentum in the disturbed setting is the same function of time and "constants" as in the unperturbed, two-body, case. Thus, we have shown, following Goldreich (1965), Brumberg et al. (1971), and Ashby & Allison (1993) that the instantaneous Keplerian ellipses (hyperbolae) defined in gauge (30) osculate the trajectory in phase space.

Not surprisingly, the generalised Lagrange gauge (30) reduces to $\mbox{{\boldmath$\vec \Phi$ }}= 0$ in the simple case of velocity-independent disturbances.

3 Gauge freedom and freedom of frame choice

3.1 Osculating ellipses described in different frames of reference

The essence of the variation-of-constants method in celestial mechanics is the following. A generic two-body-problem solution expressed by

$$\mbox{{$\vec r$ }}= {{\mbox{\vec{f}}}}\left(C,~t \right),$$ (37)

$$\left(\frac{\partial{\mbox{\vec{f}}} }{\partial t}\right)_{C}= {{\vec{g}}}\left(C,~t \right),$$ (38)

$$\left(\frac{\partial {{\vec{g}}}}{\partial t}\right)_{C}= - \frac{\mu}{f^2}\;\frac{{\mbox{\vec{f}}}}{f},$$ (39)

is employed as an ansatz to solve the disturbed problem:

 

$$% {{\vec{r}}}={{{{\mbox{$\vec{f}$ }}}}}(C(t),~t),$$ (40)

 

$$% {{\dot{\vec r}}}=\frac{\partial{{{{\mbox{$\vec{f}$ }}}}}}{\part... ...\;\frac{{\rm d} C_i}{{\rm d}t} = {{\vec g}} + {\mbox{{\boldmath$\vec \Phi$ }}},$$ (41)

 

$$% {{\ddot{\vec{r}}}} \frac{\partial {{\vec{g}}}}{\partial t} + \frac{\partial {{\vec g... ...\rm d} C_i}{{\rm d}t} + \frac{{\rm d} \mbox{{\boldmath$\vec \Phi$ }}}{{\rm d}t}$$ =

$${\left. ~ \right.}^{\left. ~ \right.} -\frac{\mu}{f^2}\;\frac{{\mbox{ {\vec{f}}}}}{f} + \frac{\partial ... ...} C_i}{{\rm d}t} + \frac{{\rm d} \mbox{{\boldmath$\vec \Phi$ }}}{{\rm d}t}\cdot$$ = (42)

As evident from (41), our choice of a particular gauge is equivalent to decomposition of the physical motion into a movement with velocity  $\vec g$ along the instantaneous ellipse (or hyperbola, in the flyby case), and a movement associated with the ellipse's (or hyperbola's) deformation that goes at the rate  $\mbox{{\boldmath$\vec \Phi$ }}$. It is then tempting to state that a choice of gauge is equivalent to a choice of an instantaneous comoving reference frame whereby we describe the motion. Such an interpretation is, however, incomplete. Beside the fact that we decouple the physical velocity in a certain proportion between $\vec g$ and  $\mbox{{\boldmath$\vec \Phi$ }}$, it also matters which physical velocity (i.e., velocity relative to what frame) is decoupled in this proportion. In other words, our choice of the gauge does not yet exhaust all freedom: we can still choose in which frame to write ansatz (40). We may write it in inertial axes or in some accelerated system. For example, in the case of a satellite orbiting an accelerated and precessing planet it is convenient to write the ansatz for the planet-related position vector.

The above kinematic formulae (40)-(42) do not yet contain information about our choice of the reference system in which we implement the variation-of-constants method. This information shows up at the next stage, when expression (42) is inserted into the dynamical equation of motion ${{\ddot{\mbox{{$\vec r$ }}}}}~=~-~(\mu\mbox{{$\vec r$ }}/r^3)~+~{\Delta \vec F}$ to yield:

$$% \frac{\partial {{\vec g}}}{\partial C_i}\;\frac{{\rm d} C_... ...l ~ \Delta {\cal L}}{\partial {{\dot {\vec r}}}}\right)\cdot$$ (43)

Complete information about the reference system in which we put the method to work (and, therefore, in which we define the orbital elements C_i) is contained in the expression for the perturbation force  $\Delta{{\vec F}}$. For example, if the operation is carried out in an inertial coordinate system, $\Delta{{\vec F}}$ contains physical forces solely. However, if we wish to implement the variation-of-constants approach in a frame moving with a linear acceleration $\bf\vec a$, then $\Delta{{\vec F}}$ also contains the inertial force  $-\bf\vec a$. In case this coordinate system rotates relative to inertial ones at a rate  $\mbox{{$\mu$ }}$, then  $\Delta{{\vec F}}$ also includes the inertial terms  $-2~{\mbox{{$\mu$ }}}$ $\times$ ${{\dot {\vec {r}}}}~-~{{\dot {\mbox{{$\mu$ }}}}}$ $\times$ ${{\vec {r}}}~-~{\mbox{{$\mu$ }}}$ $\times$ $({\mbox{{$\mu$ }}}$ $\times$ ${{\vec {r}}})$. In considering the motion of a satellite orbiting an oblate precessing planet it is most reasonable, though not obligatory, to apply the method (i.e., to define the time derivative) in axes that precess with the planet. However, this reasonable choice of coordinate system still leaves us with the freedom of gauge nomination. This will become clear in the example considered below in Sects. 3.2-3.5.

3.2 Relevant example

Gauge freedom of the perturbation equations of celestial mechanics finds an immediate practical implementation in the description of test particle motion around an precessing oblate planet (Goldreich 1965). It is trivial to extend this to account for acceleration of the planet's centre of mass.

Our starting point is the equation of motion in the inertial frame:

$$% {{{\vec r}~}}'' = \frac{\partial U}{\partial {\vec r}},$$ (44)

where U is the total gravitational potential and time derivatives in the inertial axes are denoted by primes. Suppose that the planet's spin axis precesses at angular rate $\mbox{{$\mu$ }}(t)$ and that the acceleration of its centre of mass is given by ${\vec{a}(t)}$. Physically, this acceleration originates both due to the circumsolar motion and due to the gravitational pull from the other planets (Krivov 1993). It is a very natural (and very common in celestial and planetary dynamics) technique to switch to a comoving or/and corotating reference frame, in order to better visualise the properties of the physical system. For example, in oceanography and atmospheric science almost all work is carried out in a frame corotating with the Earth. In our case, our preference will be to use not a corotating but rather a coprecessing frame, i.e., a coordinate system attached to the planet's centre of mass and precessing (but not spinning) with the planet. In the new coordinate frame the inertial forces modify the equation of motion so that it assumes the form:

$$% {{\ddot{\vec {r}}}} = \frac{\partial U}{\partial {{\vec{r}... ...} \times ({\mbox{{$\mu$ }}} \times {{\vec{r}}}) - {\vec{a}},$$ (45)

time derivatives in the accelerated frame being denoted by dots.

To implement the variation-of-constants approach in terms of the orbital elements defined in the accelerated frame, we note that the disturbing force on the right-hand side of (45) is generated according to (21) by:

$$% \Delta {\cal L}\left({{\vec r}},~{{\dot{\vec r}}},~t \righ... ...x{{$\mu$ }}} \times {{\vec{r}}}) - {\vec{a}\cdot{{\vec r}}},$$ (46)

where we denote by $R({\vec r},~t)$ the gravitational-potential perturbation (which is the perturbation of the overall gravitational potential U). Since

$$% \frac{\partial ~\Delta {\cal L}}{\partial {{\dot{r}}}} = {\mbox{{$\mu$ }}}\times{\vec{r}},$$ (47)

the corresponding Hamiltonian perturbation reads:

 

$$% \Delta {\cal H} = -\left[\Delta {\cal L} + \frac{1}{2}\left(\fr... ...\vec{r}}\times{\vec{p}})\cdot\mbox{{$\mu$ }}- {\vec{a}\cdot {\vec{ r}}}\right],$$ (48)

with vector ${\vec{J}} = {\vec{r}} \times {\vec{p}}$ being the satellite's orbital angular momentum in the inertial frame.

According to (35) and (47), the momentum can be written as:

$$% {\vec{p}} = {{\vec{g}}} + \mbox{{\boldmath$\vec \Phi$ }}+ {\mbox{{$\mu$ }}}\times{\mbox{{$\vec f$ }}},$$ (49)

whence the Hamiltonian perturbation becomes

 

$$% \Delta {\cal H} = - \left[R + \left({\mbox{{$\vec f$ }}}\times ... ...times{\mbox{{$\vec f$ }}} \right) - {\vec{a}\cdot {\mbox{{$\vec f$ }}}}\right].$$ (50)

3.3 Elements defined in an accelerated, rotating frame, that osculate in the comoving inertial frame

In this subsection we recall a calculation carried out by Goldreich (1965), Brumberg et al. (1971), and Ashby & Allison (1993) and demonstrate that it may be interpreted as an example of nontrivial gauge fixing.

Let us implement the variation-of-constants method in a frame that is accelerating at rate $\vec a$ and rotating at angular rate  $\mbox{{$\mu$ }}$ relative to some inertial system S. This means that, in the variation-of-constants Eq. (43), $\Delta {\cal L}$ is given by formula (46) and  $\Delta {\cal H}$ by (50).

We now choose to describe the motion in the generalised Lagrange gauge (30), so the expression $\left( \mbox{{\boldmath$\vec \Phi$ }} + {\mbox{{$\mu$ }}}\times{ \vec{r}} \right)$ on the right-hand side of (50) vanishes (as follows from (47)), and the expression for $\Delta {\cal H}$ in terms of ${{\mbox{{$\vec f$ }}}}$ and $\vec g$ has the form:

$$% \Delta {\cal H} = - \left[R({{\mbox{{$\vec f$ }}}},~t)+ \mb... ...{{\vec{g}}}) - {\vec{{a}}}\cdot{{\mbox{{$\vec f$ }}}}\right].$$ (51)

At the same time, the generic expression for the variation-of-constants method given by (25) simplifies to (31). Insertion of (51) therein leads us to

$$% [C_r\;C_i]\;\frac{{\rm d}C_i}{{\rm d}t} = \frac{\partial}{\... ... {{\vec g}}) - {\vec{a}}\cdot{{\mbox{{$\vec f$ }}}}\right].$$ (52)

Interestingly, this equation does not contain $\dot{\mbox{{$\mu$ }}}$ even though it is valid for nonuniform precession.

As explained in Sect. 2.3, in the generalised Lagrange gauge the vector $\vec g$ is equal to the canonical momentum ${\vec p} = {\dot{\vec{r}}} + {\partial~ \Delta {\cal L}}/{\partial {\dot{\vec{r}}}}$. In the case when the velocity dependence of $\Delta {\cal L}$ is called into being by inertial forces, the momentum is, according to (47),

$$% {\vec p} = {\dot{\vec{r}}} + \frac{\partial~ \Delta {\cal ... ...ec{r}}}} = {\dot{\vec{r}}} + \mbox{{$\mu$ }}\times {\vec{r}},$$ (53)

which is the particle's velocity relative to the inertial frame comoving with the accelerated, rotating frame. In this sense we may say that our elements are defined in the accelerated, rotating frame, but osculate in the comoving inertial one.

In the Appendix we provide explicit expression for each of the partial derivatives of $\mbox{{$\mu$ }}\cdot{\vec{J}}$ that appears in the planetary Eq. (52).

3.4 Elements defined in the accelerated, rotating frame, that osculate in this frame

Here we not only define the elements in the accelerated, rotating frame, but we also make them osculate in this system, i.e., we make them satisfy $\mbox{{\boldmath$\vec \Phi$ }}= 0$. In this gauge, expression (50) takes the following form:

 

$$% \Delta {\cal H} = - \left[R({{\mbox{{$\vec f$ }}}},~t) + \mbox{... ...\mu$ }}\times\mbox{{$\vec f$ }}) - {\vec{a}}\cdot{{\mbox{{$\vec f$ }}}}\right],$$ (54)

while Eq. (25), after some algebra, looks like this:

 

$$% [C_n\;C_i]\;\frac{{\rm d}C_i}{{\rm d}t} - \frac{\partial~\Delta {\cal H}}{\partial C_n}$$ =

$$+\mbox{{$\mu$ }}\cdot \left(\frac{\partial{\mbox{{$\vec f$ }}}}{\... ... g}} - {\mbox{{$\vec f$ }}}\times \frac{\partial{\vec{g}}}{\partial C_n}\right)$$

$$- {\bf {\dot{\mbox{{$\mu$ }} }}}\cdot\left(\mbox{{$\vec f$ }}\tim... ...c{\partial }{\partial C_n}\left(\mbox{{$\mu$ }}\times\mbox{{$\vec f$ }}\right).$$ (55)

When substituting (54) into (55), it is convenient to rent the expression for $\Delta {\cal H}$ apart and to group the term $(\mbox{{$\mu$ }}\times{\mbox{$\vec{f}$ }})\cdot(\mbox{{$\mu$ }}\times\mbox{{$\vec f$ }})$ with the last term on the right-hand side of (55):

  $$\begin{eqnarray*}[C_n\;C_i]\;\frac{{\rm d}C_i}{{\rm d}t} &=& \frac{\partial}{\... ...rtial C_n}\left(\mbox{{$\mu$ }}\times\mbox{{$\vec f$ }}\right). \end{eqnarray*}$$

In so writing (56) we have deliberately cast it into a form that eases comparison with (52).

In the Appendix we set up an apparatus from which the partial derivatives of the inertial terms with respect to the orbital elements may be obtained. We also show that some of these derivatives vanish. However, a complete evaluation of the inertial input to the planetary equations in the ordinary Lagrange gauge involves a long and tedious calculation which we do not carry out.

3.5 Comparison of the two gauges

One of the powers of gauge freedom lies in the availability of gauge choices that simplify the planetary equations, as we can see from contrasting (52) with (56). While the latter equation is written under the customary Lagrange constraint (i.e., for elements osculating in the frame where they are defined), the former equation is written under a nontrivial constraint called the "generalised Lagrange gauge''. The simplicity of (52), in contrast with (56), is evident.

By identifying the parameters C_i with the Delaunay variables, one arrives from (52) and (56) to the appropriate Delaunay-type equations (see Appendix I to Efroimsky & Goldreich 2003). The Lagrange equations corresponding to (52) and to (56) may be easily derived from each of these two equations by choosing C_i as the Kepler elements and using the appropriate Lagrange brackets.

Although the planetary equations are much simpler in the generalised Lagrange gauge than in the ordinary Lagrange gauge, some of these differences are less important than others. In many physical situations, though not always, the ${{\mbox{{$\mu$ }}}^2}$and $\dot{\mbox{{$\mu$ }}}$ terms in (56) are of a higher order of smallness compared to those linear in $\mbox{{$\mu$ }}$, and therefore may be neglected, at least for sufficiently short times.

3.6 Further applications

In the above example of a satellite orbiting a wobbling planet, an evident simplification of the planetary equations (both in the Lagrange and Delaunay forms) was achieved through imposition of the generalised Lagrange gauge. This optimal gauge differed from the standard Lagrange constraint, because in the said example the Lagrangian perturbation depended upon velocities.

A similar situation emerges in the relativistic two-body problem. In the relativistic dynamics, even the two-body problem is perturbed, the role of perturbation being played by the relativistic correction to the Newton law of gravity. Interestingly, this correction depends not only upon the positions, but also upon the velocities of the binary components (Brumberg 1992). For this reason, to simplify the orbit integration of a binary, the generalised Lagrange gauge (not the customary Lagrange constraint) should be imposed. In this gauge the calculations will be very considerably simplified (and, for example, it is in the generalised Lagrange gauge that the equations for the Delaunay elements will retain their canonicity).

Another simple example is a non-relativistic reduced two-body problem with a variable mass. In this case, too, the Lagrangian acquires a velocity-dependent correction. Hence, in this case, the orbital elements will be convenient to introduce in the generalised Lagrange gauge, not in the customary Lagrange gauge.

4 Planetary equations and gauges in the hamilton-Jacobi approach

In this section we demonstrate that the derivation of planetary equations in the N-particle ($N \geq 3$) case, performed through the medium of Hamilton-Jacobi method, implicitly contains a gauge-fixing condition not visible to the naked eye. We present a compact account of our study; a comprehensive description containing technical details may be found in Efroimsky & Goldreich (2003).

The Hamilton-Jacobi analysis rests on the availability of different canonical descriptions of the same physical process. Any two such descriptions, $(q,~p,~{\cal H}(q,p))$ and $(Q,~P,~{\cal H}^{*}(Q,P))$, correspond to different parametrisations of the same phase flow, and both obey Hamilton's equations. Due to the latter circumstance the infinitesimally small variations

$$% {\rm d} \theta = p~{\rm d}q~-~{\cal H}~{\rm d}t$$ (56)

and

$$% {\rm d} \tilde{\theta} = P~{\rm d}Q~-~{\cal H}^{*}~{\rm d}t$$ (57)

are perfect differentials, and so is their difference

$$% -{\rm d}W \equiv {\rm d} \tilde{\theta} - {\rm d} \theta =... ...- p~{\rm d}q - \left({\cal H}^{*} - {\cal H}\right)~{\rm d}t.$$ (58)

Here, vectors q, p, Q, and P each contain Ncomponents. Given a phase flow parametrised by a set $(q,~p,~{\cal H}(q,p,t))$, it is always useful to simplify the description by a canonical transformation to a new set $(Q,~P,~{\cal H}^{*}(Q,P,t))$, with the new Hamiltonian ${\cal H}^{*}$ being constant in time. Most advantageous are transformations that nullify the new Hamiltonian ${\cal H}^{*}$, because then the new canonical equations render the variables (Q, P) constant. A powerful method of generating such transformations stems from (58) being a perfect differential. It is sufficient to consider W to be a function of the time and only two other canonical variables, for example q and Q. Then (58) may be written as:

 

$$% -\frac{\partial W}{\partial t}\;{\rm d}t - \frac{\partial W}{\partial Q}\;{\rm d}Q - \frac{\partial W}{\partial q}\;{\rm d}q = P~{\rm d}Q - p~{\rm d}q + \left({\cal H} - {\cal H}^{*}\right)~{\rm d}t$$ (59)

from which it follows that

 

$$P = -\frac{\partial W}{\partial Q},\;\;\;\;\;\;\;p = \frac{\partial W}{\partial q},$$

$${\cal H}(q,~p,~t) + \frac{\partial W}{\partial t} = {\cal H}^{*}(Q,~P,~t).$$ (60)

Inserting the second equation into the third and assuming that ${\cal H}^{*}(Q,~P,~t)$ is simply a constant, we get the famous Jacobi equation:

$$% {\cal H}\left(q,~\frac{\partial W}{\partial q},~t\right) + \frac{\partial W}{\partial t} = {\cal H}^{*}$$ (61)

whose solution furnishes the transformation-generating function W. The elegant power of the method becomes most visible if the constant  ${\cal H}^{*}$ is set to zero. Under this assumption the reduced two-body problem is easily resolved. Starting with the three spherical coordinates and their canonical momenta as (q, p), one arrives to canonically conjugate constants (Q, P) that coincide with the Delaunay elements (27):
$(Q_1,\;P_1) = (L, - M_0); \;\;(Q_2,\;P_2) = (G,\;-\omega);\\ (Q_3,\;P_3) = ({H},\;-\Omega)$.

Extension of this approach to the N-particle problem begins with consideration of a disturbed two-body setting. The number of degrees of freedom is still the same (three coordinates qand three conjugate momenta p), but the initial Hamiltonian is perturbed:

$$% \dot q = \frac{\partial ({\cal H} + \Delta {\cal {H}})}{\pa... ...rac{\partial ({\cal H} + \Delta {\cal {H}})}{\partial q}\cdot$$ (62)

While in (59)-(61) one begins with the initial Hamiltonian ${\cal H}$ and ends up with ${\cal H}^{*} = 0$, the method may be extended to the perturbed setting by accepting that now we start with a disturbed initial Hamiltonian ${\cal H} + \Delta {\cal {H}}$ and arrive, through the same canonical transformation, to an equally disturbed eventual Hamiltonian ${\cal H}^{*} + \Delta {\cal {H}} = \Delta {\cal {H}}$. Plugging these new Hamiltonians into (59) leads to cancellation of the disturbance  $\Delta {\cal H}$ on the right-hand side, whereafter one arrives to the same equation for W(q, Q, t) as in the unperturbed case. Now, however, the new canonical variables are no longer conserved but obey the dynamical equations:

$$% \dot Q = \frac{\partial ~\Delta {\cal {H}}}{\partial P},\;... ...\dot P = - \frac{\partial ~\Delta {\cal H}}{\partial Q}\cdot$$ (63)

Because the same generating function is used in the perturbed and unperturbed cases, the new, perturbed, solution (q, p) is expressed through the perturbed "constants" Q(t) and P(t) in the same manner as the old, undisturbed, q and p were expressed through the old constants Q and P. This form-invariance provides the key to the N-particle problem: one should choose the transformation-generating function W to be additive over the particles and repeat this procedure for each of the bodies, separately.

Armed with this preparation, we can proceed to uncover the implicit gauge choice made in using the Hamilton-Jacobi method to derive evolution equations for the orbital elements. To do this we substitute the equalities:

$$% \dot{ Q} = \frac{\partial ~\Delta {\cal {H}} }{ \partial P ... ...ta {\cal {H}}}{\partial p }\;\frac{\partial p }{ \partial P }$$ (64)

and

$$% \dot{P } = - \frac{\partial ~\Delta {\cal {H}}}{\partial Q}... ...{\cal {H}} }{ \partial p } \;\frac{\partial p }{ \partial P }$$ (65)

into the expression for the velocity:

$$% \dot{ q } = \frac{\partial q}{\partial t} + \frac{\partial ... ...tial Q }\;\dot{Q} + \frac{\partial q }{\partial P }\;\dot{P}.$$ (66)

This leads to:

 

$$% \dot{ q} \frac{\partial q}{\partial t} + \left( \frac{\partial q}{\partial... ...partial q }{\partial Q }\right)\;\frac{\partial ~\Delta {\cal {H}}}{\partial q}$$ =

$$+ \left(\frac{\partial q}{\partial Q }\;\frac{\partial p}{\partia... ...partial p }{\partial Q }\right)\;\frac{\partial ~\Delta {\cal {H}}}{\partial p}$$

$$g + \left( \frac{ \partial ~ \Delta {\cal {H}} }{ \partial p} \right)_{q,~t},\;\;\;\;\;\;\;\;\;\;\;\;\;g \equiv \frac{\partial q}{\partial t},$$ = (67)

where we have taken into account that the Jacobian of the canonical transformation is unity:

$$% \frac{\partial q}{\partial Q}\;\frac{\partial p}{\partial P... ...partial q }{\partial P}\;\frac{\partial p }{\partial Q } = 1.$$ (68)

To establish the link between the regular variation-of-constsnts method and the canonical treatment, compare (67) with (41). We see that the symplectic description necessarily imposes a particular gauge $\Phi = {\partial ~\Delta {\cal {H}}}/{\partial p}$.

It can be easily demonstrated that this special gauge coincides with the generalised Lagrange gauge (30) discussed in Sect. 2.2. To that end one has to compare the Hamilton equation for the perturbed Hamiltonian (19),

$$% \dot{q} = \frac{\partial ~\left({\cal H} + \Delta {\cal {H}... ...al p} = p + \frac{\partial ~\Delta {\cal {H}}}{\partial p},$$ (69)

with the definition of momentum from the Lagrangian (17),

$$% p \equiv \frac{\partial~\left({\cal L}(q,~\dot{q},~t) + \De... ...{q} + \frac{\partial ~\Delta {\cal L}}{\partial \dot{q}}\cdot$$ (70)

Equating the above two expressions immediately yields:

$$% \Phi \equiv \left( \frac{ \partial ~ \Delta {\cal {H}} }{ \... ...ac{\partial\;\Delta {\cal L}}{\partial \dot{q}}\right)_{q,~t}$$ (71)

which coincides with (30). Thus, the transformation generated by W(q, Q, t) is canonical only if the physical velocity $\dot q$ is split in a fashion prescribed by (67), i.e., if (71) is fulfilled. This is exactly what our Theorem from Sect. 2.2 says.

To summarise, the generalised Lagrange constraint, $\Phi = - {\partial ~\Delta {\cal L}}/{\partial \dot q}$, is tacitly instilled into the Hamilton-Jacobi method. Simply by employing this method (at least, in its straightforward form), we automatically fix the gauge. By sticking to the Hamiltonian description we sacrifice gauge freedom.

Above, in Sect. 2.3, we established that in the generalised Lagrange gauge the momentum coincides with $\vec g$. We now can get to the same conclusion from (67), (70) and (71):

$$% p \equiv \frac{\partial~\left({\cal L}(q,~\dot{q},~t) + \De... ...q,~\dot{q},~t)\right)}{\partial \dot q} = \dot{q} - \Phi = g.$$ (72)

Thus, implementation of the Hamilton-Jacobi theory in celestial mechanics demands the orbital elements to osculate in phase space. Naturally, this demand reduces to that of regular osculation in the simple case of velocity-independent $\Delta {\cal L}$.

5 Conclusions

In the article thus far we have studied the topic recently raised in the literature: the planetary equations' internal symmetry that stems from the freedom of supplementary condition's choice. The necessity of making such a choice constrains the trajectory to a 9-dimensional submanifold of the 12-dimensional space spanned by the orbital elements and their time derivatives. Similarly to the field theory, the choice of the constraint (=the choice of gauge) is vastly ambiguous and reveals a hidden symmetry instilled in the description of the N-body problem in the language of orbital elements.

We addressed the issue of internal freedom in a sufficiently general setting where a perturbation to the two-body problem depends not only upon positions but also upon velocities. Such situations emerge when relativistic corrections to Newton's law are taken into account or when the variation-of-constants method is employed in rotating systems of reference.

We derived the most general form of the gauge-invariant perturbation equation of celestial mechanics, written in terms of a disturbing force. Then we transformed it into the most general gauge-invariant perturbation equation expressed through the Lagrangian disturbance.

Just as a choice of an appropriate gauge simplifies solution of the equations of motion in electrodynamics, an alternative (to that of Lagrange) choice of gauge in the celestial mechanics can simplify orbit calculations. We provided one such example, a satellite orbiting a precessing planet. In this example, the choice of the generalised Lagrange gauge considerably simplifies matters. To achieve this simplification, we not only exercised our right to choose a convenient gauge, but we also chose a preferred coordinate system in which to implement the variation-of-constants method. This interplay of the two types of freedom enabled us to eliminate some of the mathematical complications associated with the inertial forces. Not surprisingly, it has turned out to be convenient to define the orbital elements in the precessing frame of the planet; however, for the sake of mathematical simplification, it also turned out to be beneficial to make these elements osculate in a different, inertial frame of reference.

We have explained where the Lagrange constraint tacitly enters the Hamilton-Jacobi derivation of the Delaunay equations. This constraint turns out to be an inseparable (though not easily visible) part of the method: in the case of momentum-independent disturbances, the N-body generalisation of the two-body Hamilton-Jacobi technique is legitimate only if we use orbital elements that are osculating. In the situation where the disturbance depends not only upon positions but also upon velocities, another constraint (which we call the "generalised Lagrange constraint") turns out to be stiffly embedded in the Hamilton-Jacobi development of the problem.

Unless a specific constraint (gauge) is imposed by hand, the planetary equations assume their general, gauge-invariant, form. In the case of a velocity-independent disturbance, any gauge different from that of Lagrange drives the Delaunay system away from its symplectic form. If we permit the disturbing force to depend also upon velocities, the Delaunay equations retain their canonicity only in the generalised Lagrange gauge. Interestingly, in this special gauge the instantaneous ellipses (hyperbolae) osculate in phase space.

Briefly speaking, N-body celestial mechanics, expressed in terms of orbital elements, is a gauge theory, but it is not strictly canonical. It becomes canonical in the generalised Lagrange gauge.

Acknowledgements

The authors are grateful to William Newman and Victor Slabinski for their help in improving the manuscript, to Sergei Klioner for drawing the authors' attention to the paper by Brumberg et al., and to Michael Nauenberg for the information on Newton's priority in developing the variation-of-constants method. Research by ME was supported by NASA grant W-19948. Research by PG was partially supported by NSF grant AST 00-98301.

Appendix

In this Appendix we set up an apparatus from which one may evaluate the partial derivatives with respect to the orbital elements of inertial terms that appear in the planetary equations derived in Sects. 3.3 and 3.4. We then show that some of these derivatives vanish. Following that, we derive explicit expressions for each derivative of $\mbox{{$\mu$ }}\cdot\left(\mbox{{$\vec f$ }}\times{\vec g}\right)$, which provides a complete analytic evaluation of the rotational input in the generalised Lagrange gauge. The topic is further developed (and the appropriate generalised Lagrange system of equations is presented) in Efroimsky & Goldreich (2003).

To find the explicit form of the dependence $\mbox{{$\vec f$ }}= \mbox{{$\vec f$ }}(C_i,~t)$, it is conventional to introduce an auxiliary set of Cartesian coordinates ${\vec q}$, with an origin at the gravitating centre, and with the first two axes located in the plane of orbit. The ${\vec q}$coordinates are easy to express through the semimajor axis a, the eccentricity e and the eccentric anomaly E:

 

$$% q_1 \equiv a\left(\cos E - e\right),\;\;q_2 \equiv a\sqrt{1 - e^2}\;\sin E,\;\;q_3 = 0,$$ (73)

where E itself is a function of the semimajor axis a, the eccentricity e, the mean anomaly at epoch, M₀, and the time, t. The time dependence is realised through the Kepler equation

$$% E - e\;\sin E = M,$$ (74)

where

$$% M \equiv M_0 + \mu^{1/2}~\int_{t_0}^t \;a^{-3/2}\;{\rm d}t.$$ (75)

The inertial-frame-related position of the body reads:

 

$$% {{\mbox{{$\vec r$ }}}} = \mbox{{$\vec f$ }} \left(\Omega,~{\it ... ...{R}}}(\Omega,~{\it i},~\omega)\;~{{\vec q}}\left(a,~e,\;E(a,~e,~M_0,~t)\right),$$ (76)

$\vec{\hat R}(\Omega,~{\it i},~\omega)$ being the matrix of rotation from the orbital-plane-related coordinate system ${\vec q}$ to the fiducial frame (x, y, z) in which the vector ${\vec r}$ is defined. This rotation is parametrised by the three Euler angles: inclination, ${\it i}$; the longitude of the node, $\Omega$; and the argument of the pericentre, $\omega$.

In the unperturbed two-body setting the velocity is expressed by

 

$$% {{\vec{g}}} = \frac{\partial }{\partial t} \;\mbox{{$\vec f$ }}... ...t i},~\omega)\;\;\left(\frac{\partial {\vec q}}{\partial E}\right)_{a,~ e}\cdot$$ (77)

One can similarly calculate partial derivatives of  ${{\mbox{{$\vec f$ }}}}$with respect to M₀:

 

$$% \frac{\partial }{\partial M_0} \;\mbox{{$\vec f$ }} \left(\Omeg... ...~{\it i},~\omega)\;\;\left(\frac{\partial {\vec q}}{\partial E}\right)_{a,~ e},$$ (78)

whence it becomes evident that $\partial \mbox{{$\vec f$ }}/\partial M_0$is parallel to $\vec g$ and, hence,

$$% {{\vec g}} \times \left(\frac{\partial \mbox{{$\vec f$ }}}{\partial M_0}\right)_{\Omega,~{\it i},~\omega,~a,~e,~t} = 0.$$ (79)

By a similar trick it is possible to demonstrate that $\partial (\mbox{{$\vec f$ }}\times{\vec g})/\partial M_0$ is proportional to $\partial (\mbox{{$\vec f$ }}\times{\vec g})/\partial E$ and, therefore, to $\partial (\mbox{{$\vec f$ }}\times{\vec g})$/ $\partial t$. Hence, this derivative vanishes (because in the two-particle case the cross product $\mbox{{$\vec f$ }}\times{\vec g}$is an integral of motion). This vanishing of $\partial (\mbox{{$\vec f$ }}\times{\vec g})/\partial M_0$, along with (A.7), implies:

$$% {{\mbox{{$\vec f$ }}}} \times \left(\frac{\partial \vec g}{\partial M_0}\right)_{\Omega,~{\it i},~\omega,~a,~e,~t} = 0.$$ (80)

In the situation when the parameters C_i are implemented by the Delaunay elements, a similar sequence of calculations leads to

$$% {{\vec g}} \times \left(\frac{\partial \mbox{{$\vec f$ }}}{\partial M_0}\right)_{\Omega,~\omega,~L,~G,~{H},~t} = 0$$ (81)

and, appropriately, to:

$$% {{\mbox{{$\vec f$ }}}} \times \left(\frac{\partial \vec g}{\partial M_0}\right)_{\Omega,~\omega,~L,~G,~{H},~t} = 0.$$ (82)

We can proceed much farther in the generalised Lagrange gauge, at least in so far as derivatives of the rotational input ${\vec J}/\mu = \mbox{{$\vec f$ }}\times\mbox{{$\vec{g}$ }}$ are concerned. (We remind that here and everywhere $\mu$ stands for the reduced mass, while $\mbox{{$\mu$ }}$ denotes the precession rate.)

As we proved above, this cross product is independent of M₀ and, hence,

$$% \mbox{{$\mu$ }}\;\cdot\;\frac{\partial(\mbox{{$\vec f$ }}\times\mbox{{$\vec{g}$ }})}{\partial M_0} = 0.$$ (83)

Since ${\vec J}$ is orthogonal to the orbit plane, it is invariant under rotations of the orbit within its plane, whence

$$% \mbox{{$\mu$ }}\;\cdot\;\frac{\partial(\mbox{{$\vec f$ }}\times\mbox{{$\vec{g}$ }})}{\partial \omega} = 0.$$ (84)

To continue, we note that in the two-body setting the ratio ${\vec J}/\mu$, is known to be equal to $\sqrt{a(1~-~e^2)}\;{\vec{\hat{w}}}$ where $\vec{\hat{w}}$ is a unit vector perpendicular to the unperturbed orbit's plane. Moreover, in planet-associated noninertial axes (x, y, z) with corresponding unit vectors $({\vec{\hat{x}}},~{\vec{\hat{y}}},~{\vec{\hat{z}}})$, the normal to the orbit is expressed by:

$$% {\vec{\hat{w}}} = {\vec{\hat{x}}}\;\sin {\it i}\;\sin \Omeg... ...\;\sin {\it i}\;\cos \Omega + {\vec{\hat{z}}}\;\cos {\it i}.$$ (85)

Hence,

$$% \mbox{{$\mu$ }}~\cdot~ \frac{\partial~\left(\mbox{{$\vec f... ...{\partial a}=\frac{1}{2}~\sqrt{\frac{1~-~e^2}{a}}~~\mu_\perp$$ (86)

and

$$% \mbox{{$\mu$ }}~\cdot~ \frac{\partial~\left(\mbox{{$\vec f... ...ight)}{\partial e}= -\sqrt{\frac{a~e^2}{1~-~e^2}}~~\mu_\perp$$ (87)

where $\mu_\perp = \mu_x\;\sin {\it i}\;\sin \Omega - \mu_y\;\sin {\it i}\;\cos \Omega + \mu_z\;\cos {\it i}$ is the orthogonal-to-orbit component of the precession rate. The remaining two derivatives look:

 

$$\mbox{{$\mu$ }}\;\cdot\;\frac{\partial~\left(\mbox{{$\vec f$ }}\t... ...t)~}\;\;\mbox{{$\mu$ }}\;\cdot\;\frac{\partial \vec{\hat w}}{\partial \Omega} =$$

$$\sqrt{a~\left(1~-~e^2 \right)~}\;\;\left\{ \mu_x\;\sin {\it i}\;\cos \Omega + \mu_y\;\sin {\it i}\;\sin \Omega \right\}$$ (88)

and

 

$$\mbox{{$\mu$ }}\;\cdot\;\frac{\partial~\left(\mbox{{$\vec f$ }}\t... ...)~}\;\;\mbox{{$\mu$ }}\;\cdot\;\frac{\partial \vec{\hat w}}{\partial {\it i}} =$$

$$\sqrt{a~\left(1-e^2 \right)}~\left\{ \mu_x~\cos {\it i}~\sin \Omega - \mu_y~\cos {\it i}~\cos \Omega - \mu_z~\sin {\it i} \right\}.$$ (89)

As for the derivatives of ${\vec a}~\cdot~{\mbox{{$\vec f$ }}}$, they may be calculated directly from the expression for $\mbox{{$\vec f$ }} (\Omega,~\omega,~{\it i},~a,~e,~M_0;\;t)$ presented above. However, the resulting expressions are cumbersome so we do not present them here.

References

 
• Ashby, N., & Allison, T. 1993, Cel. Mech. Dyn. Astron., 57, 537 In the text NASA ADS  
• Borderies, N., & Longaretti, P. Y. 1987, Icarus, 72, 593 NASA ADS  
• Brouwer, D., & Clemence, G. M. 1961, Methods of Celestial Mechanics (NY & London: Academic Press), Chapter XI In the text  
• Brumberg, V. A., Evdokimova, L. S., & Kochina, N. G. 1971, Cel. Mech., 3, 197 In the text  
• Brumberg, V. A. 1992, Essential Relativistic Celestial Mechanics (Bristol: Adam Hilger) In the text  
• Efroimsky, M. 2002, Equations for the orbital elements. Hidden symmetry. Preprint No. 1844 of the Institute of Mathematics and its Applications, University of Minnesota. In the text http://www.ima.umn.edu/preprints/feb02/1844.pdf  
• Efroimsky, M. 2003, [astro-ph/0212245] In the text  
• Efroimsky, M., & Goldreich, P. 2003, J. Math. Phys., 44, 5958 In the text  
• Euler, L. 1748, Recherches sur la question des inégalites du mouvement de Saturne et de Jupiter, sujet proposé pour le prix de l'année, Berlin. For modern edition see: Euler L. Opera mech. et astron. (Birkhauser-Verlag, Switzerland, 1999) In the text  
• Euler, L. 1753, Theoria motus Lunae exhibens omnes ejus inaequalitates etc. Impensis Academiae Imperialis Scientarum Petropolitanae (St. Petersburg, Russia). For modern edition see: Euler L. Opera Mech. Astron. (Switzerland: Birkhauser-Verlag, 1999) In the text  
• Goldreich, P. 1965, AJ, 70, 5 In the text NASA ADS  
• Jacobi, C. G. J. 1866, Vorlesungen $\ddot{\rm u}$ber Dynamik (Berlin: G. Reimer) In the text  
• Kovalevsky, J. 1967, Introduction to Celestial Mechanics (Dordrecht: Reidel)  
• Krivov, A. 1993, Astron. Rep. (Sov. Astron.), 37, 109 In the text NASA ADS  
• Lagrange, J.-L. 1778, Sur le Problème de la détermination des orbites des comètes d'après trois observations, 1 $^{{\rm er}}$ et 2 $^{\rm i\grave{e}me}$ mémoires (Nouveaux Mémoires de l'Académie de Berlin), 1778, later edition: in uvres de Lagrange. IV (Paris: Gauthier-Villars, 1869) In the text  
• Lagrange, J.-L. 1783, Sur le Problème de la détermination des orbites des comètes d'après trois observations, 3 $^{i\grave{e}me}$ mémoire (Ibidem), 1783, later edition: in uvres de Lagrange. IV (Paris: Gauthier-Villars, 1869) In the text  
• Lagrange, J.-L. 1808, Sur la théorie des variations des éléments des planètes et en particulier des variations des grands axes de leurs orbites, lu le 22 août 1808 à l'Institut de France, later edition: in uvres de Lagrange. VI, 713 (Paris: Gauthier-Villars, 1877) In the text  
• Lagrange, J.-L. 1809, Sur la théorie générale de la variation des constantes arbitraires dans tous les problèmes de la mécanique, lu le 13 mars 1809 à l'Institut de France, later edition: in uvres de Lagrange. VI, 771 (Paris: Gauthier-Villars, 1877) In the text  
• Lagrange, J.-L. 1810, Second mémoire sur la théorie générale de la variation des constantes arbitraires dans tous les problèmes de la mécanique, lu le 19 février 1810 à l'Institut de France, later edition: in uvres de Lagrange. VI, 809 (Paris: Gauthier-Villars, 1877) In the text  
• Newman, W., & Efroimsky, M. 2003, Chaos, 13, 476 In the text  
• Plummer, H. C. 1918, An Introductory Treatise on Dynamical Astronomy (UK: Cambridge University Press) In the text