text:
# ^{𝑎𝑠𝑡𝑟𝑜-𝑝ℎ/0212245} Implicit gauge symmetry emerging in the N-body problem of celestial mechanics.

## 1  Introduction

### 1.1 The Variation-of-Parameters Method

The method of variation of parameters (VOP) was invented in the middle of XVIII-th century by Euler (1748, 1753) as a tool for treating highly-nonlinear problems emerging in the context of celestial mechanics. Though Euler himself eventually lost interest in this approach, his line of research was intensively continued by Lagrange who employed it for deriving his celebrated system of equations describing evolution of the orbital parameters (which in the astronomical parlance are called ”orbital elements”).

In the modern textbooks, the VOP method is normally introduced as one of the means by which one can solve an inhomogeneous linear differential equation: one first finds all solutions of the appropriate linear homogeneous equation, and then instills time dependence into the coefficients in the linear combination of these solutions. Here follows the easiest example:

$`\ddot{x}+p(t)\dot{x}+q(t)x=F(t).`$ (1)
To solve this inhomogeneous equation, one starts out with the homogeneous one:

$$\ddot{x}+p(t)\dot{x}+q(t)x=\mathrm{\hspace{0.33em}0}.$$
(2)
A linear combination of its two fundamental solutions will read: $`Sx_1(t)+Cx_2(t)`$. The recipe has it that at this point one should look for a solution to (1) in ansatz

$$x=S(t)x_1(t)+C(t)x_2(t).$$
(3)
Since the functions $`x_{1,2}(t)`$ are known, what one has to do is just to find $`S(t)`$ and $`C(t)`$. Equation (1) will, by itself, be insufficient for determining two independent functions. The excessive freedom can be removed through a by-hand imposure of an extra equality often chosen as

$$\dot{S}x_1+\dot{C}x_2=\mathrm{\hspace{0.33em}0}.$$
(4)
It greatly simplifies the expressions for $`\dot{x}`$ and $`\ddot{x}`$:

$$\dot{x}=S\dot{x}_1+C\dot{x}_2,\ddot{x}=\dot{S}\dot{x}_1+\dot{C}\dot{x}_2+C\ddot{x}_1+C\ddot{x}_2,$$
(5)
substitution whereof in (1) entails:

$$\dot{S}\dot{x}_1+\dot{C}\dot{x}_2=F.$$
(6)
Together with (4), the latter yields:

$$S=^t\frac{F(t^{})x_2(t^{})}{W[x_1,x_2](t^{})}dt^{},C=^t\frac{F(t^{})x_1(t^{})}{W[x_1,x_2](t^{})}dt^{}$$
(7)
$$\text{where}W[x_1,x_2](t)x_1(t)\dot{x}_2(t)x_2(t)\dot{x}_1(t).$$
(8)
This traditional way of introducing the method of variation of parameters (VOP) is pedagogically flawed because it does not illustrate the full might and generality of this approach. What is important is that the initial equation, whose solution(s) is (are) assumed to be known, does not necessarily need to be linear. Moreover, the parameters to be varied should not necessarily be the coefficients in the linear combination of solutions. Historically, Euler and Lagrange developed this approach in order to solve the nonlinear equation (27), so the parameters to vary (the orbital Keplerian elements) were not coefficients of a linear combination of solutions to the homogeneous equation. Rather, these were the constants of the unperturbed motion (i.e., the quantities which were conserved in the homogeneous (2-body) case but no longer conserved in the inhomogeneous (N-body) case).

A pretty evident aspect of this method is that any choice of the supplementary condition different from (4) will render the same solution $`y(x)`$, even though the expressions for the “constants” $`C_{1,\mathrm{\hspace{0.17em}2}}`$ will differ from (7). In other words, the solution $`y(x)`$ remains invariant under a certain family of transformations of $`C_{1,\mathrm{\hspace{0.17em}2}}`$. This internal symmetry is inherent in any problem that can be approached through the VOP method, including the N-body problem of celestial mechanics where the parameters $`C_i`$ are the constants of the unperturbed motion. To illustrate this, let us consider a trivial example kindly offered to the author by William Newman.

### 1.2 Newman’s Example

As a simplest instance, consider a forced harmonic oscillator

$`H(x,p)={\displaystyle \frac{p^2}{2}}+{\displaystyle \frac{x^2}{2}}+xF(t)`$
that leads to the well known initial-condition problem

$$\ddot{x}+x=F(t),\text{with}x(0)\text{and}\dot{x}(0)\text{known}.$$
(9)
As prescribed by the method of variation of parameters, we seek a solution in ansatz

$$x=S(t)\mathrm{sin}t+C(t)\mathrm{cos}t.$$
(10)
which yields

$$\dot{x}=\left\{\dot{S}(t)\mathrm{sin}t+\dot{C}(t)\mathrm{cos}t\right\}+S(t)\mathrm{cos}tC(t)\mathrm{sin}t.$$
(11)
The standard procedure implies that we put $`\dot{S}(t)\mathrm{sin}t+\dot{C}(t)\mathrm{cos}t=\mathrm{\hspace{0.17em}0}`$, in order to get rid of the ambiguity. The by-hand imposure of this equality is convenient but not necessarily required. Any other way of fixing the ambiguity, like for example,

$$\dot{S}(t)\mathrm{sin}t+\dot{C}(t)\mathrm{cos}t=\mathrm{\Phi }(t)$$
(12)
will be equally good. Then

$$\ddot{x}=\dot{\mathrm{\Phi }}+\dot{S}(t)\mathrm{cos}t\dot{C}(t)\mathrm{sin}tS(t)\mathrm{sin}tC(t)\mathrm{cos}t,$$
(13)
whence

$$\ddot{x}+x=\dot{\mathrm{\Phi }}+\dot{S}(t)\mathrm{cos}t\dot{C}(t)\mathrm{sin}t$$
(14)
Thus one faces the system

$`\dot{\mathrm{\Phi }}+\dot{S}(t)\mathrm{cos}t\dot{C}(t)\mathrm{sin}t=F(t)`$
(15)
$`\dot{S}(t)\mathrm{sin}t+\dot{C}(t)\mathrm{cos}t=\mathrm{\Phi }(t),`$
the first line being the equation of motion (obtained through combining (14) with (9)), and the second line being identity (12).111This identity becomes a constraint if one fixes the functional form of function $`\mathrm{\Phi }`$ which thus far has remained arbitrary. The system trivially resolves to

$`\dot{S}=F\mathrm{cos}t{\displaystyle \frac{d}{dt}}\left(\mathrm{\Phi }\mathrm{cos}t\right)`$
(16)
$`\dot{C}=F\mathrm{sin}t+{\displaystyle \frac{d}{dt}}\left(\mathrm{\Phi }\mathrm{sin}t\right)`$
The function $`\mathrm{\Phi }`$ still remains arbitrary,222Mind that this arbitrariness of $`\mathrm{\Phi }`$ stays, no matter what the initial conditions are to be. Indeed, for fixed $`x(0)`$ and $`\dot{x}(0)`$, the system $`C(0)=x(0),\mathrm{\Phi }(0)+S(0)=\dot{x}(0)`$ solves for $`S(0)`$ and $`C(0)`$ for an arbitrary choice of $`\mathrm{\Phi }(0)`$. as can be easily seen either from the above derivation or from direct substitution of (16) in (9). Integration of (16) trivially yields

$`S={\displaystyle ^t}F\mathrm{cos}t^{}dt^{}\mathrm{\Phi }\mathrm{cos}t+c_1`$
(17)
$`C={\displaystyle ^t}F\mathrm{sin}t^{}dt^{}+\mathrm{\Phi }\mathrm{sin}t+c_2`$
inclusion whereof into (10) entails

$`x=S\mathrm{sin}t+C\mathrm{cos}t=\mathrm{cos}t{\displaystyle ^t}F\mathrm{sin}t^{}dt^{}+\mathrm{sin}t{\displaystyle ^t}F\mathrm{cos}t^{}dt^{}+c_1\mathrm{sin}t+c_2\mathrm{cos}t`$ (18)
with the $`\mathrm{\Phi }`$-terms cancelled out. Analogy of this simple example with the osculating-element formalism is evident. Variable $`x`$ considered in this example is analogous to the Cartesian coordinates and velocities in Section III, while parameters $`S`$ and $`C`$ are counterparts of the osculating elements. Equations (16) are analogues to (66) - (71). These expressions for $`\dot{S}`$ and $`\dot{C}`$ through $`\mathrm{\Phi }`$ show that it is impossible (at least, in this particular example) to pick up such $`\mathrm{\Phi }`$ that the right-hand sides in both equations (16) simplify. Newman’s example also thwarts one’s hope to separate the time scales: as can be seen from (17), even when the perturbation $`F(t)`$ is a slow function of time, the time evolution of the “osculating elements” $`S`$ and $`C`$ is determined by the “fast” time scale associated with $`\mathrm{sin}t`$ and $`\mathrm{cos}t`$, i.e., with the solutions for the homogeneous equation.

The above setting was much simpler than the N-body case where one is faced with nonlinearity and the parameters to be varied are not coefficients of a linear combination of solutions to a homogeneous equations. However, Newman’s example well illustrates the origin of internal symmetry in the VOP context. This symmetry is a generic feature that shows itself in all implementations of the VOP approach, including the N-particle dynamics.

Another good thing about Newman’s example is that it illustrates the impossibility of time-scales’ separation in the VOP method: in the right-hand side of (17) the expressions standing under the integral contain both the free-oscillator and disturbing frequencies. No matter what gauge $`\mathrm{\Phi }`$ one chooses, at least one of these integrals will contain a mixture of frequencies. This way, if we have a fast oscillator subject to a slow force, we still cannot reduce its motion to a slow evolution of the amplitudes $`S(t)`$ and $`C(t)`$; instead, the time dependence of these amplitudes will contain fast terms (which is very counter-intuitive).

