© ESO 2020

1. Introduction

Perturbed Keplerian motion is a multiscale problem, in which the orbital elements evolve slowly when compared to the change with time of the ephemeris, whose fast evolution is determined by the rate of variation of the mean anomaly. Usual integration schemes that look for the separation of fast and slow frequencies of the motion may be superior to the simpler integration with the Cowell method, in which the Newtonian acceleration is directly integrated in rectangular coordinates, and formulations based on Hansen’s ideal frame concept (Hansen 1857) closely approach this decoupling (see Lara 2017, and references therein). Moreover, the use of symplectic integrators is widely adopted in the case of Hamiltonian perturbation problems (Laskar & Robutel 2001; Blanes et al. 2013).

On the other hand, the effective decoupling of short- and long-period effects is achieved with perturbation methods. After removing short-period effects in an averaging process, the long-term behavior of the (mean) orbital elements is efficiently integrated from the equations of variation of parameters. The parameters, which are constant in the pure Keplerian motion, can take different representations, such as classical Keplerian elements or nonsingular variables (see Hintz 2008 for a survey), but are always related to the two fundamental vectors of the orbital motion: the angular momentum vector, which determines the orientation of the orbital plane, and the eccentricity vector, which determines the shape and orientation of the reference ellipse in that plane. This last is a nondimensional version of the classical Laplace vector, which has dimensions of the gravitational parameter, or of the Runge-Lenz vector, which has dimensions of angular momentum. While vectorial formulations normally introduce redundancy by increasing the dimension of the differential system, they commonly admit a compact and symmetric formulation of the right side of the variation equations that renders faster evaluation (Allan 1962; Musen 1963; Roy & Moran 1973)¹, and may even disclose integrability (Deprit 1984; see also Mignard & Henon 1984; Richter & Keller 1995).

The vectorial formulation has been shown to be useful in the case of third-body perturbations (Breiter & Ratajczak 2005; Correia et al. 2011; Katz et al. 2011). In particular, it may be an efficient alternative to classical formulations when the orbits are highly elliptic, which is a common case in extrasolar planetary systems (Lee & Peale 2003; Migaszewski & Goździewski 2008), in artificial satellite theory (Lara et al. 2012, 2018), and in hierarchical N-body systems in general (Hamers et al. 2015; Will 2017). These types of orbits may need higher degrees in the expansion of the third-body disturbing function to provide a reasonable approximation of the dynamics even in the simplifications offered by the secular dynamics (Beaugé & Michtchenko 2003; Libert & Sansottera 2013; Andrade-Ines et al. 2016; Andrade-Ines & Robutel 2018; Sansottera & Libert 2019). This fact has motivated the recent appearance of general expressions of the expansion of the disturbing function in powers of the ratio of the disturbing and disturbed semimajor axes, including the secular terms (Laskar & Boué 2010; Mardling 2013; Palacián et al. 2017). These expressions rely on classical elements of the Keplerian motion and apply to any eccentricity and inclination.

We constrain ourselves to the restricted approximation, in which the lower mass body does not affect the motion of the primaries. This is a reasonable approximation in artificial satellite theory as well as in some problems of Solar System dynamics, but it does not apply to the dynamics of extrasolar planetary systems. For the restricted case, we make use of the apsidal frame formulation, which effectively displays the fast and slow components of the orbital motion, and we use perturbation theory to remove the short-period components of the third-body disturbing function. The third-body direction is assumed to be a known function of time, given by an ephemeris, but we avoid time-related issues in the averaging by constraining ourselves to the usual case in which the position of the third body can be taken as fixed during one orbital period. Instead of relying on the classical Hansen expansions or related eccentricity functions (Hansen 1855; Kaula 1962; Giacaglia 1974; Lane 1989; Celletti et al. 2017), we average the Legendre polynomial expansion of the third-body potential in closed form of the eccentricity after the usual reformulation in terms of the eccentric anomaly (Deprit 1983; Kelly 1989).

The decoupling of short-period effects is achieved after the standard Delaunay normalization (Deprit 1982), which we extend only to the first order in the Hamiltonian perturbation approach; that is, we do not consider the possible coupling between the different terms in which the third-body perturbation is expanded. However, while the normalization is properly carried out in Delaunay canonical variables, which are the action-angle variables of the Kepler problem, the canonical variables related to the Keplerian elements are collected in vectorial form throughout the whole normalization procedure. In this way, we obtain alternative, general expressions for the expansion of the third-body disturbing function in mean vectorial elements up to an arbitrary degree. For completeness, the generating function of the infinitesimal contact transformation from mean to osculating elements is also provided in vectorial elements up to an arbitrary degree of the expansion.

The Hamilton equations of the mean elements Hamiltonian are then computed to obtain analogous general expressions for the variation equations of the flow in Delaunay (mean) elements, which, as expected, are flawed by the appearance of the eccentricity and the sine of the inclination as divisors. Singularities are avoided by reformulating the flow in different sets of canonical and noncanonical variables. In particular, we provide vectorial, nonsingular expressions for the variation equations of the eccentricity and angular momentum vectors in mean elements. These expressions can be used to extend to an arbitrary degree existing lower-order truncations in the literature (Allan 1962; Katz et al. 2011).

2. Third-body disturbing potential

Let (O, i, j, k) be an orthonormal frame with origin O in the center of mass of a central attracting body and fixed directions defined by the unit vectors i, j, k. Let r and r_⋆ be the position vectors of a massless and a massive body, respectively, of modulus r = ∥r∥ and r_⋆ = ∥r_⋆∥. When r ≪ r_⋆, the third-body disturbing potential

$$\begin{array}{c}\hfill {\mathcal{V}}_{\star }=-\frac{{\mu }_{\star }}{r_{\star }}(\frac{r_{\star }}{\Vert \mathit{r}-{\mathit{r}}_{\star }\Vert }-\frac{\mathit{r}·{\mathit{r}}_{\star }}{r_{\star }^2}),\end{array}$$(1)

where μ_⋆ is the third-body gravitational parameter, is customarily replaced by the Legendre polynomial expansion

$$\begin{array}{c}\hfill {\mathcal{V}}_{\star }=-\frac{n_{\star }^2{a}_{\star }^3}{r_{\star }}{\chi }_{\star }{\displaystyle \sum _{i\ge 2}}\frac{r^i}{r_{\star }^i}{P}_i(cos{\psi }_{\star }),\end{array}$$(2)

in which μ_⋆ has been rewritten in terms of the mean motion n_⋆ and the semimajor axis a_⋆ of the third-body orbit relative to the central body, χ_⋆ is the mass ratio of the system, $cos{\psi }_{\star }=\widehat{\mathit{r}}·{\widehat{\mathit{r}}}_{\star }$, with $\widehat{\mathit{r}}=\mathit{r}/r$ and ${\widehat{\mathit{r}}}_{\star }={\mathit{r}}_{\star }/r_{\star }$, and the Legendre polynomials P_i are given by the usual expansion of Rodrigues’ formula,

$$\begin{array}{c}\hfill P_i(cos{\psi }_{\star })=\frac{1}{2^i}{\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}{(-1)}^l\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){cos}^{i-2l}{\psi }_{\star },\end{array}$$(3)

with ⌊i/2⌋ denoting the integer part of the division i/2. We recall that as usual, the term −μ_⋆/r_⋆ is neglected in Eq. (2) because it does not contribute to the disturbing acceleration of the massless body.

We define now the apsidal moving frame $(O,\widehat{\mathit{e}},\mathit{b},\mathit{n})$ with the same origin as before and the unit vectors n in the direction of the (instantaneous) angular momentum vector (per unit of mass),

$$\begin{array}{c}\hfill \mathit{G}=G\mathit{n}=\mathit{r}\times \mathit{v},\end{array}$$(4)

where G = ∥G∥ and v = dr/dt, with t denoting time. The unit vector $\widehat{\mathit{e}}$ has the direction of the (instantaneous) eccentricity vector,

$$\begin{array}{c}\hfill \mathit{e}=e\widehat{\mathit{e}}=\frac{\mathit{v}\times \mathit{G}}{\mu }-\frac{\mathit{r}}{r},\end{array}$$(5)

where e = ∥e∥ is the (instantaneous) eccentricity of the orbit of the massless body, and μ is the gravitational parameter of the central body. Finally, the unit vector

$$\begin{array}{c}\hfill \mathit{b}=\mathit{n}\times \widehat{\mathit{e}}\end{array}$$(6)

defines the binormal direction, in this way completing a direct orthonormal frame.

For convenience, we define the unit vector ℓ in the direction of the ascending node, given by

$$\begin{array}{c}\hfill \mathit{k}\times \mathit{n}=\mathit{\ell }sinI,\end{array}$$(7)

where

$$\begin{array}{c}\hfill I=arccos(\mathit{k}·\mathit{n}),\end{array}$$(8)

is the inclination angle of the (instantaneous) orbital plane with respect to the (i, j) plane. In addition,

$$\begin{array}{c}\hfill \mathit{\ell }·\widehat{\mathit{e}}=cos\omega ,\mathit{\ell }·\mathit{b}=-sin\omega ,\end{array}$$(9)

where ω is the argument of the pericenter. Then, it is simple to check that

$$\begin{array}{cc}\hfill \widehat{\mathit{e}}·\mathit{k}& =sinIsin\omega =-\mathit{\ell }·\mathit{b}sinI,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill \mathit{b}·\mathit{k}& =sinIcos\omega =\mathit{\ell }·\widehat{\mathit{e}}sinI.\hfill \end{array}$$(11)

Finally, by denoting with h the longitude of the node, the components of the apsidal frame in the inertial frame are given by the rotations

$$\begin{array}{c}\hfill (\widehat{\mathit{e}},\mathit{b},\mathit{n})=R_3(-h)R_1(-I)R_3(-\omega ),\end{array}$$(12)

where

$$\begin{array}{cc}\hfill R_1(\alpha )& =\left(\begin{array}{ccc}1& 0& 0\\ 0& cos\alpha & sin\alpha \\ 0& -sin\alpha & cos\alpha \end{array}\right),\hfill \\ \hfill R_3(\alpha )& =\left(\begin{array}{ccc}cos\alpha & sin\alpha & 0\\ -sin\alpha & cos\alpha & 0\\ 0& 0& 1\end{array}\right)\hfill \end{array}$$

are standard rotation matrices.

When the direction of the massless body is given by its components in the apsidal frame,

$$\begin{array}{c}\hfill \widehat{\mathit{r}}=\widehat{\mathit{e}}cosf+\mathit{b}sinf,\end{array}$$(13)