## 2  Background

Solar-System dynamics is, largely, variations of the old theme, the gravity law anticipated by Hook and derived from Kepler’s laws by Newton:

$$m_i\ddot{\stackrel{\mathbf{}}{𝝆}}_i=m_i\left\{\underset{ji}{}Gm_j\frac{\stackrel{\mathbf{}}{𝒓}_{ij}}{r_{ij}^3}\right\}=m_i\frac{U_i}{\stackrel{\mathbf{}}{𝝆}_i},\stackrel{\mathbf{}}{𝒓}_{ij}\stackrel{\mathbf{}}{𝝆}_j\stackrel{\mathbf{}}{𝝆}_i,i,j=1,\mathrm{},N,$$
(19)
$`m_i`$ and $`\stackrel{\mathbf{}}{𝝆}_i`$ being the masses and their inertial-frame-related positions, $`U_i`$ being the overall potential acting on $`m_i`$:

$$U_i\underset{ji}{}Gm_j\frac{1}{\rho _{ij}},$$
(20)
the sign convention chosen as in the astronomical, not as in the physical literature. The equations of motion may be conveniently reformulated in terms of the relative locations

$$\stackrel{\mathbf{}}{𝒓}_i\stackrel{\mathbf{}}{𝒓}_{is}\stackrel{\mathbf{}}{𝝆}_i\stackrel{\mathbf{}}{𝝆}_s,$$
(21)
$`\stackrel{\mathbf{}}{𝝆}_s`$ standing for the position of Sun. The difference between

$$\ddot{\stackrel{\mathbf{}}{𝝆}}_i=\underset{ji,s}{}G\frac{m_j\stackrel{\mathbf{}}{𝒓}_{ij}}{r_{ij}^3}+G\frac{m_s\stackrel{\mathbf{}}{𝒓}_{is}}{r_{is}^3}$$
(22)
and

$$\ddot{\stackrel{\mathbf{}}{𝝆}}_s=\underset{ji,s}{}G\frac{m_j\stackrel{\mathbf{}}{𝒓}_{sj}}{r_{sj}^3}+G\frac{m_i\stackrel{\mathbf{}}{𝒓}_{si}}{r_{si}^3}$$
(23)
amounts to:

$$\ddot{\stackrel{\mathbf{}}{𝒓}}_i=\underset{ji,s}{}G\frac{m_j\stackrel{\mathbf{}}{𝒓}_{ij}}{r_{ij}^3}\underset{ji,s}{}G\frac{m_j\stackrel{\mathbf{}}{𝒓}_j}{r_j^3}G\frac{\left(m_i+m_s\right)\stackrel{\mathbf{}}{𝒓}_i}{r_i^3}=\frac{\stackrel{~}{U}_i}{\stackrel{\mathbf{}}{𝒓}_i}$$
(24)
$`\stackrel{~}{U}`$ being the new potential:

$$\stackrel{~}{U}_i\frac{G\left(m_i+m_s\right)}{r_i}+R_i,$$
(25)
with the disturbing function

$$R_i\underset{ji}{}Gm_j\left\{\frac{1}{r_{ij}^3}\frac{\stackrel{\mathbf{}}{𝝆}_i\stackrel{\mathbf{}}{𝝆}_j}{\rho _j^3}\right\}$$
(26)
singled out.

Formulae (24) - (25) become trivial in the case of the two-body problem where only $`m_i`$ and $`m_s`$ are present. In this situation the disturbing function vanishes and the motion is, mathematically, equivalent to rotation about a nailed-down body of mass $`m_i+m_s`$ located at some fixed point $`O`$:

$$\ddot{\stackrel{\mathbf{}}{𝒓}}+\frac{\mu }{r^2}\frac{\stackrel{\mathbf{}}{𝒓}}{r}=\mathrm{\hspace{0.33em}0},\stackrel{\mathbf{}}{𝒓}\stackrel{\mathbf{}}{𝝆}_{planet}\stackrel{\mathbf{}}{𝝆}_s,\mu G(m+m_s).$$
(27)
In here $`\stackrel{\mathbf{}}{𝒓}\stackrel{\mathbf{}}{𝒓}_1\stackrel{\mathbf{}}{𝒓}_i`$, because the subscript $`i`$ runs through one value solely: $`i=\mathrm{\hspace{0.17em}1}`$.

This setting permits exact analytical treatment that leads to the famous Newtonian result: the orbit is elliptic and has the gravitating centre in one of its foci. This enables a transition from Cartesian to Keplerian coordinates. For our further study this transition will be very important, so we shall recall it in detail.

At any instant of time, the position $`𝒓`$ and velocity $`\dot{𝒓}`$ of an orbiting body can be determined by its coordinates $`(x,y,z)`$ and derivatives $`(\dot{x},\dot{y},\dot{z})`$ in an inertial frame with origin located in point $`O`$ where the mass $`m_i+m_s`$ rests. The position of the orbital ellipse may be fully defined by the longitude of the node, $`\mathrm{\Omega }`$; the inclination, $`i`$; and the argument of pericentre, $`\omega `$ (instead of the latter, one can introduce the longitude of pericentre, $`\stackrel{~}{\omega }\mathrm{\Omega }+\omega `$). The shape of the ellipse is parametrised by its eccentricity, $`e`$, and semimajor axis, $`a`$. Position of a point on the ellipse may be charachterised, for example, by the eccentric anomaly, $`E`$. As well known,

$$Ee\mathrm{sin}E=ntB,$$
(28)
$`B`$ being a constant of integration, and $`n`$ being the mean motion defined as

$$n\mu ^{1/2}a^{3/2}.$$
(29)
One can then introduce, following Kepler, the mean anomaly, $`M`$ as

$$MEe\mathrm{sin}E.$$
(30)
Let $`t_o`$ be the fiducial time. Then, by putting $`B=M_o+nt_o`$, we can introduce, instead of $`B`$, another integration constant, $`M_o`$. Hence, (28) will read:

$$M=M_o+n\left(tt_o\right),$$
(31)
the meaning of $`M_o`$ being self-evident: it is the value of $`M`$ at the reference epoch $`t_o`$. So introduced the mean anomaly provides another parameterisation of the position of a planet on the ellipse. In the disturbed case the latter formula naturally becomes

$$M=M_o+_{t_o}^tn(t^{})𝑑t^{}.$$
(32)
One more convenient parameter often employed in the literature is the mean longitude $`\lambda `$ defined by

$$\lambda \stackrel{~}{\omega }+M=\mathrm{\Omega }+\omega +M.$$
(33)
Unless the inclination is zero, neither the longitude of the pericentre, $`\stackrel{~}{\omega }`$, nor the mean longitude, $`\lambda `$, is a true angle. They are sums of angles in two different planes that meet at the node.

Planetary dynamics is based on application of the above, 2-body, formalism to the N-body case. Naively speaking, since the mutual disturbances of planets are very weak compared to the solar gravity, it seems natural to assume that the planets move along ellipses which are slowly evolving. Still, the weakness of perturbations is, by itself, a very shaky foundation for the varying-ellipse method. This so physically-evident circumstance has a good illustrative power but is of no help when the following questions arise:

(1) To what degree of rigour can an orbit curve be modelled by a family of instantaneous ellipses having the Sun in one of their foci? Can this be performed exactly?

(2) Is this representation of the curve by a family of ellipses unique?

These two questions will not seem anecdotal, if we recall that the concept of evolving instantaneous ellipses had been introduced into practice (and that major developments of the disturbing-function theory had been accomplished) long before Frenet and Serait developed the theory of curves.333I am grateful to William Newman who drew my attention to this circumstance. (This historical paradox explains the reason why words ”helicity” and ”torsion” are still absent from the astronomy textbooks.)

Fortunately, Lagrange, who authored the idea of instantaneous ellipses, fortified it with such powerful tools of calculus, that in this case they surpassed the theory of curves. Moreover, these tools in no way relied on the weakness of the disturbances. Hence, Lagrange’s treatment of the problem already contained an affirmative answer to the first question.

Below we shall demonstrate that the answer to the second question is negative. Moreover, it turns out that the question calls into being a rich, though not new, mathematical structure. We shall show that the Lagrange system of equations for the instantaneous orbital elements possesses a hidden symmetry not visible with the naked eye. This symmetry is very similar to the gauge symmetry, one well known from electrodynamics. A careful analysis shows that the Lagrange system, as we know it, is written in some specific gauge: all trajectories constrained to some 9-dimensional submanifold in the 12-dimensional space constituted by the Keplerian elements and their time derivatives.

Beside the possible practical relevance to orbit computation, the said symmetry unveils a fiber bundle structure hidden behind Lagrange’s system of equations for the Keplerian elements. The symmetry is absent in the 2-body case, but comes into being in the N-body setting ($`N\mathrm{\hspace{0.17em}3}`$) where each orbiting body follows a ellipse of varying shape, but the time evolution of the ellipse contains an inherent ambiguity.