with f being the true anomaly, we obtain

$$\begin{array}{c}\hfill cos{\psi }_{\star }=(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })cosf+(\mathit{b}·{\widehat{\mathit{r}}}_{\star })sinf,\end{array}$$(14)

which immediately discloses the short-period terms affecting Eq. (3).

However, the short-period terms of the third-body disturbing potential in Eq. (2) are not limited to those contributed by the Legendre polynomials, and they are better handled when written in terms of the eccentric anomaly u, in contrast to the true anomaly. This is done using the geometric relations

$$\begin{array}{c}\hfill sinf=\frac{a}{r}\eta sinu,cosf=\frac{a}{r}(cosu-e),\end{array}$$(15)

where the eccentricity function η = (1 − e²)^1/2 is used for convenience. We recall also that

$$\begin{array}{c}\hfill r=a(1-ecosu),\end{array}$$(16)

and after some rearrangement, we obtain

$$\begin{array}{c}\hfill {\mathcal{V}}_{\star }=-\frac{\mu }{2a}{\chi }_{\star }\frac{n_{\star }^2}{n^2}\frac{a_{\star }^3}{r_{\star }^3}\frac{a}{r}{\displaystyle \sum _{i\ge 2}}\frac{1}{2^{i-1}}\frac{a^{i-2}}{r_{\star }^{i-2}}{V}_i,\end{array}$$(17)

in which the nondimensional functions V_i take the form

$$\begin{array}{cc}\hfill V_i=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}{(-1)}^l\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right)\frac{r^{2l+1}}{a^{2l+1}}[(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })(cosu-e)\hfill \\ & +{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })\eta sinu]}^{i-2l}.\hfill \end{array}$$(18)

The reasons for the particular arrangement of Eq. (17), in which we excluded the coefficient a/r from the summation, is that the removal of short-period terms from 𝒱_⋆, which is carried out by Eq. (24), is easily achieved in closed form of the eccentricity by taking advantage of the differential relation between the mean and true anomalies.

Next, we substitute Eq. (16) into Eq. (18) and use the binomial expansion to obtain

$$\begin{array}{c}\hfill \frac{r^{2l+1}}{a^{2l+1}}={(1-ecosu)}^{2l+1}={\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){(-1)}^j{e}^j{cos}^{j}u.\end{array}$$(19)

Analogously, the binomial expansion is applied to the last term in the right side of Eq. (18), to yield

$$\begin{array}{cc}& {[(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })(cosu-e)+(\mathit{b}·{\widehat{\mathit{r}}}_{\star })\eta sinu]}^{i-2l}\hfill \\ \hfill & ={\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(cosu-e)}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{i-2l-k}{\eta }^{i-2l-k}{sin}^{i-2l-k}u,\hfill \end{array}$$(20)

where

$$\begin{array}{c}\hfill {(cosu-e)}^k={\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right){(-1)}^m{e}^m{cos}^{k-m}u.\end{array}$$(21)

Plugging Eqs. (19)–(21) into Eq. (18), we finally obtain

$$\begin{array}{cc}\hfill V_i=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right)e^{j+m}\hfill \\ & \times {(-1)}^{j+l+m}{\eta }^{i-2l-k}{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{i-2l-k}\hfill \\ & \times {cos}^{j+k-m}u{sin}^{i-2l-k}u.\hfill \end{array}$$(22)

3. Short-period averaging

Because third-body perturbations derive from the potential in Eq. (2), we can take advantage of the Hamiltonian formalism. Thus

$$\begin{array}{c}\hfill \mathcal{H}=-\frac{\mu }{2a}+{\mathcal{V}}_{\star },\end{array}$$(23)

where ℋ must be expressed in some set of canonical variables. In particular, we use Delaunay variables (ℓ, g, h, L, G, H), where the coordinates are the mean anomaly ℓ, the argument of the periapsis g, and the longitude of the node h, and their conjugate momenta are the Delaunay action $L=\sqrt{\mu a}$, the modulus of the angular momentum vector G, and its projection along the k direction H, respectively. A modern derivation of this useful set of variables can be found, for instance, in Lara (2016).

Disregarding short-period effects, the orbit evolution is customarily studied in mean elements (ℓ′, g′,h′,L′,G′,H′), which aim to represent the average value of the true osculating variables. The transformation from mean to osculating variables is found using the tools of perturbation theory. To avoid time-dependency issues introduced by the third-body ephemeris ${\widehat{\mathit{r}}}_{\star }\equiv {\widehat{\mathit{r}}}_{\star }(t)$, we used the extended phase-space formulation. In particular, we relied on the Lie transforms method (Hori 1966; Deprit 1969), which is considered standard and is thoroughly described in the literature (see, e.g., Meyer & Hall 1992; Boccaletti & Pucacco 2002). The details of the transformation will be presented elsewhere, and we focus here on the averaging of the disturbing potential. However, because its computation is immediate once the disturbing potential has been averaged, we also provide the generating function from which the short-period corrections are derived.

The short-period elimination is carried out in closed form of the eccentricity. It is effectively achieved with the help of the differential relation between the eccentric and true anomalies dℓ=(1 − ecosu) du, which is obtained by differentiating Kepler’s equation. Thus, because of Eq. (16),

$$\begin{array}{c}\hfill {⟨{\mathcal{V}}_{\star }⟩}_{\ell }=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathcal{V}}_{\star }\mathrm{d}\ell =\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{\mathcal{V}}_{\star }\frac{r}{a}\mathrm{d}u,\end{array}$$(24)