Here follows a crude illustration of this point. Imagine two coplanar ellipses sharing one focus. Let one ellipse slowly rotate within its plane, about the shared focus. Let the other ellipse rotate faster, also in its plane, in the same direction, and about that same shared focus. Suppose a planet is at one of the points of these ellipses’ intersection. One observer may state that the planet is rapidly moving along a slowly rotating ellipse, while another observer may insist that the planet is slowly describing the fast-moving ellipse. Both descriptions will be equally legitimate, for there exists an infinite number of ways of dividing the actual motion of the planet into its motion along some orbit and simultaneous evolution of the orbit itself. Needless to say, the real, physical trajectory is unique. However, its description (parameterisation 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 on that point. However, in his treatment he passingly introduced a convenient mathematical condition similar to (4), which removed the said ambiguity. This condition and possible alternatives to it will be the topic of the further sections.

## 3 Keplerian coordinates in the 2-body and N-body problems: Osculating Elements vs Orbital Elements

Although not widely recognised, the perturbation equations of celestial mechanics possess a gauge freedom. It is probable that this was already noticed by Euler and Lagrange in the middle of the XVIII century. However, although the existence of this freedom did not entirely escape attention,444 In 2002 a veteran of the Russian astronomy Yurii Batrakov mentioned to the author that back in his student years his lecturer at the St.Petersburg (then Leningrad) State University, Mihail Subbotin, emphasised in his course the possibility of making alternative choices. Interestingly, though, this fact was not developed in the published version of that lecture course (Subbotin 1968) where it is only passingly mentioned that the Lagrange constraint is imposed merely for convenience. its consequences have yet to be fully explored.

Perhaps the easiest way to gain an appreciation of this freedom is to follow the derivation of the perturbation equations by application of the variation of parameters (VOP) technique as invented by Euler and Lagrange and shaped into its final form in Lagrange (1808, 1809, 1810). Lagrange put the planetary equations in a closed form in which temporal derivatives of the orbital elements were expressed in terms of partial derivatives of the disturbing function with respect to the orbital elements. A closely related development is presented in the textbook by Brouwer and Clemence (1961). We shall start in the spirit of Lagrange but will soon deviate from it in two points. First, in this subsection we shall not assume that the disturbing force is conservative and that it depends upon the positions solely, but shall permit it to depend also upon velocities. Second, neither in this subsection nor further shall we impose the Lagrange constraint.

Before addressing the N-particle case, Lagrange referred to the reduced 2-body problem,

$`\ddot{\stackrel{}{𝐫}}+{\displaystyle \frac{\mu }{r^2}}{\displaystyle \frac{\stackrel{}{𝐫}}{r}}=\mathrm{\hspace{0.33em}0},`$
(34)
$`\stackrel{\mathbf{}}{𝒓}\stackrel{\mathbf{}}{𝒓}_{planet}\stackrel{\mathbf{}}{𝒓}_{sun},\mu G(m_{planet}+m_{sun}).`$
whose generic solution, a Keplerian ellipse or a hyperbola, can be expressed, in some fixed Cartesian frame, as

$`x=f_1(C_1,\mathrm{},C_6,t),\dot{x}=g_1(C_1,\mathrm{},C_6,t),`$
$`y=f_2(C_1,\mathrm{},C_6,t),\dot{y}=g_2(C_1,\mathrm{},C_6,t),`$ (35)
$`z=f_3(C_1,\mathrm{},C_6,t),\dot{z}=g_3(C_1,\mathrm{},C_6,t),`$
or, shortly:

$$\stackrel{\mathbf{}}{𝒓}=\stackrel{}{𝐟}(C_1,\mathrm{},C_6,t),\dot{\stackrel{\mathbf{}}{𝒓}}=\stackrel{}{𝐠}(C_1,\mathrm{},C_6,t),$$
(36)
the functions $`g_i`$ being, by definition, partial derivatives of $`f_i`$ with respect to the last argument:

$$\stackrel{}{𝐠}\left(\frac{\stackrel{}{𝐟}}{t}\right)_{C=const}.$$
(37)
Naturally, the general solution contains six adjustable constants, $`C_i`$, since the problem (34) is constituted by three, second order, differential equations. To find the explicit form of the dependence (36), one can employ an auxiliary set of Cartesian coordinates $`\stackrel{}{𝐪}`$, with an origin at the gravitating centre, and with the first two axes located in the plane of orbit. In terms of the true anomaly $`f`$, these coordinates will read:

$`q_1r\mathrm{cos}f,q_2r\mathrm{sin}f,q_3=\mathrm{\hspace{0.33em}0}.`$ (38)
In the two-body case, their time derivatives can be easily computed and expressed through the major semiaxis, $`a`$, the eccentricity, $`e`$, and the true anomaly, $`f`$:

$`\dot{q}_1={\displaystyle \frac{na\mathrm{sin}f}{\sqrt{1e^2}}},\dot{q}_2={\displaystyle \frac{na(e+\mathrm{cos}f)}{\sqrt{1e^2}}},\dot{q}_3=\mathrm{\hspace{0.33em}0}.`$ (39)
(true anomaly $`f`$ itself being a function of $`a`$, $`e`$, and of the mean anomaly $`MM_o+_{t_o}^tn𝑑t`$, $`n\mu ^{1/2}a^{3/2}`$). In the two-body setting, the inertial-frame-related position and velocity will appear as:

$`\stackrel{\mathbf{}}{𝒓}=𝐑(\mathrm{\Omega },i,\omega )\stackrel{}{𝐪}(a,e,M_o,t),`$
(40)
$`\dot{\stackrel{\mathbf{}}{𝒓}}=𝐑(\mathrm{\Omega },i,\omega )\dot{\stackrel{}{𝐪}}(a,e,M_o,t),`$
$`𝐑(\mathrm{\Omega },i,\omega )`$ being the matrix of rotation from the orbital-plane-related axes $`𝐪`$ to the fiducial frame $`(x,y,z)`$ in which the vector $`\stackrel{}{𝐫}`$ is defined. The rotation is parametrised by the three Euler angles: inclination, $`i`$; the longitude of the node, $`\mathrm{\Omega }`$; and the argument of the pericentre, $`\omega `$.

This is one possible form of the general solution (36). It has been obtained under the convention that a particular ellipse is parametrised by the Lagrange set of orbital elements, $`C_i(e,a,M_o,\omega ,\mathrm{\Omega },i)`$. A different functional form of the same solution is achieved in terms of the Delaunay set, $`D_i(L,G,H,\omega ,\mathrm{\Omega },M_o)`$. Still another possibility is to express the general solution through the initial conditions: then the constants $`(x_o,y_o,z_o,\dot{x}_o,\dot{y}_o,\dot{z}_o)`$ are the six parameters defining a particular orbit. The latter option is natural when the integration is carried out in Cartesian components, but is impractical otherwise.

At this point it is irrelevant which particular set of the adjustable parameters is employed. Hence we shall leave, for a while, the solution in its most general form (35 \- 37) and shall, following Lagrange (1808, 1809, 1810), employ it as an ansatz for a solution of the N-particle problem where the disturbing force acting at a particle is denoted by $`\mathrm{\Delta }\stackrel{}{𝐅}`$:555Our treatment covers disturbing forces $`\mathrm{\Delta }\stackrel{}{𝐅}`$ that are arbitrary vector-valued functions of position and velocity.

$$\ddot{\stackrel{}{𝐫}}+\frac{\mu }{r^2}\frac{\stackrel{}{𝐫}}{r}=\mathrm{\Delta }\stackrel{}{𝐅},$$
(41)
the ”constants” now being time dependent:

$$\stackrel{}{𝐫}=\stackrel{}{𝐟}(C_1(t),\mathrm{},C_6(t),t),$$
(42)
and the functional form of $`\stackrel{}{𝐟}`$ remaining the same as in (36). Now the velocity

$$\frac{d\stackrel{\mathbf{}}{𝒓}}{dt}=\frac{\stackrel{}{𝐟}}{t}+\underset{i}{}\frac{\stackrel{}{𝐟}}{C_i}\frac{dC_i}{dt}=\stackrel{}{𝐠}+\underset{i}{}\frac{\stackrel{}{𝐟}}{C_i}\frac{dC_i}{dt},$$
(43)
will contain a new input, $`(\stackrel{}{𝐟}/C_i)(dC_i/dt)`$, while the first term, $`\stackrel{}{𝐠}`$, will retain the same functional form as it used to have before.

Substitution of $`\stackrel{}{𝐟}(C_1(t),\mathrm{},C_6(t),t)`$ into the perturbed equation of motion (41) gives birth to three independent scalar differential equations of the second order. These three equations contain one independent parameter, time, and six time-dependent variables $`C_i(t)`$ whose evolution is to be determined. Evidently, this cannot be done in a single way because the number of variables exceeds, by three, the number of equations. This means that, though the  physical  orbit (comprised by the locus of points in the Cartesian frame and by the values of velocity in each of these points) is unique, its parameterisation in terms of the orbital elements is ambiguous. Lagrange, in his treatment, noticed that the system was underdefined, and decided to amend it with exactly three independent conditions which would make it solvable.

Before moving on with the algebra, let us look into the mathematical nature of this ambiguity. A fixed Keplerian ellipse (36), which is the solution to the 2-body problem (34), gives birth to a time-dependent one-to-one (within one revolution period) mapping

$$(C_1,\mathrm{},C_6)(x(t),y(t),z(t),\dot{x}(t),\dot{y}(t),\dot{z}(t)).$$
(44)
In the N-body case, the new ansatz (42) is incompatible with (36). This happens because now the time derivatives of coordinates $`C_i`$ come into play in (43). Hence, instead of (44), one gets a time-dependent mapping between a 12-dimensional and a 6-dimensional spaces:

$$(C_1(t),\mathrm{},C_6(t),\dot{C}_1(t),\mathrm{},\dot{C}_6(t))(x(t),y(t),z(t),\dot{x}(t),\dot{y}(t),\dot{z}(t)).$$
(45)
This brings up two new issues. One is the multiple scales: while the physical motion along an instantaneous ellipse parametrised by $`C_i`$ is associated with the ”fast time scale”, the evolution of the osculating elements $`C_i(t)`$ represents the ”slow time scale”. The quotation marks are necessary here, because in reality these time scales are inseparably connected, so that the ”slow time scale” is not always slower than the ”fast time scale”. What is important here is that, in general, ansatz (42) gives birth to two separate time scales. 666In practice, the mean longitude $`\lambda =\lambda _o+_{t_o}^tn(t)𝑑t`$ is often used instead of its fiducial-epoch value $`\lambda _o`$. Similarly, those authors who prefer the mean anomaly to the mean longitude, often use $`M=M_o+_{t_o}^tn(t)𝑑t`$ instead of $`M_o`$. While $`M_o`$ and $`\lambda _o`$ are orbital elements, the quantities $`M`$ and $`\lambda `$ are not. Still, the time-dependent change of variables from $`\lambda _o`$ to $`\lambda `$ (or from $`M_o`$ to $`M`$) is perfectly legitimate. Being manifestly time-dependent, this change of variables intertwines two different time scales: for example, $`M`$ carries a “fast” time dependence through the upper limit of the integral in $`M=M_o+_{t_o}^tn(t)𝑑t`$, and it also carries a “slow” time-dependence due to the adiabatic evolution of the osculating element $`M_o`$. The same concerns $`\lambda `$. The second important issue is that mapping (45) cannot be one-to-one. As already mentioned, the three equations of motion (41) are insufficient for determination of six functions $`C_1,\mathrm{}C_6`$ and, therefore, one has a freedom to impose, by hand, three extra constraints upon these functions and their derivatives777A more accurate mathematical discussion of this freedom should be as follows. The dynamics, in the form of first-order differential equations for the orbital coordinates $`C_i(t)`$ and their derivatives $`H_i(t)\dot{\mathrm{C}}_i(t)`$, will include six evident first-order identities for these twelve functions: $`H_i(t)=dC_i(t)/dt`$. Three more differential equations will be obtained by plugging $`\stackrel{\mathbf{}}{𝒓}=\stackrel{\mathbf{}}{𝒇}(C_1,\mathrm{},C_6,t)`$ into (41). These equations will be of the second order in $`C_i(t)`$. However, in terms of both $`C_i(t)`$ and $`H_i(t)`$ these equations will be of the first order only. Altogether, we have nine first-order equations for twelve functions $`C_i(t)`$ and $`H_i(t)`$. Hence, the problem is underdefined and permits three extra conditions to be imposed by hand. The arbitrariness of these conditions reveals the ambiguity of the representation of an orbit by instantaneous Keplerian ellipses. Mappings between different representations reveal an internal symmetry (and a symmetry group) underlying this formalism..

Though Lagrange did notice that the system was underdefined, he never elaborated on the underlying symmetry. He simply imposed three convenient extra conditions

$$\underset{i}{}\frac{\stackrel{}{𝐟}}{C_i}\frac{dC_i}{dt}=\mathrm{\hspace{0.33em}0},$$
(46)
and went on, to derive (in this particular gauge, which is often called ”Lagrange constraint”) his celebrated system of equations for orbital elements. Now we can only speculate on why Lagrange did not bother to explore this ambiguity and its consequences. One possible explanation is that he did not have the concept of continuous groups and symmetries in his arsenal (though it is very probable that he knew the concept of discrete group888In his paper on solution, in radicals, of equations of degrees up to four, Réflexions sur la résolution algébrique des équations, dated by 1770, Lagrange performed permutations of roots. Even though he did not consider compositions of permutations, his technique reveals that, most likely, he was aware of, at least, the basic idea of discrete groups and symmetries.). Another possibility is that Lagrange did not expect that exploration of this ambiguity would reveal any promising tools for astronomical calculations.

Lagrange’s choice of the supplementary constraints was motivated by both physical considerations and the desire to simplify calculations. Since, physically, the time-dependent set $`(C_1(t),\mathrm{},C_6(t))`$ can be interpreted as an instantaneous ellipse, in a bound-orbit case, or an instantaneous hyperbola, in a fly-by situation, Lagrange decided to make the instantaneous orbital elements $`C_i`$ osculating, i.e., he postulated that the instantaneous ellipse (or hyperbola) must always be tangential to the physical trajectory. This means that the physical trajectory defined by $`(C_1(t),\mathrm{},C_6(t))`$ must, at each instant of time, coincide with the unperturbed (two-body) orbit that the body would follow if perturbations were to cease instantaneously. This can be achieved only in case the velocities depend upon the elements, in the N-body problem, in the same manner as they did in the 2-body case. This, in turn, can be true only if one asserts that the extra condition (46) is fulfilled. That condition, also called Lagrange constraint, consists of three scalar equations which, together with the three equations of motion (41), constitute a well-defined system of six equations with six variables $`(C_1(t),\mathrm{},C_6(t))`$.

For the reasons explained above, this constraint, though well motivated, remains, from the mathematical viewpoint, completely arbitrary: it considerably simplifies the calculations but does not influence the shape of the physical trajectory and the rate of motion along that curve. One could as well choose some different supplementary condition

$$\underset{i}{}\frac{\stackrel{}{𝐟}}{C_i}\frac{dC_i}{dt}=\stackrel{}{𝚽}(C_{1,\mathrm{},6},t),$$
(47)
$`\stackrel{\mathbf{}}{𝚽}`$ now being an arbitrary function of time and parameters $`C_i`$.999In principle, one may endow $`\stackrel{\mathbf{}}{𝚽}`$ also with dependence upon the parameters’ time derivatives of all orders. This would yield higher-than-first-order time derivatives of the $`C_i`$ in subsequent developments requiring additional initial conditions, beyond those on $`\stackrel{}{𝐫}`$ and $`\dot{\stackrel{}{𝐫}}`$, to be specified in order to close the system. We avoid this unnecessary complication by restricting $`\stackrel{\mathbf{}}{𝚽}`$ to be a function of time and the $`C_i`$. Substitution of (46) by (47) would leave the physical motion unchanged but would alter the subsequent mathematics and, most importantly, would eventually yield different solutions for the orbital elements. Such invariance of a physical theory under a change of parameterisation is an example of gauge symmetry (Efroimsky 2002). This freedom is a particular case of a more general type of invariance that emerges in applications of the VOP method to a vast variety of ODE problems.

The importance of this gauge freedom is determined by two circumstances which parallel similar circumstances in electrodynamics. One consequence of the gauge invariance is its non-conservation in the course of orbit computation. This, purely numerical, phenomenon may be called gauge drift. Another relevant issue is that a clever choice of gauge often simplifies the solution of the equations of motion. In application to the theory of orbits, this means that a deliberate choice of non-osculating orbital elements (i.e., of a set $`C_i`$ obeying some condition (47) different from (46)) can sometimes simplify the equations for these elements’ evolution. As explained in Efroimsky & Goldreich (2003b), some earlier developments by Goldreich (1965), Brumberg et al (1971) can be interpreted as calculations in non-Lagrange gauges.

The standard derivation of both Delaunay and Lagrange systems of planetary equations by the VOP method rests on the assumption that the Lagrange constraint (46) is fulfilled. Both systems get altered under a different gauge choice. Derivation of the gauge-invariant (i.e., valid in an arbitrary gauge (47)) Lagrange system can be found in Efroimsky (2002). Derivation of the gauge-invariant Delaunay system is given in Efroimsky (2002). In Efroimsky & Goldreich (2003b) both systems are amended with extra terms which emerge when the disturbance happens to depend not only upon the positions but also upon velocities. This situation is encountered when one is working in a non-inertial coordinate frame or when relativistic corrections to the Newtonian force are accounted for.

Without going into excessive details, that may be found in (Efroimsky 2002), we give the essence of that derivation so that we can refer to it in one of the subsequent sections.

¿From (43), it follows that the formula for the acceleration reads:

$$\frac{d^2\stackrel{\mathbf{}}{𝒓}}{dt^2}=\frac{\stackrel{}{𝐠}}{t}+\underset{i}{}\frac{\stackrel{}{𝐠}}{C_i}\frac{dC_i}{dt}+\frac{d\stackrel{}{𝚽}}{dt}=\frac{^2\stackrel{}{𝐟}}{^2t}+\underset{i}{}\frac{\stackrel{}{𝐠}}{C_i}\frac{dC_i}{dt}+\frac{d\stackrel{}{𝚽}}{dt},$$
(48)
Together with the equation of motion (41), it leads to:

$$\frac{^2\stackrel{}{𝐟}}{t^2}+\frac{\mu }{r^2}\frac{\stackrel{}{𝐟}}{r}+\underset{i}{}\frac{\stackrel{}{𝐠}}{C_i}\frac{dC_i}{dt}+\frac{d\stackrel{}{𝚽}}{dt}=\mathrm{\Delta }\stackrel{}{𝐅},r|\stackrel{\mathbf{}}{𝒓}|=|\stackrel{}{𝐟}|.$$
(49)
As $`\stackrel{}{𝐟}`$ is, by definition, a Keplerian solution to the two-body problem (and, therefore, obeys the unperturbed equation (34) ), the above formula becomes:

$$\underset{i}{}\frac{\stackrel{}{𝐠}}{C_i}\frac{dC_i}{dt}=\mathrm{\Delta }\stackrel{}{𝐅}\frac{d\stackrel{}{𝚽}}{dt}.$$
(50)
This is the equation of disturbed motion, written in terms of the orbital elements. Together with constraint (47) it constitutes a well-defined system of equations that can be solved with respect to $`dC_i/dt`$ for all $`i`$’s. The easiest way of doing this is to employ the elegant technique offered by Lagrange: to multiply the equation of motion by $`\stackrel{}{𝐟}/C_n`$ and to multiply the constraint by $`\stackrel{}{𝐠}/C_n`$. These operations yield the following equalities

$$\frac{\stackrel{}{𝐟}}{C_n}\left(\underset{j}{}\frac{\stackrel{}{𝐠}}{C_j}\frac{dC_j}{dt}\right)=\frac{\stackrel{}{𝐟}}{C_n}\mathrm{\Delta }\stackrel{}{𝐅}\frac{\stackrel{}{𝐟}}{C_n}\frac{d\stackrel{}{𝚽}}{dt}$$
(51)
and

$$\frac{\stackrel{}{𝐠}}{C_n}\left(\underset{j}{}\frac{\stackrel{}{𝐟}}{C_j}\frac{dC_j}{dt}\right)=\frac{\stackrel{}{𝐠}}{C_n}\stackrel{}{𝚽},$$
(52)
summation whereof results in:

$$\underset{j}{}[C_nC_j]\frac{dC_j}{dt}=\frac{\stackrel{}{𝐟}}{C_n}\mathrm{\Delta }\stackrel{}{𝐅}\frac{\stackrel{}{𝐟}}{C_n}\frac{d\stackrel{}{𝚽}}{dt}\frac{\stackrel{}{𝐠}}{C_n}\stackrel{}{𝚽},$$
(53)
where the symbol $`[C_nC_j]`$ denotes the unperturbed (i.e., defined as in the two-body case) Lagrange brackets:

$$[C_nC_j]\frac{\stackrel{}{𝐟}}{C_n}\frac{\stackrel{}{𝐠}}{C_j}\frac{\stackrel{}{𝐟}}{C_j}\frac{\stackrel{}{𝐠}}{C_n}.$$
(54)
Above we specified that $`\stackrel{\mathbf{}}{𝚽}`$ is a function of time and the parameters $`C_i`$, but not of their derivatives. Under this restriction, (53) may be conveniently rewritten as

$$\underset{j}{}\left([C_nC_j]+\frac{\stackrel{}{𝐟}}{C_n}\frac{\stackrel{}{𝚽}}{C_j}\right)\frac{dC_j}{dt}=\frac{\stackrel{}{𝐟}}{C_n}\mathrm{\Delta }\stackrel{}{𝐅}\frac{\stackrel{}{𝐟}}{C_n}\frac{\stackrel{}{𝚽}}{t}\frac{\stackrel{}{𝐠}}{C_n}\stackrel{}{𝚽}.$$
(55)
This expresses the most generic form of the gauge-invariant perturbation equations of celestial mechanics, that follows from the VOP method.101010 I am grateful to Peter Goldreich who recommended me to present the equations in this concise and so general form.

We can obtain an immediate solution for the individual $`dC_i/dt`$ in the Lagrange gauge by exploiting the well known expression for the inverse Lagrange-bracket, or Poisson-bracket, matrix that is offered in the literature for the two-body problem. The presence of the term proportional to $`\stackrel{}{\mathrm{\Phi }}/C_j`$ on the lhs of (65) complicates the solution for the individual $`dC_i/dt`$ in an arbitrary gauge, but only to the extent of requiring the resolution of a set of six, simulataneous, linear, algebraic equations.

All different choices of three (compatible and sufficient) gauge conditions expressed by the vector $`\stackrel{\mathbf{}}{𝚽}`$ will lead to physically equivalent results. This equivalence means the following. Suppose we solve the equations of motion for $`C_{1,\mathrm{},6}`$, with some gauge condition $`\stackrel{}{𝚽}`$ enforced. This will give us the solution, $`C_{1,\mathrm{},6}(t)`$. If, though, we choose to integrate the equations of motion with another gauge $`\stackrel{~}{\stackrel{\mathbf{}}{𝚽}}`$ enforced, then we shall arrive at a solution $`\stackrel{~}{C}_i(t)`$ of a different functional form. Stated alternatively, in the first case the integration in the 12-dimensional space $`(C_{1,\mathrm{},6},H_{1,\mathrm{},6})`$ will be restricted to one 9-dimensional time-dependent submanifold, whereas in the second case it will be restricted to another submanifold. Despite this, both solutions, $`C_i(t)`$ and $`\stackrel{~}{C}_i(t)`$, will give, when substituted back in (42), the same orbit $`(x(t),y(t),z(t))`$ with the same velocities $`(\dot{x}(t),\dot{y}(t),\dot{z}(t))`$. This is a fiber-bundle-type structure, and it gives birth to a 1-to-1 map of $`C_i(t)`$ onto $`\stackrel{~}{C}_i(t)`$. This map is merely a reparameterisation. In physicists’ parlance it will be called gauge transformation. All such reparameterisations constitute a group of symmetry, which would be called, by a physicist, gauge group. The real orbit is invariant under the reparameterisations which are permitted by the ambiguity of gauge-condition choice. This physical invariance implements itself, technically, as form-invariance of the expression (42) under the afore mentioned map. This is analogous to Maxwell’s electrodynamics: the components $`x,y,`$ and $`z`$ of vector $`\stackrel{\mathbf{}}{𝒓}`$ , and their time derivatives, play the role of the physical fields $`\stackrel{}{𝐄}`$ and $`\stackrel{}{𝐁}`$, while the Keplerian coordinates $`C_1,\mathrm{},C_6`$ play the role of the 4-potential $`A^\mu `$. This analogy can go even further.111111Suppose one is solving a problem of electromagnetic wave proliferation, in terms of the 4-potential $`A^\mu `$ in some fixed gauge. An analytic calculation will render the solution in that same gauge, while a numerical computation will furnish the solution in a slightly different gauge. This will happen because of numerical errors’ accumulation. In other words, numerical integration will slightly deviate from the chosen submanifold. This effect may become especially noticeable in long-term orbit computations. Another relevant topic emerging in this context is comparison of two different solutions of the N-body problem: just as in the field theory, in order to compare solutions, it is necessary to make sure if they are written down in the same gauge. Otherwise, the difference between them may be, to some extent, not of a physical but merely of a gauge nature.

## 4 The hidden symmetry of the Lagrange system

If we impose, following Lagrange, the gauge condition (46), then the equation of motion (65) will simplify:

$$\underset{j}{}[C_nC_j]\frac{dC_j}{dt}=\frac{\stackrel{}{𝐟}}{C_n}\mathrm{\Delta }\stackrel{}{𝐅}.$$
(56)
In assumption of the disturbing force being dependent only upon positions and being conservative,121212 I am grateful to Peter Goldreich who drew my attention to the fact that this assertion is right in inertial axes solely. Consider an observer who associates himself with a non-inertial system of reference (say, with a frame fixed on a precessing planet) and defines the orbital elements in this system. He will then have to take into account the non-inertial contribution to the disturbing function. This contribution, denoted in (Goldreich 1965) by $`R_I`$, will depend not only upon $`\stackrel{}{𝐫}`$ but also upon $`\dot{\stackrel{}{𝐫}}`$. In this situation, the terms $`R/C_r`$ in the Lagrange equations should be substituted by $`R/C_r(\dot{\stackrel{}{𝐫}}/C_r)(R/\dot{\stackrel{}{𝐫}})`$. The gauge approach to the problem of a satellite orbiting a wobbling planet will be addressed in one of our subsequent articles (Efroimsky & Goldreich 2003b). In regard to the velocity-dependence of the disturbing function, we would also mention that a similar situation emerges in the relativistic mechanics, because the relativistic correction to the Newton law of gravity bears dependence upon the velocity. we may substitute in the above equation the disturbing force by the gradient of a (position-dependent) disturbing potential:

$$\mathrm{\Delta }\stackrel{}{𝐅}=\frac{R(\stackrel{\mathbf{}}{𝒓})}{\stackrel{\mathbf{}}{𝒓}},$$
(57)
which will result in:

$$\underset{j}{}[C_nC_j]\frac{dC_j}{dt}=\frac{R}{C_n}.$$
(58)
Expressions for the Lagrange brackets are known (Brouwer & Clemence 1961, p. 284), and their insertion into (58) equation will entail the well-known Lagrange system of planetary equation:

$$\frac{da}{dt}=\frac{2}{na}\frac{R}{M_o}$$
(59)
$$\frac{de}{dt}=\frac{1e^2}{na^2e}\frac{R}{M_o}\frac{(1e^2)^{1/2}}{na^2e}\frac{R}{a}$$
(60)
$$\frac{d\omega }{dt}=\frac{\mathrm{cos}i}{na^2(1e^2)^{1/2}\mathrm{sin}i}\frac{R}{i}+\frac{(1e^2)^{1/2}}{na^2e}\frac{R}{e}$$
(61)
$$\frac{di}{dt}=\frac{\mathrm{cos}i}{na^2(1e^2)^{1/2}\mathrm{sin}i}\frac{R}{\omega }\frac{1}{na^2(1e^2)^{1/2}\mathrm{sin}i}\frac{R}{\mathrm{\Omega }}$$
(62)
$$\frac{d\mathrm{\Omega }}{dt}=\frac{1}{na^2(1e^2)^{1/2}\mathrm{sin}i}\frac{R}{i}$$
(63)
$$\frac{dM_o}{dt}=\frac{1e^2}{na^2e}\frac{R}{e}\frac{2}{na}\frac{R}{a}$$
(64)
If analytical integration of this system were possible, it would render a correct orbit, in the fixed gauge (46). A numerical integrator, however, may cause drift from the chosen submanifold (46). Even if the drift is not steady, some deviation from the submanifold is unavoidable. One may wish $`𝚽`$ to be as close to zero as possible, but in reality $`𝚽`$ will be some unknown function whose proximity to zero will be determined by the processor’s error and by the number of integration steps. Even if we begin with (46) fulfilled exactly, the very first steps will give us such values of $`C_i`$ that, being substituted into the lhs of (46), they will give some new value of $`𝚽`$ slightly different from zero. Thus, the Lagrange gauge will no longer be observed.

In case we relax the Lagrange constraint and accept an arbitrary gauge (47) then, under the simplifying assertion (57), equation (65) will become

$$\underset{j}{}\left([C_nC_j]+\frac{\stackrel{}{𝐟}}{C_n}\frac{\stackrel{}{𝚽}}{C_j}\right)\frac{dC_j}{dt}=\frac{R}{C_n}\frac{\stackrel{}{𝐟}}{C_n}\frac{\stackrel{}{𝚽}}{t}\frac{\stackrel{}{𝐠}}{C_n}\stackrel{}{𝚽}.$$
(65)
The Lagrange brackets depend exclusively on the functional form of $`x,y,z=f_{1,2,3}(C_{1,\mathrm{},6},t)`$ and $`g_{1,2,3}f_{1,2,3}/t`$, and are independent from the gauge and from the time evolution of $`C_i`$. Hence the gauge-invariant generalisation of the Lagrange system will emerge:

$$\frac{da}{dt}=\frac{2}{na}\left[\frac{R}{M_o}\stackrel{}{𝚽}\frac{\stackrel{}{𝐠}}{M_o}\frac{\stackrel{}{𝐟}}{M_o}\frac{d\stackrel{}{𝚽}}{dt}\right],$$
(66)
$`{\displaystyle \frac{de}{dt}}={\displaystyle \frac{1e^2}{na^2e}}\left[{\displaystyle \frac{R}{M_o}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{M_o}}{\displaystyle \frac{\stackrel{}{𝐟}}{M_o}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right]`$
$`{\displaystyle \frac{(1e^2)^{1/2}}{na^2e}}\left[{\displaystyle \frac{R}{a}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{a}}{\displaystyle \frac{\stackrel{}{𝐟}}{a}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right],`$ (67)
$`{\displaystyle \frac{d\omega }{dt}}={\displaystyle \frac{\mathrm{cos}i}{na^2(1e^2)^{1/2}\mathrm{sin}i}}\left[{\displaystyle \frac{R}{i}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{i}}{\displaystyle \frac{\stackrel{}{𝐟}}{i}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right]+`$
$`{\displaystyle \frac{(1e^2)^{1/2}}{na^2e}}\left[{\displaystyle \frac{R}{e}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{e}}{\displaystyle \frac{\stackrel{}{𝐟}}{e}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right],`$ (68)
$`{\displaystyle \frac{di}{dt}}={\displaystyle \frac{\mathrm{cos}i}{na^2(1e^2)^{1/2}\mathrm{sin}i}}\left[{\displaystyle \frac{R}{\omega }}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{\omega }}{\displaystyle \frac{\stackrel{}{𝐟}}{\omega }}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right]`$
$`{\displaystyle \frac{1}{na^2(1e^2)^{1/2}\mathrm{sin}i}}\left[{\displaystyle \frac{R}{\mathrm{\Omega }}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{\mathrm{\Omega }}}{\displaystyle \frac{\stackrel{}{𝐟}}{\mathrm{\Omega }}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right],`$ (69)
$$\frac{d\mathrm{\Omega }}{dt}=\frac{1}{na^2(1e^2)^{1/2}\mathrm{sin}i}\left[\frac{R}{i}\stackrel{}{𝚽}\frac{\stackrel{}{𝐠}}{i}\frac{\stackrel{}{𝐟}}{i}\frac{d\stackrel{}{𝚽}}{dt}\right],$$
(70)
$`{\displaystyle \frac{dM_o}{dt}}={\displaystyle \frac{1e^2}{na^2e}}\left[{\displaystyle \frac{R}{e}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{e}}{\displaystyle \frac{\stackrel{}{𝐟}}{e}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right]{\displaystyle \frac{2}{na}}\left[{\displaystyle \frac{R}{a}}\stackrel{}{𝚽}{\displaystyle \frac{\stackrel{}{𝐠}}{a}}{\displaystyle \frac{\stackrel{}{𝐟}}{a}}{\displaystyle \frac{d\stackrel{}{𝚽}}{dt}}\right].`$ (71)
These gauge-invariant equations reveal the potential possibility of simplification of orbit integration. One can deliberately choose gauges different from (46). In principle, it is possible to pick up the gauge so as to nullify the right-hand sides in three out of six equations (66 \- 71). This possibility is worth probing (we know from electrodynamics that a clever choice of gauge considerably simplifies solution of the equations of motion). Slabinski (2003) recently offered an example of gauge that nullifies the rhs of equation (69).

Another tempting possibility may be to pick up the gauge so that the $`\stackrel{\mathbf{}}{𝚽}`$-terms in (66 \- 71) fully compensate the short-period terms of the disturbing functions, leaving only the secular and resonant ones. However, the afore discussed Newman’s example strongly speaks against this possibility, at least in the case of bound orbits. The case of flybys seems to be more favourable, because in that case we do not have two different time scales. Hence in that case a choice of some nonvanishing $`\stackrel{\mathbf{}}{𝚽}`$ may, potentially, lead to an even larger simplification of calculations. We shall address this matter in a separate work.

Finally, we would point out that, in case the gauge $`\stackrel{\mathbf{}}{𝚽}`$ depends not only upon time but also upon the ”constants” $`C_i`$, the right-hand sides of (66) - (71) implicitly contain time derivatives of $`C_i`$. Then, in order to continue with analytic developments, it is necessary to transfer such terms to the left-hand side, as we did in (65). This alteration will be dwelled upon in (Efroimsky & Goldreich 2003a,b).

## 5 Delaunay’s variables

As well known (Brouwer & Clemence 1961, p. 290), it is possible to choose as the parameters $`C_i`$ not the six Keplerian elements $`C_i=(e,a,M_o,\omega ,\mathrm{\Omega },i)`$ but the set $`D_i=(L,G,H,M_o,\omega ,\mathrm{\Omega })`$, new variables $`L`$, $`G`$, and $`H`$ being defined as

$`L\mu ^{1/2}a^{1/2},G\mu ^{1/2}a^{1/2}\left(1e^2\right)^{1/2},H\mu ^{1/2}a^{1/2}\left(1e^2\right)^{1/2}\mathrm{cos}i,`$ (72)
where $`\mu G(m_{sun}+m_{planet})`$.

The advantage of these, Delaunay, variables lies in the diagonality of their Lagrange-bracket matrix. Inversion thereof yields, similarly to (66) - (71), the so-called Delaunay system:

$`{\displaystyle \frac{dL}{dt}}={\displaystyle \frac{R}{M_o}},{\displaystyle \frac{dM_o}{dt}}={\displaystyle \frac{R}{L}},`$
$`{\displaystyle \frac{dG}{dt}}={\displaystyle \frac{R}{\omega }},{\displaystyle \frac{d\omega }{dt}}={\displaystyle \frac{R}{G}},`$ (73)
$`{\displaystyle \frac{dH}{dt}}={\displaystyle \frac{R}{\mathrm{\Omega }}},{\displaystyle \frac{d\mathrm{\Omega }}{dt}}={\displaystyle \frac{R}{H}},`$
provided these parameters obey the Lagrange-type gauge condition analogous to (46):

$$\underset{i}{}\frac{\stackrel{~}{\stackrel{}{𝐟}}}{D_i}\frac{dD_i}{dt}=\mathrm{\hspace{0.33em}0}.$$
(74)
where, similarly to (36 \- 37),

$$\stackrel{\mathbf{}}{𝒓}=\stackrel{~}{\stackrel{}{𝐟}}(D_{1,\mathrm{},6},t),\dot{\stackrel{\mathbf{}}{𝒓}}=\stackrel{~}{\stackrel{}{𝐠}}(D_{1,\mathrm{},6},t),\frac{\stackrel{~}{\stackrel{}{𝐟}}}{t}=\stackrel{~}{\stackrel{}{𝐠}},$$
(75)
and the tilde symbol emphasises that the functional dependencies of the position and velocity upon $`D_i`$ differ from their dependencies upon $`C_i`$. We, thus, must keep in mind that the system of equations for the Delaunay elements only pretends to exist in a 6-dimensional phase space. In reality, it lives on a 9-dimensional submanifold (74) of a 12-dimensional manifold spanned by the Delaunay elements and their time derivatives. In the case of analytical calculations this, of course, makes no difference. Not in the case of numerical computation, though.

If instead of the gauge condition (74) we impose some alternative gauge

$$\underset{i}{}\frac{\stackrel{~}{\stackrel{}{𝐟}}}{D_i}\frac{dD_i}{dt}=\stackrel{~}{\stackrel{}{𝚽}}(t,D_{1,\mathrm{},6}),$$
(76)
the generalised Delaunay-type equations will read:131313 Similarly to our comment in the end of Section 4 we would mention that if the gauge $`\stackrel{~}{\stackrel{\mathbf{}}{𝚽}}`$ depends not only on time but also on the parameters $`C_i`$ then the right-hand sides of (77) contain time derivatives of $`C_i`$. A further analytic treatment of (77) will then demand that we transfer those terms to the left-hand side, as in (65). This will be discussed in (Efroimsky & Goldreich 2003a,b).

$`{\displaystyle \frac{dL}{dt}}={\displaystyle \frac{R}{M_o}}\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{M_o}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{M_o}}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},{\displaystyle \frac{dM_o}{dt}}={\displaystyle \frac{R}{L}}+\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{L}}+{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{L}}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},`$
$`{\displaystyle \frac{dG}{dt}}={\displaystyle \frac{R}{\omega }}\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{\omega }}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{\omega }}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},{\displaystyle \frac{d\omega }{dt}}={\displaystyle \frac{R}{G}}+\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{G}}+{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{G}}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},`$ (77)
$`{\displaystyle \frac{dH}{dt}}={\displaystyle \frac{R}{\mathrm{\Omega }}}\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{\mathrm{\Omega }}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{\mathrm{\Omega }}}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},{\displaystyle \frac{d\mathrm{\Omega }}{dt}}={\displaystyle \frac{R}{H}}+\stackrel{~}{\stackrel{}{𝚽}}{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐠}}}{H}}+{\displaystyle \frac{\stackrel{~}{\stackrel{}{𝐟}}}{H}}{\displaystyle \frac{d\stackrel{~}{\stackrel{}{𝚽}}}{dt}},`$
and the $`\stackrel{}{𝚽}`$ terms should not be ignored, because they account for the trajectory’s deviation from the submanifold (74) of the ambient 12-dimensional space $`(D_{1,\mathrm{},6},\dot{D}_{1,\mathrm{},6})`$.

The meaning of $`\stackrel{~}{\stackrel{}{𝐟}}`$ and $`\stackrel{~}{\stackrel{}{𝐠}}`$ in the above formulae is different than that of $`\stackrel{}{𝐟}`$ and $`\stackrel{}{𝐠}`$ in Sections 3 and 4. In those sections $`\stackrel{}{𝐟}`$ and $`\stackrel{}{𝐠}`$ denoted the functional dependencies (36) of $`x,y,z`$ and $`\dot{x},\dot{y},\dot{z}`$ upon parameters $`C_i=(e,a,M_o,\omega ,\mathrm{\Omega },i)`$. Here, $`\stackrel{~}{\stackrel{}{𝐟}}`$ and $`\stackrel{~}{\stackrel{}{𝐠}}`$ stand for the dependencies of $`x,y,z`$ and $`\dot{x},\dot{y},\dot{z}`$ upon the different set $`D_i=(L,G,H,M_o,\omega ,\mathrm{\Omega })`$. Despite the different functional forms, the values of $`\stackrel{~}{\stackrel{}{𝐟}}`$ and $`\stackrel{~}{\stackrel{}{𝐠}}`$ coincide with those of $`\stackrel{}{𝐟}`$ and $`\stackrel{}{𝐠}`$:

$$\stackrel{~}{\stackrel{}{𝐟}}\left(D_{1,\mathrm{},6}\right)=\stackrel{\mathbf{}}{𝒓}=\stackrel{}{𝐟}\left(C_{1,\mathrm{},6}\right)\text{and}\stackrel{~}{\stackrel{}{𝐠}}\left(D_{1,\mathrm{},6}\right)=\dot{\stackrel{\mathbf{}}{𝒓}}=\stackrel{}{𝐠}\left(C_{1,\mathrm{},6}\right).$$
(78)
Similarly, $`\stackrel{~}{\stackrel{}{𝚽}}\left(D_{1,\mathrm{},6}\right)`$ and $`\stackrel{}{𝚽}\left(C_{1,\mathrm{},6}\right)`$ are different functional dependencies. It is, though, easy to show (using the differentiation chain rule) that their values do coincide:

$$\stackrel{~}{\stackrel{}{𝚽}}\left(D_{1,\mathrm{},6}\right)=\stackrel{}{𝚽}\left(C_{1,\mathrm{},6}\right)$$
(79)
which is analogous to the covariance of Lorentz gauge in electrodynamics. This means that, for example, analytical calculations carried out by means of the Lagrange system (59 \- 64) are indeed equivalent to those performed by means of the Delaunay system (73), because imposure of the Lagrange gauge $`\stackrel{}{𝚽}=\mathbf{\hspace{0.33em}0}`$ is equivalent to imposure of $`\stackrel{~}{\stackrel{}{𝚽}}=\mathrm{\hspace{0.33em}0}`$.

Can one make a similar statement about numerical integrations? This question is nontrivial. In order to tackle it, we should recall that in the computer calculations the Lagrange condition $`\stackrel{}{𝚽}=\mathbf{\hspace{0.33em}0}`$ cannot be imposed exactly, for the numerical error will generate $`some`$ nonzero $`\stackrel{}{𝚽}`$. In other words, the orbit will never be perfectly constrained to the submanifold $`\stackrel{}{𝚽}=\mathbf{\hspace{0.33em}0}`$. Thereby, some nonzero $`\stackrel{}{𝚽}`$ will, effectively, emerge in (66 \- 71). Similarly, a small nonzero $`\stackrel{~}{\stackrel{}{𝚽}}`$ will, effectively, appear in (74), and the Delaunay system will no longer be canonical. In other words, we get not just an error in integration of the canonical system, but we get an error that drives the system of equation away from canonicity. This effect is not new: it is well known that not every numerical method preserves the Hamiltonian structure. Therefore, the unavoidable emergence of a nonzero numerical-error-caused $`\stackrel{~}{\stackrel{}{𝚽}}`$ in the Delaunay system may, potentially, be a hazard. This topic needs further investigation.

The gauge-invariant Delaunay-type system (77) is no longer Hamiltonian, unless the gauge is that of Lagrange, $`\stackrel{\mathbf{}}{𝚽}=\mathrm{\hspace{0.17em}0}`$. It is possible to demonstrate that by introducing a velocity-dependent disturbing force (which can be done, for example, by defining the orbital elements in a non-inertial reference frame) we further alter the Lagrange-type and Delaunay-type gauge-invariant systems. The gauge-invariant Delaunay-type system will remain non-canonical, except in one special gauge which is called in (Efroimsky & Goldreich 2003a) ”the generalised Lagrange gauge” and which simplifies to the regular Lagrange gauge when the velocity dependence is removed. Applications of this construction will be explained in Efroimsky & Goldreich (2003b).

To draw to a close, we would mention that in most books the Lagrange and Delaunay systems of equations are derived not by a straightforward application of the VOP method but through the medium of Hamilton-Jacobi technique (Plummer 1918, Smart 1953, Pollard 1966, Kovalevsky 1967, Stiefel and Scheifele 1971). At first glance it may seem that such a derivation does not exploit any sort of supplementary constraint and is, therefore, spared of the gauge invariance. In our upcoming publication Efroimsky & Goldreich (2003a) we shall explain that the Lagrange constraint is implicitly imposed within the Hamilton-Jacobi treatment of this problem. Another method of derivation, different from both the VOP and the Hamilton-Jacobi methods, was offered in the book by Kaula (1968) and, once again, it may seem that his development is devoid of whatever extra conditions. In the next section we explain how Kaula’s method, too, tacitly exploits the Lagrange constraint.

## 6 Brouwer and Kaula

As we saw above, the gauge-invariant generalisation of the Delaunay system is no longer canonical. In celestial mechanics, equations of type

$`\dot{r}={\displaystyle \frac{H}{p}}+X(r,p)`$
(80)
$`\dot{p}={\displaystyle \frac{H}{r}}+P(r,p)`$
were introduced by Dirk Brouwer who used them to include the atmospheric-drag forces into the canonical picture. The quantities $`X`$ and $`P`$ are called ”canonical forces”. The appropriate formalism is comprehensively explained in Stiefel & Scheifele (1971). We see now that such ”forces” can emerge not only due to real physical interactions but also as purely mathematical artifacts called into being by a nontrivial gauge choice. This in no sense means that such terms may be omitted. Their relevance to the stability of numerical methods will be discussed in Murison & Efroimsky (2003).

Above we saw that for a disturbance of type (57) there exists one special gauge (the Lagrange one) that makes the equations for Delaunay’s variables canonical. Still, in the case of an arbitrary gauge the symplectic structure is destroyed.141414 It will be proven in our next paper (Efroimsky & Goldreich 2003a) that in the more general case of a disturbing force dependent on positions and velocities the situation is similar: the equations for the Delaunay elements will become symplectic in a so-called generalised Lagrange gauge. To better understand the anatomy of this phenomenon, let us explore an alternative derivation of Delaunay’s equations, one based on a direct change of variables. This method is akin to the VOP technique, but has a better illustrative power. To simplify things, we shall assume that the disturbing force depends on the positions solely. Following Kaula (1968), we shall begin with the canonical equations for the position and velocity in some inertial Cartesian frame:

$`{\displaystyle \frac{d\stackrel{\mathbf{}}{𝒓}}{dt}}=\dot{\stackrel{\mathbf{}}{𝒓}},`$
(81)
$`{\displaystyle \frac{d\dot{\stackrel{\mathbf{}}{𝒓}}}{dt}}={\displaystyle \frac{U}{\stackrel{\mathbf{}}{𝒓}}}.`$
A solution to these will depend upon seven quantities $`S_0t,S_1,\mathrm{},S_6`$. One of these will be time, another six will be the parameters which can be devoid of or imparted with a time dependence of their own. This paves our way to the new, modified, system of equations:

$`{\displaystyle \frac{x_i}{S_k}}{\displaystyle \frac{dS_k}{dt}}=\dot{r_i},`$
(82)
$`{\displaystyle \frac{\dot{r}_i}{S_k}}{\displaystyle \frac{dS_k}{dt}}={\displaystyle \frac{U}{r_i}},`$
summation over repeating indices, from $`\mathrm{\hspace{0.33em}0}`$ through $`\mathrm{\hspace{0.33em}6}`$, being implied. After we multiply the upper and the lower equations by $`\dot{r}_i/S_l`$ and $`r_i/S_l`$, correspondingly, we should sum up both lines and arrive to:

$`[S_l,S_k]{\displaystyle \frac{dS_k}{dt}}={\displaystyle \frac{}{S_l}}\left(U{\displaystyle \frac{1}{2}}\dot{r}_i\dot{r}_i\right),[S_l,S_k]{\displaystyle \frac{r_i}{S_l}}{\displaystyle \frac{\dot{r}_i}{S_k}}{\displaystyle \frac{\dot{r}_i}{S_l}}{\displaystyle \frac{r_i}{S_k}}.`$ (83)
What remains now is to calculate the Lagrange-bracket matrix $`[S_l,S_k]`$ and to invert the result. This will yield the desired expressions for $`dS_k/dt`$. The expression for $`k=\mathrm{\hspace{0.17em}0}`$ will be merely an identity $`\mathrm{\hspace{0.33em}1}=\mathrm{\hspace{0.17em}1}`$, because $`S_0t`$. The other six will be the planetary equations. In case the parameters $`S_1,\mathrm{},S_6`$ make a Delaunay set, these equations will be expected to have a Hamiltonian form. Have they really? To check this, one has to calculate the appropriate Lagrange brackets. Or to accurately examine how these are calculated in the literature. The most direct way of computing the brackets is to differentiate formulae (40) with respect to the orbital elements. This will give the well-known time-independent expressions for the Lagrange brackets. These often-quoted standard expressions will, though, be valid in the two-body case solely. The latter circumstance is missing in (Kaula 1968) where the author implies that all will work equally well in the case of N bodies. In reality, two types of additional inputs will show themselves in the above formula for $`\dot{\stackrel{}{𝐫}}`$ in the N-body problem: First of all, extra items will appear in the expression for $`\dot{\stackrel{}{𝐪}}`$:

$`\dot{\stackrel{}{𝐪}}={\displaystyle \frac{\stackrel{}{𝐪}}{t}}+{\displaystyle \frac{\stackrel{}{𝐪}}{a}}{\displaystyle \frac{da}{dt}}+{\displaystyle \frac{\stackrel{}{𝐪}}{e}}{\displaystyle \frac{de}{dt}}+{\displaystyle \frac{\stackrel{}{𝐪}}{M_o}}{\displaystyle \frac{dM_o}{dt}}.`$ (84)
Beside this, the contribution $`\dot{𝐑}\stackrel{}{𝐪}`$ will now enter the expression for $`\dot{\stackrel{}{𝐱}}`$. Here

$`\dot{𝐑}={\displaystyle \frac{𝐑}{\mathrm{\Omega }}}{\displaystyle \frac{d\mathrm{\Omega }}{dt}}+{\displaystyle \frac{𝐑}{\omega }}{\displaystyle \frac{d\omega }{dt}}+{\displaystyle \frac{𝐑}{i}}{\displaystyle \frac{di}{dt}}.`$ (85)
Altogether, these will result in the following expression for the velocity:

$`\dot{\stackrel{\mathbf{}}{𝒓}}=𝐑(\mathrm{\Omega },i,\omega )\dot{\stackrel{}{𝐪}}(a,e,M)+\dot{𝐑}\stackrel{}{𝐪}`$
$`=𝐑{\displaystyle \frac{\stackrel{}{𝐪}}{t}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{a}}{\displaystyle \frac{da}{dt}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{e}}{\displaystyle \frac{de}{dt}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{M_o}}{\displaystyle \frac{dM_o}{dt}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{\mathrm{\Omega }}}{\displaystyle \frac{d\mathrm{\Omega }}{dt}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{\omega }}{\displaystyle \frac{d\omega }{dt}}+{\displaystyle \frac{(𝐑\stackrel{}{𝐪})}{i}}{\displaystyle \frac{di}{dt}}`$
$`=𝐑{\displaystyle \frac{\stackrel{}{𝐪}}{t}}+\stackrel{}{𝚽}.`$ (86)
In Kaula (1968) the latter term is ignored, so that the author implicitly imposes the Lagrange gauge. In this gauge, differentiation of the above expression with respect to the elements will indeed entail the standard time-independent expressions for the Lagrange brackets. In a non-Lagrange gauge, though, the situation will be different, and we shall not end up with a canonical Delaunay system. Instead, we shall arrive to the non-canonical gauge-invariant Delaunay-type system (77).

Very similarly, the Lagrange constraint enters all the methods by which the Gauss system of equations is derived in the literature, and the imposure of this constraint is camouflaged in a manner similar to what we saw above in regard to the Delaunay equations. We shall not engage in a comprehensive discussion of this issue, but shall rather provide a typical example. In section 11.5 of (Danby 1988) an auxiliary vector $`\widehat{𝐡}`$ is introduced as a unit vector aimed along the instantaneous orbital momentum of the body, relative to the gravitating centre. Then $`\widehat{𝐡}`$ gets resolved along the inertial axes $`(x,,y,z)`$, where $`\widehat{𝐱}`$ points toward the vernal equinox and $`\widehat{𝐳}`$ toward the north of the ecliptic pole:

$`\widehat{𝐡}=\widehat{𝐱}\mathrm{sin}\mathrm{\Omega }\mathrm{sin}i\widehat{𝐲}\mathrm{cos}\mathrm{\Omega }\mathrm{sin}i+\widehat{𝐳}\mathrm{cos}i.`$
This expression is certainly correct in the two-body case. It remains valid also in the N-body problem, only if the orbital elements $`\mathrm{\Omega }`$ and $`i`$ are osculating, i.e., only if the instantaneous inertial velocity is tangential to the ellipse parametrised by the Keplerian set that includes these $`\mathrm{\Omega }`$ and $`i`$. Thus, the Lagrange constraint is implied.

## 7 Is it worth it?

At this point one may ask if it is at all worth taking the nonzero $`\stackrel{~}{\stackrel{}{𝚽}}`$ into account in the Delaunay equations. After all, one can simply restrict himself to the 6-dimensional phase space defined by $`D_i`$, and postulate that the six unwanted extra dimensions $`\dot{D}_i`$ do not exist (i.e., postulate that $`\stackrel{~}{\stackrel{}{𝚽}}=\mathrm{\hspace{0.33em}0}`$). This, of course, can be done, but only at some cost: a certain type of accumulating numerical error will be ignored (not eliminated), and it will keep accumulating. As explained in the end of the previous section, the overall integration error of a Hamiltonian system consists of an error that leaves the system canonical (like, for example, an error in calculation of $`R`$ in (77)) and an error that drives the system away from its canonicity (like the error reflected in the accumulated nonzero value of $`\stackrel{}{𝚽})`$. The $`\stackrel{}{𝚽}`$ terms in (77) play the role of correctors: they fully compensate for the errors of the second type (i.e., for what in numerical electrodynamics is called gauge shift).

Similarly, in the case of Lagrange system, one may enquire if it is worth introducing the 12-dimensional space spanned by the orbital elements $`C_i`$ and their time derivatives $`H_i\dot{C}_i`$. Why not simply consider a trajectory in the 6-dimensional space of $`C_i`$ and assume that the Lagrange gauge is fixed exactly? Indeed, if we are solving the problem $`\ddot{\stackrel{}{𝐫}}=\stackrel{}{𝐟}(\stackrel{}{𝐫})`$, is it worth introducing an extra entity $`\stackrel{}{𝐯}\dot{\stackrel{}{𝐫}}`$ and considering the orbit integration in the space of a larger dimension, spanned by the components of both $`\stackrel{}{𝐫}`$ and $`\stackrel{}{𝐯}`$? Will this new entity $`\stackrel{}{𝐯}`$ add anything? The answer to this question will be affirmative if we take into account the fact that a trajectory is not merely a locus of points visited by the body: the notion of trajectory also comprises the rate at which the body was travelling. Appropriately, the accumulated numerical error will consist of two parts: distortions of the orbit shape, and distortions in the time-dependence of the speed at which the orbit was described. This explains why the events are taking place not just in the space of orbital elements but in the larger space of the elements and their time derivatives. This issue is best of all explained by Hagihara (1970). In subsection 1.6 of the first volume of his treatise he contrasts the cases of exact and apparent equivalences of dynamical systems. (The latter case is that of the orbital curves being geometrically, not dynamically, identical.)

Still, there is more to it, because a convenient choice of gauge may simplify the solution of the equations of motion (Slabinski 2003; Efroimsky & Goldreich 2003b).

## 8 Conclusions

We have demonstrated a previously unrecognised aspect of the Lagrange and Delaunay systems of planetary equations. Due to the Lagrange gauge condition (46), the motion is, in both cases, constrained to a 9-dimensional submanifold of the ambient 12-dimensional space spanned by the orbital elements and their time derivatives. Similarly to the field theory, the choice of gauge is vastly ambiguous and reveals a hidden symmetry (and an appropriate symmetry group) inherent in the description of the N-body problem in terms of the orbital elements. Just as a choice of a particular gauge simplifies solution of the equations of motion in electrodynamics, an alternative (to that of Lagrange) choice of gauge can simplify orbit calculations. We have written down the generalised Lagrange-type (66 \- 71) and Delaunay-type (77) equations in their gauge-invariant, form (with no specific gauge imposed). We have pointed out that neither the Lagrange gauge (46) nor any other constraint is exactly preserved in the course of numerical computation. This may be a source of numerical instability.

Acknowledgements
Debts of gratitude are always pleasure to pay. I am grateful to William Newman for our numerous conversations on planetary dynamics. I also wish to acknowledge the critical questions and comments by Yurii Batrakov, Robert Burridge, Sergei Klioner, Mark Levi, Robert O’Malley, Scot Tremaine, and Jack Wisdom, which greatly helped me to improve the paper. My truly special thanks go to Peter Goldreich who kindly discussed with me some relevant concepts, and offered very important corrections and amendments to the manuscript. This research was supported by NASA grant W-19948.