from which, after replacing Eq. (17), we obtain

$$\begin{array}{c}\hfill {⟨{\mathcal{V}}_{\star }⟩}_{\ell }=-{\chi }_{\star }\frac{\mu }{2a}\frac{n_{\star }^2}{n^2}\frac{a_{\star }^3}{r_{\star }^3}{\displaystyle \sum _{i\ge 2}}\frac{1}{2^{i-1}}\frac{a^{i-2}}{r_{\star }^{i-2}}{⟨V_i⟩}_u,\end{array}$$(25)

where

$$\begin{array}{c}\hfill {⟨V_i⟩}_u=\frac{1}{2\pi }{{\displaystyle \int }}_0^{2\pi }{V}_i\mathrm{d}u.\end{array}$$(26)

In view of Eq. (22), the only terms that remain below the integral symbol of Eq. (26) are of the form sin^a u cos^bu with a and b integers. Terms of this type are easily integrated when they are expressed as Fourier series in u. Straightforward computations using standard relations between exponentials and circular functions (see, e.g., Kaula 1966), yield

$$\begin{array}{cc}\hfill {cos}^{j+k-m}u{sin}^{i-2l-k}u& =\frac{{(-\mathbf{i})}^{i-k-2l}}{2^{i+j-m-2l}}{\displaystyle \sum _{q=0}^{i-k-2l}}{\displaystyle \sum _{t=0}^{j+k-m}}\left(\begin{array}{c}i-k-2l\\ q\end{array}\right){(-1)}^q\hfill \\ \hfill & \times \left(\begin{array}{c}j+k-m\\ t\end{array}\right)\hfill \\ \hfill & \times [cos2(\stackrel{\sim }{q}-q)u+\mathbf{i}sin2(\stackrel{\sim }{q}-q)u],\hfill \end{array}$$(27)

where i notes the imaginary unit, and we abbreviated

$$\begin{array}{c}\hfill \stackrel{\sim }{q}=\frac{i+j-m}{2}-l-t.\end{array}$$(28)

Now, it becomes obvious that the only nonperiodic terms of Eq. (27)

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{j+k-m}^{i-2l-k}={\langle {cos}^{j+k-m}u{sin}^{i-2l-k}u\rangle }_u,\end{array}$$(29)

are those such that $q=\stackrel{\sim }{q}$, a condition that is only accomplished when i + j − m is even, which in turn implies that i + j + m is also even. In addition, because Eq. (27) must be free from imaginary terms, i − k, the exponent of the imaginary unit, must be even in the case of nonperiodic terms, yielding the numeric coefficient

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{2a}^{2b}=\frac{1}{2^{2a+2b}}{\displaystyle \sum _{t=0}^{2a}}\left(\begin{array}{c}2b\\ b+a-t\end{array}\right){(-1)}^{a-t}\left(\begin{array}{c}2a\\ t\end{array}\right),\end{array}$$(30)

where $a=\frac{1}{2}(j+k-m)$ and $b=\frac{1}{2}(i-k)-l$ are integer numbers.

On the other hand, i − k even allows us to replace η^{i − k − 2l} by the binomial expansion of ${(1-e^2)}^{\frac{1}{2}(i-k-2l)}$. Thus, replacing

$$\begin{array}{c}\hfill e^{m+j}{(1-e^2)}^{\frac{i-k}{2}-l}={\displaystyle \sum _{q=0}^{\frac{i-k}{2}-l}}\left(\begin{array}{c}\frac{i-k}{2}-l\\ q\end{array}\right){(-1)}^q{e}^{m+j+2q},\end{array}$$(31)

into, Eq. (26), we finally obtain

$$\begin{array}{cc}\hfill {⟨V_i⟩}_u=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right){\mathrm{\Xi }}_{j+k-m}^{i-2l-k}\hfill \\ & \times {\displaystyle \sum _{q=0}^{\frac{i-k}{2}-l}}\left(\begin{array}{c}\frac{i-k}{2}-l\\ q\end{array}\right){(-1)}^{i-l+q}{e}^{m+j+2q}{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{i-2l-k},\hfill \end{array}$$(32)

thus completing the calculation of Eq. (25).

The generating function of the infinitesimal contact transformation leading to the averaging is computed from the usual relation

$$\begin{array}{c}\hfill \mathcal{W}=\frac{1}{n}{\displaystyle \int }({\mathcal{V}}_{\star }-{⟨{\mathcal{V}}_{\star }⟩}_{\ell })\mathrm{d}\ell .\end{array}$$(33)

That is,

$$\begin{array}{cc}\hfill \mathcal{W}& =-{\chi }_{\star }L\frac{n_{\star }^2}{n^2}\frac{a_{\star }^3}{r_{\star }^3}{\displaystyle \sum _{i\ge 2}}\frac{1}{2^i}\frac{a^{i-2}}{r_{\star }^{i-2}}{\displaystyle \int }(\frac{a}{r}{V}_i-{⟨V_i⟩}_u)\mathrm{d}\ell \hfill \\ \hfill & =-{\chi }_{\star }L\frac{n_{\star }^2}{n^2}\frac{a_{\star }^3}{r_{\star }^3}{\displaystyle \sum _{i\ge 2}}\frac{1}{2^i}\frac{a^{i-2}}{r_{\star }^{i-2}}(-{⟨V_i⟩}_u\ell +{\displaystyle \int }{V}_i\mathrm{d}u),\hfill \end{array}$$(34)

which can be written in the form

$$\begin{array}{c}\hfill \mathcal{W}=-{\chi }_{\star }L\frac{n_{\star }^2}{n^2}\frac{a_{\star }^3}{r_{\star }^3}{\displaystyle \sum _{i\ge 2}}\frac{1}{2^i}\frac{a^{i-2}}{r_{\star }^{i-2}}[{⟨V_i⟩}_u(u-\ell )+{\displaystyle \int }(V_i-{⟨V_i⟩}_u)\mathrm{d}u],\end{array}$$(35)

where u − ℓ = e sin u, from the Kepler equation, and the integrand of the last term in the square brackets is only composed of periodic terms in the eccentric anomaly. Then, after replacing Eq. (27) into Eq. (22), we solve the indefinite integral by keeping all the trigonometric terms except for those that make i + j − m − 2l − 2q − 2t = 0. Namely,

$$\begin{array}{cc}\hfill {\displaystyle \int }(V_i-{⟨V_i⟩}_u)\mathrm{d}u=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right)\hfill \\ & \times {\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right)e^{m+j}{\eta }^{i-2l-k}{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{i-2l-k}\hfill \\ & \times {(-1)}^{j+l+m}\frac{{(-\mathbf{i})}^{i-2l-k}}{2^{i-2l+j-m}}{\displaystyle \sum _{\begin{array}{c}q=0\\ q\ne \stackrel{\sim }{q}\end{array}}^{i-2l-k}}{\displaystyle \sum _{t=0}^{j+k-m}}\left(\begin{array}{c}i-2l-k\\ q\end{array}\right)\hfill \\ & \times \left(\begin{array}{c}j+k-m\\ t\end{array}\right){(-1)}^q\hfill \\ & \times \frac{sin2(\stackrel{\sim }{q}-q)u-\mathbf{i}cos2(\stackrel{\sim }{q}-q)u}{2(\stackrel{\sim }{q}-q)},\hfill \end{array}$$(36)

where $sin2(\stackrel{\sim }{q}-q)u-\mathbf{i}cos2(\stackrel{\sim }{q}-q)u=-\mathbf{i}{(cosu+\mathbf{i}sinu)}^{2(\stackrel{\sim }{q}-q)}$ from which, using Eq. (15),

$$\begin{array}{c}\hfill sin2(\stackrel{\sim }{q}-q)u-\mathbf{i}cos2(\stackrel{\sim }{q}-q)u=-\mathbf{i}{(e+\frac{\mathit{r}·\widehat{\mathit{e}}}{a}+\mathbf{i}\frac{\mathit{r}·\mathit{b}}{a\eta })}^{2(\stackrel{\sim }{q}-q)}.\end{array}$$(37)

Finally, the Hamiltonian in mean elements, up to higher order effects, is obtained by replacing osculating by mean elements into both the Keplerian −μ/(2a) and the third-body averaged potential, Eq. (25), to yield

$$\begin{array}{c}\hfill \mathcal{K}=-\frac{\mu }{2a}+{⟨{\mathcal{V}}_{\star }⟩}_{\ell },\end{array}$$(38)

where now all the symbols appearing in 𝒦 are functions of the Delaunay prime elements.

4. Hamilton equations of the averaged flow

The flow in mean (prime) elements is obtained from the Hamilton equations. Because the mean anomaly has been removed by the perturbation approach up to the truncation order of the perturbation solution, the prime Delaunay action is a constant that decouples the reduced flow

$$\begin{array}{cc}& \frac{\mathrm{d}{g}^{\prime }}{\mathrm{d}t}=\frac{\partial \mathcal{K}}{\partial G^{\prime }}=\frac{\partial {⟨{\mathcal{V}}_{\star }⟩}_{\ell }}{\partial G^{\prime }},\hfill \end{array}$$(39)

$$\begin{array}{cc}& \frac{\mathrm{d}{G}^{\prime }}{\mathrm{d}t}=-\frac{\partial \mathcal{K}}{\partial g^{\prime }}=-\frac{\partial {⟨{\mathcal{V}}_{\star }⟩}_{\ell }}{\partial g^{\prime }},\hfill \end{array}$$(40)

$$\begin{array}{cc}& \frac{\mathrm{d}{h}^{\prime }}{\mathrm{d}t}=\frac{\partial \mathcal{K}}{\partial H^{\prime }}=\frac{\partial {⟨{\mathcal{V}}_{\star }⟩}_{\ell }}{\partial H^{\prime }},\hfill \end{array}$$(41)

$$\begin{array}{cc}& \frac{\mathrm{d}{H}^{\prime }}{\mathrm{d}t}=-\frac{\partial \mathcal{K}}{\partial h^{\prime }}=-\frac{\partial {⟨{\mathcal{V}}_{\star }⟩}_{\ell }}{\partial h^{\prime }},\hfill \end{array}$$(42)

from the integration of the variation of the prime mean anomaly,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\ell }^{\prime }}{\mathrm{d}t}=\frac{\partial \mathcal{K}}{\partial L^{\prime }}=n+\frac{\partial {⟨{\mathcal{V}}_{\star }⟩}_{\ell }}{\partial L^{\prime }},\end{array}$$(43)

which is computed by indefinite integration after solving the reduced system of Eqs. (39)–(42).

In fact, in view of Eqs. (25) and (32), computing the right sides of the variation Eqs. (39)–(42) only requires the corresponding partial derivatives of the function

$$\begin{array}{c}\hfill {\alpha }_{\sigma ,k,p}=e^p{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k},\sigma =i-2l,p=m+j+2q.\end{array}$$(44)

In Eq. (44), both direction vectors $\widehat{\mathit{e}}$ and b depend on the orbital inclination I, the argument of the periapsis g, and the longitude of the node h, as shown in Eq. (12).

The required partial derivatives of α_σ, k, p with respect to the Delaunay variables are computed using the chain rule:

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial G}& =\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial e}\frac{\partial e}{\partial G}+\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}\left(\frac{\partial \widehat{\mathit{e}}}{\partial I}\frac{\partial I}{\partial G}\right)·{\widehat{\mathit{r}}}_{\star }\hfill \\ \hfill & +\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\mathit{b}·{\widehat{\mathit{r}}}_{\star })}\left(\frac{\partial \mathit{b}}{\partial I}\frac{\partial I}{\partial G}\right)·{\widehat{\mathit{r}}}_{\star },\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial H}& =\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}\left(\frac{\partial \widehat{\mathit{e}}}{\partial I}\frac{\partial I}{\partial H}\right)·{\widehat{\mathit{r}}}_{\star }\hfill \\ \hfill & +\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\mathit{b}·{\widehat{\mathit{r}}}_{\star })}\left(\frac{\partial \mathit{b}}{\partial I}\frac{\partial I}{\partial H}\right)·{\widehat{\mathit{r}}}_{\star },\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial g}& =\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}\frac{\partial \widehat{\mathit{e}}}{\partial g}·{\widehat{\mathit{r}}}_{\star }+\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\mathit{b}·{\widehat{\mathit{r}}}_{\star })}\frac{\partial \mathit{b}}{\partial g}·{\widehat{\mathit{r}}}_{\star },\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial h}& =\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}\frac{\partial \widehat{\mathit{e}}}{\partial h}·{\widehat{\mathit{r}}}_{\star }+\frac{\partial {\alpha }_{\sigma ,k,p}}{\partial (\mathit{b}·{\widehat{\mathit{r}}}_{\star })}\frac{\partial \mathit{b}}{\partial h}·{\widehat{\mathit{r}}}_{\star }.\hfill \end{array}$$(48)

Because c ≡ cos I = H/G, s ≡ sin I = (1 − c²)^1/2, and e = (1 − G²/L²)^1/2, we easily obtain

$$\begin{array}{c}\hfill \frac{\partial e}{\partial G}=-\frac{1}{G}\frac{{\eta }^2}{e},\frac{\partial I}{\partial H}=-\frac{1}{G}\frac{1}{s},\frac{\partial I}{\partial G}=\frac{1}{G}\frac{c}{s}·\end{array}$$(49)

In addition, because the effect of a differential rotation about the axis n (resp. ℓ, k) is an infinitesimal increase of the angle g (resp. I, h), we find

$$\begin{array}{c}\hfill \frac{\partial \widehat{\mathit{e}}}{\partial g}=\mathit{n}\times \widehat{\mathit{e}}=\mathit{b},\frac{\partial \widehat{\mathit{e}}}{\partial h}=\mathit{k}\times \widehat{\mathit{e}},\frac{\partial \widehat{\mathit{e}}}{\partial I}=\mathit{\ell }\times \widehat{\mathit{e}}=\frac{1}{s}(\widehat{\mathit{e}}·\mathit{k})\mathit{n},\end{array}$$(50)

and

$$\begin{array}{c}\hfill \frac{\partial \mathit{b}}{\partial g}=\mathit{n}\times \mathit{b}=-\widehat{\mathit{e}},\frac{\partial \mathit{b}}{\partial h}=\mathit{k}\times \mathit{b},\frac{\partial \mathit{b}}{\partial I}=\mathit{\ell }\times \mathit{b}=\frac{1}{s}(\mathit{b}·\mathit{k})\mathit{n}.\end{array}$$(51)

Thus, straightforward computations yield

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial G}& =-\frac{1}{G}\frac{{\eta }^2}{e}\left[p{e}^{p-1}\right]{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k}\hfill \\ \hfill & +\frac{1}{G}\frac{c}{s}{e}^p\left[k{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^{k-1}\right]{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k}(\mathit{n}·{\widehat{\mathit{r}}}_{\star })\frac{\widehat{\mathit{e}}·\mathit{k}}{s}\hfill \\ \hfill & +\frac{1}{G}\frac{c}{s}{e}^p{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k\left[(\sigma -k){(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k-1}\right](\mathit{n}·{\widehat{\mathit{r}}}_{\star })\frac{\mathit{b}·\mathit{k}}{s},\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial H}& =-\frac{1}{G}\frac{1}{s}{e}^p\left[k{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^{k-1}\right]{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k}(\mathit{n}·{\widehat{\mathit{r}}}_{\star })\frac{\widehat{\mathit{e}}·\mathit{k}}{s}\hfill \\ \hfill & -\frac{1}{G}\frac{1}{s}{e}^p{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k\left[(\sigma -k){(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k-1}\right](\mathit{n}·{\widehat{\mathit{r}}}_{\star })\frac{\mathit{b}·\mathit{k}}{s},\hfill \end{array}$$(53)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial g}& =e^p\{\left[k{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^{k-1}\right]{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k+1}\hfill \\ \hfill & -{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^{k+1}\left[(\sigma -k){(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k-1}\right]\},\hfill \end{array}$$(54)

$$\begin{array}{cc}\hfill \frac{\partial {\alpha }_{\sigma ,k,p}}{\partial h}& =e^p\{\left[k{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^{k-1}\right]{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k}(\widehat{\mathit{e}}\times {\widehat{\mathit{r}}}_{\star })\hfill \\ \hfill & +{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k\left[(\sigma -k){(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^{\sigma -k-1}\right](\mathit{b}\times {\widehat{\mathit{r}}}_{\star })\}·\mathit{k},\hfill \end{array}$$(55)

where the dispensable square brackets in these equations are used to emphasize that corresponding enclosed terms never introduce denominators.

5. Long-term flow in mean vectorial elements

As expected from the singularities of Delaunay variables, Eqs. (52) and (53) are singular in the case of circular and equatorial orbits. However, these singularities are of virtual nature (Henrard 1974) and may be removed when the equations of the flow are represented in other sets of canonical or noncanonical variables. In particular, the variations experienced by the angular momentum vector and the eccentricity vector under third-body perturbations are free from singularities and provide a general, compact, and elegant way of presenting the variation of the mean elements.

When the mean angular momentum vector is scaled by L′, the variations of the mean vectorial elements h = G/L′=ηn and $\mathit{e}=e\widehat{\mathit{e}}$ take the neat symmetric form (Milankovitch 1941; Allan & Ward 1963; Allan & Cook 1964; Rosengren & Scheeres 2014)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{h}}{\mathrm{d}t}& =\mathit{h}\times \frac{\partial \mathcal{R}}{\partial \mathit{h}}+\mathit{e}\times \frac{\partial \mathcal{R}}{\partial \mathit{e}},\hfill \end{array}$$(56)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{e}}{\mathrm{d}t}& =\mathit{e}\times \frac{\partial \mathcal{R}}{\partial \mathit{h}}+\mathit{h}\times \frac{\partial \mathcal{R}}{\partial \mathit{e}},\hfill \end{array}$$(57)

in which

$$\begin{array}{c}\hfill \mathcal{R}=-\frac{1}{L^{\prime }}{⟨{\mathcal{V}}_{\star }⟩}_{\ell },\end{array}$$(58)

and ⟨𝒱_⋆⟩_ℓ is given in Eq. (25). Furthermore, derivation of Eqs. (56)–(57) is trivial after minor rearrangement of Eq. (32) to replace the directions $\widehat{\mathit{e}}$ and n by the nondimensional magnitudes e and h, respectively.

Thus, instead of replacing Eq. (31) in Eq. (32), the latter is written in the more convenient form

$$\begin{array}{cc}\hfill {⟨V_i⟩}_u=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right)\hfill \\ & \times {\mathrm{\Xi }}_{j+k-m}^{i-2l-k}{(-1)}^{j+l+m}{e}^{j+m-k}{(e\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^k{\left[{\eta }^2{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^2\right]}^{\frac{i-k}{2}-l}.\hfill \end{array}$$(59)

When now we make use of the identity

$$\begin{array}{c}\hfill {(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^2+{(\mathit{b}·{\widehat{\mathit{r}}}_{\star })}^2+{(\mathit{n}·{\widehat{\mathit{r}}}_{\star })}^2=1,\end{array}$$(60)

which gives the square of the modulus of the third-body direction $\Vert {\widehat{\mathit{r}}}_{\star }\Vert =1$ when it is computed in the apsidal frame, the dependence of the averaged potential on the unit vector b in the binormal direction is replaced by a corresponding dependence on the normal direction n. That is

$$\begin{array}{cc}\hfill {⟨V_i⟩}_u=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right)\hfill \\ & \times {\mathrm{\Xi }}_{j+k-m}^{i-2l-k}{(-1)}^{j+l+m}{(\mathit{e}·\mathit{e})}^{\frac{j+m-k}{2}}{(\mathit{e}·{\widehat{\mathit{r}}}_{\star })}^k{[X-{(\widehat{\mathit{e}}·{\widehat{\mathit{r}}}_{\star })}^2]}^{\frac{i-k}{2}-l},\hfill \end{array}$$(61)

with the scalar function

$$\begin{array}{c}\hfill X=X(\mathit{e},\mathit{h},{\widehat{\mathit{r}}}_{\star })\equiv 1-e^2+{\xi }^2-{\zeta }^2,\end{array}$$(62)

where we introduced the abbreviations $\xi =\mathit{e}·{\widehat{\mathit{r}}}_{\star }$, $\zeta =\mathit{h}·{\widehat{\mathit{r}}}_{\star }$.

Now, applying the binomial expansion again, we find

$$\begin{array}{cc}\hfill {⟨V_i⟩}_u=& {\displaystyle \sum _{l=0}^{\lfloor i/2\rfloor }}\left(\begin{array}{c}i\\ l\end{array}\right)\left(\begin{array}{c}2i-2l\\ i\end{array}\right){\displaystyle \sum _{j=0}^{2l+1}}\left(\begin{array}{c}2l+1\\ j\end{array}\right){\displaystyle \sum _{k=0}^{i-2l}}\left(\begin{array}{c}i-2l\\ k\end{array}\right){\displaystyle \sum _{m=0}^k}\left(\begin{array}{c}k\\ m\end{array}\right){\displaystyle \sum _{s=0}^{\frac{i-k}{2}-l}}\left(\begin{array}{c}\frac{i-k}{2}-l\\ s\end{array}\right)\hfill \\ & \times {\mathrm{\Xi }}_{j+k-m}^{i-2l-k}{(-1)}^{j+l+m+s}{(\mathit{e}·\mathit{e})}^{\frac{j+m-k}{2}-s}{(\mathit{e}·{\widehat{\mathit{r}}}_{\star })}^{2s+k}{X}^{\frac{i-k}{2}-l-s},\hfill \end{array}$$(63)

which only depends on e and h, as desired, as well as on the direction ${\widehat{\mathit{r}}}_{\star }$ of the disturbing body.

Explicit expressions for the first terms of the Legendre polynomial expansion yield the compact expressions

$$\begin{array}{cc}\hfill {⟨V_2⟩}_u& =1-3(2e^2-5{\xi }^2+{\zeta }^2),\hfill \end{array}$$(64)

$${\langle V_3\rangle }_u=-\frac{5}{2}\xi \left[3\left(1-8^2\right)-15{\zeta }^2+35{\xi }^2\right],$$(65)

$$\begin{array}{cc}\hfill {⟨V_4⟩}_u& =\frac{3}{4}\{3-20e^2+80e^4\hfill \\ \hfill & +10[7(1-10e^2){\xi }^2-(3-10e^2){\zeta }^2]\hfill \\ \hfill & +35(21{\xi }^4-14{\zeta }^2{\xi }^2+{\zeta }^4)\},\hfill \end{array}$$(66)

$$\begin{array}{cc}\hfill {⟨V_5⟩}_u& =-\frac{3}{4}\xi \{35(1-8e^2+40e^4)\hfill \\ \hfill & +490[(1-12e^2){\xi }^2-(1-4e^2){\zeta }^2]\hfill \\ \hfill & +147(33{\xi }^4-30{\zeta }^2{\xi }^2+5{\zeta }^4)\},\hfill \end{array}$$(67)

$$\begin{array}{cc}\hfill {⟨V_6⟩}_u& =\frac{1}{4}\{5(5-42e^2+168e^4-560e^6)\hfill \\ \hfill & +105[3(3-28e^2+168e^4){\xi }^2\hfill \\ \hfill & -(5-28e^2+56e^4){\zeta }^2]\hfill \\ \hfill & +315[(5-14e^2){\zeta }^4-18(3-14e^2){\zeta }^2{\xi }^2\hfill \\ \hfill & +33(1-14e^2){\xi }^4]\hfill \\ \hfill & +231(429{\xi }^6-495{\zeta }^2{\xi }^4+135{\zeta }^4{\xi }^2-5{\zeta }^6)\},\hfill \end{array}$$(68)

$$\begin{array}{cc}\hfill {⟨V_7⟩}_u& =-\frac{9}{16}\xi \{35(5-48e^2+224e^4-896e^6)\hfill \\ \hfill & +105[11(3-32e^2+224e^4){\xi }^2\hfill \\ \hfill & -3(15-96e^2+224e^4){\zeta }^2]\hfill \\ \hfill & +231[143(1-16e^2){\xi }^4-110(3-16e^2){\zeta }^2{\xi }^2\hfill \\ \hfill & +15(5-16e^2){\zeta }^4]\hfill \\ \hfill & +429(715{\xi }^6-1001{\zeta }^2{\xi }^4+385{\zeta }^4{\xi }^2-35{\zeta }^6)\},\hfill \end{array}$$(69)

$$\begin{array}{cc}\hfill {⟨V_8⟩}_u& =\frac{5}{64}\{7(35-360e^2+1728e^4-5376e^6+16128e^8)\hfill \\ \hfill & +252[11(5-54e^2+288e^4-1344e^6){\xi }^2\hfill \\ \hfill & -(35-270e^2+864e^4-1344e^6){\zeta }^2]\hfill \\ \hfill & +1386[143(1-12e^2+96e^4){\xi }^4\hfill \\ \hfill & -66(5-36e^2+96e^4){\zeta }^2{\xi }^2\hfill \\ \hfill & +(35-180e^2+288e^4){\zeta }^4]\hfill \\ \hfill & +12012[143(1-18e^2){\xi }^6\hfill \\ \hfill & +33(5-18e^2){\xi }^2{\zeta }^4-429(1-6e^2){\xi }^4{\zeta }^2\hfill \\ \hfill & -(7-18e^2){\zeta }^6]\hfill \\ \hfill & +6435(2431{\xi }^8-4004{\zeta }^2{\xi }^6\hfill \\ \hfill & +2002{\zeta }^4{\xi }^4-308{\zeta }^6{\xi }^2+7{\zeta }^8)\}.\hfill \end{array}$$(70)