© ESO, 2010

1. Introduction

In the three-body problem, there are two classical ways to compute the principal part of the disturbing function  −Gmm′/Δ. The first approach is to expand it in powers of eccentricity and inclination, with coefficients that are expressed in term of Laplace coefficients (e.g.  Laplace 1785; Abu-El-Ata & Chapront 1975; Laskar & Robutel 1995), but this approach, which is well suited to the study of the Solar System, has some limitations for some extra-solar planetary systems, where the eccentricity can reach very high values. Another drawback, is that without some special truncation corrections, the angular momentum will not be conserved exactly in the truncated system.

In the second approach, the expansion is made with respect to the ratio of the semi-major axes of the two bodies α = a/a′, where a′ is related to the external body. Although it may be less efficient for low eccentricity and large values of α, as for the inner planets in the Solar system, the advantage is that this expansion allows us to obtain finite expressions for arbitrary eccentricities in the secular system, while both approaches allow expansions for arbitrary inclinations. The most important contribution to this problem was made in the XIXth century (Hansen 1855; Hill 1875; Tisserand 1899). With the development of space technology, expansions of the disturbing function in term of Legendre polynomials have been rejuvenated for the construction of satellite theories in the vicinity of the Earth (Kozai 1959; Kaula 1962; Brumberg 1967; Brumberg et al. 1971; Giacaglia 1974; Abu-El-Ata & Chapront 1975).

The discovery of new planetary systems, starting with 51Pegb (Mayor & Queloz 1995), has raised the need to revise these methods as many systems have planets with very high eccentricity, as GL 581, HD 217107, HD 69830, HD 74156, HD 168443, HD 102272, HD 169830, HD 202206, HD 183263, or even HD 80606, where the eccentricity reaches 0.93¹. With the discovery of these numerous new planetary systems, the previous analytical expansions in Laplace coefficients developed for the Solar System are no longer the most appropriate, and we are faced with the need to understand more globally the dynamics of these systems. On the other hand, as the parameters of these extrasolar planetary systems remain not very well known, there is no necessity for very precise analytical approximations. This has led several authors to use the Legendre expansion of the potential for the dynamical study of the secular planetary system, following the previous studies of stellar systems (Krymolowski & Mazeh 1999; Ford et al. 2000; Blaes et al. 2002), where the secular spatial three-body system was expanded up to the octupolar order (α³). In particular, Lee & Peale (2003) used the secular planar system at octupolar order to study the secular dynamics of the HD 168443 and HD 12661 systems. Migaszewski & Goździewski (2008) computed the planar secular system to high order using computer algebra to average over the mean anomaly.

As it appears that there is a growing interest in analytical studies of extra solar planetary systems, we considered it interesting to present a derivation of the planetary (or stellar) disturbing function in a very simple and explicit way that does not already appear in the existing literature. We aimed to write a self-contained paper that allows one to construct explicitly the planetary Hamiltonian to high order, with minimal additional computation. In particular, in Sect. 3 we show that the planar secular system can be obtained explicitly at any order without the need for computer algebra. The spatial case is treated in Sect. 4 with minimal computation when expressed with respect to the mutual inclination J. In the presence of more than two planets, it may be more convenient to use a fixed reference frame, and the derivation of the spatial Hamiltonian is provided in this case in Sect. 4.2. Our derivation is close to the original formulations of Hansen (1855) and Tisserand (1899), but with emphasis on the secular system and the direct derivation of explicit expressions. We also present in Appendix A a new derivation of the Hansen coefficients for the secular terms.

2. Expansion of the disturbing function

To simplify the notations, and although everything can be generalized to an N-body system, we consider here a three-body problem with a central body of mass M and two other celestial bodies of masses m,m′. If we write the Hamiltonian of the Newtonian interactions among these three bodies in Poincaré canonical heliocentric variables, we obtain (Laskar & Robutel 1995, Eq. (15))

(1)

where

(2)

is the Keplerian interaction with elliptical elements a,e,i,M,ω, and Ω that denote, respectively, the semi-major axis, eccentricity, inclination, mean anomaly, argument of perihelion, and longitude of the node (with primes for the external body of mass m′). The principal part U₁ part of the perturbation and indirect part T₁ are then

(3)

where r,r′ are the radius vectors of the inner and outer planets, with norms r,r′, unit vectors u = r/r and u′ = r′/r′, and conjugate momentum . We focus first on the principal part of the Hamiltonian U₁, which is the most difficult part to compute, while the computation of the indirect part will be made in Sect. 5. With ρ = r/r′, and , we have

(4)

where P_n(x) are the Legendre polynomials that can be written as

(5)

with

(6)

With α = a/a′,γ = r/a,γ′ = r′/a′, we have ρ = αγ/γ′, and thus, as γ,γ′,u,u′ do not depend on a,a′, the expansion of ℱ in powers of α is

(7)

and

(8)

3. Planar case

We first study the planar case that leads to some simplifications in the expansions. We define v,v′ to be the true anomalies of u,u′, and ω,ω′ their argument of perihelion, u = v + ω, u′ = v′ + ω′, and x = u − u′. We have then u·u′ = cosx and, following a classical computation (e.g. Whittaker & Watson 1927, p. 303), we have for all z ∈  [ 0,1 [ ,

(9)

and, with n = p + q, and after changing q to n − q,

(10)

Thus

(11)

with

(12)

for 0 ≤ q ≤ n. If we write x = v − v′ + ω − ω′, ℱ_n can now be expressed in the form

(13)

with

(14)

The quantities ℱ_n,q can then be expressed in term of Hansen coefficients defined for n,m ∈ Z as

(15)

For convenience, we denote and . We have thus

(16)

For arbitrary values of k ∈ Z, the Hansen coefficients can be expressed in an explicit manner as an infinite series involving Bessel functions and hypergeometric functions (Hansen 1855; Tisserand 1899), or in terms of generalized Laplace coefficients (Laskar 2005), but for k = 0, reduces to a finite polynomial expression in e, 1/e, , and (see Appendix A). We have thus a very compact expression for the coefficient of any argument kM + k′M′ in explicit form at all orders n in α. Indeed, if we denote this coefficient , that is

(17)

we have

(18)

In particular, for all n ∈ N, the secular part can even be simplified, using the classical relation among Hansen coefficients and the relation f_n,q = f_n,n − q from Eq. (12)

(19)

where ϵ_n = 0 if n is odd, and ϵ_n = 1 if n is even.

3.1. Practical algorithm

Equations (12) and (19) provide an explicit algorithm for the computation of the principal part of the disturbing function ℱ, and in particular for the computation of the secular Hamiltonian at any order N in α

(20)

for which the involved Hansen coefficients reduce to finite expressions (see Appendix A). As there are no secular terms in the indirect part of the Hamiltonian (see Laskar & Robutel 1995), the computation of the secular Hamiltonian reduces to the computation of the rational constants f_n,k given by Eq. (12). In practice, for finite order n, it is even easier to use the expression of F_n (11) (Table 1) and to translate it into ℱ_n (19) by the simple transformation

(21)

which allows an immediate computation of any argument in terms of Hansen coefficients.

3.2. Computation of the secular part

The computation of the secular part of the Hamiltonian is quite simple, as we will just have to make the transformation

(22)

in the expression of F_n. Moreover, we have for m ≥ n ≥ 1 (see Appendix A). The term in cosnx can thus be discarded in the expression of F_n (Table 1) which simplifies the expression of the Hamiltonian. The secular Hamiltonian is thus expressed in finite form, using the values of the Hansen coefficients given in Appendix A. We have

These results are equivalent to the expression obtained by Migaszewski & Goździewski (2008) using computer algebra.

4. Spatial case

The spatial case is more complicated as it involves additional variables. Our goal is to derive explicit formulae that are as compact as possible. We thus expand in terms of the mutual inclination J. For each orbit, we use a reference frame (i,j,k) associated with the orbit, with first vector i in the direction of the ascending node of r′ over r. With u = v + ω and u′ = v′ + ω′, we have

(23)

with the same notations as Tisserand (1885)

(24)

As in the planar case, we have F_n = P_n(u·u′), but now u·u′ is given by the slightly more complex expression (23). For all n, we have

(25)

where the are polynomials in μ,ν of degree n that are called the Tisserand polynomials² as a recognition of the work of Tisserand (1885), although these expressions are already present in (Hansen 1855). As (u,u′) −  → (− u, − u′) leaves u·u′ unchanged, we have , thus F_n can be expressed as a trigonometric polynomial in cos(mu + m′u′). Although they can be computed explicitly for all n (see Appendix B), for a given n, it is often more efficient to make a direct computation of F_n on a computer algebra system. For example, F₂₀ is computed in less than 1.5 s in exact rational arithmetics with TRIP (Gastineau & Laskar 2009) on an average laptop computer using the simple expression (8). We can then express ℱ_n = (γⁿ/γ′^n + 1)F_n in term of Hansen coefficients. We thus have

(26)

with

(27)

More practically, starting from the finite expressions in cosine polynomials (Table 2) of the Tisserand functions F_n, the expression of ℱ_n is obtained as in the planar case by means of the more general transformation

(28)

4.1. Computation of the secular part

As in the planar case, the computation of the secular part of the Hamiltonian is just a straightforward translation of F_n in Table 2, with the transformation

(29)

Moreover, as for n ≥ 1, all terms in cos(mu ± nu′) can be discarded in F_n. We thus have,

We thus observe here an important simplification in the quadrupolar secular Hamiltonian , as all terms involving the external planet longitude of perihelion ω′ vanish and we are left with an integrable Hamiltonian. This is what Lidov & Ziglin (1976) called a happy coincidence (see also Farago & Laskar 2010). This is no longer the case at higher orders.

4.2. Expression in a fixed reference frame

The above expressions are given with respect to the mutual inclination to shorten the algebraic expansions. This is especially useful in a three-body problem, but to obtain expressions valid in a fixed reference frame, then one needs to substitute into u·u′ its expression in terms of the elliptical elements of the two bodies. We can then generalize the expression (23) and write it now as (Brumberg 1967; Abu-El-Ata & Chapront 1975; Laskar & Robutel 1995)

(30)

with as before x = u − u′, y = u + u′, and

(31)

with c = cos(i/2), s = sin(i/2), and the same for the primes. We write

(32)

with

(33)

With these notations, we have

(34)

and

(35)

The remaining part is identical to the previous case. The expressions of the Tisserand functions F_n are given for 0 ≤ n ≤ 3 in Table 3 and for all n in Appendix B. In practice, for small values of n we can use explicit expressions of F_n, using the straightforward approach consisting of computing F_n directly with computer algebra, using the relation (8), and then to translate it as previously in order to obtain the expression of any argument kM + k′M′ by means of the relations

(36)

We thus have for the secular part ,

and

5. Indirect part

Until now, we considered only the principal part of the perturbing Hamiltonian. The computation of the indirect part T₁ (3) is more straightforward, as it compares with the computation of r·r′. Indeed, we have (Laskar & Robutel 1995, Eq. (24)),

(37)

where are the velocities in the corresponding Keplerian problem. With classical computation, we obtain

(38)

where E is the eccentric anomaly and J_k(x) are the Bessel functions. The coordinates (Ẋ,Ẏ) of the velocity in the reference frame of the orbit with origin at perihelion are then easily expressed in Fourier series of the mean anomaly, as

(39)

If we denote 𝒵 = Ẋ + iẎ, we then have for the spatial problem in the fixed reference frame, with the same notations as in the previous section (Laskar & Robutel 1995, Eq. (37)),

(40)

If one considers the mutual inclination J, as in Sect. 4, this expression simplifies since Ω = Ω′, s′ = 0, c′ = 1, and μ _* ,ν _*  are then real, with μ _*  = cos²(J/2), ν _*  = sin²(J/2). In the planar case (Sect. 3), μ _*  = 1, ν _*  = 0.

We considered here heliocentric coordinates. One could also use Jacobi coordinates for the three-body problem. In this case, the indirect part does not require additional computations as it is expressed in term of u·u′ (see Laskar 1990). It should be noted that in terms of both heliocentric coordinates or Jacobi coordinates, the indirect part does not contribute to the secular system. We have provided here the expression of the indirect part for the computation of non-secular inequalities.

6. Conclusion

We have presented a self-contained exposition of the expansion of the three-body Hamiltonian in canonical heliocentric coordinates in term of Legendre polynomials at any order in the ratio of semi-major axes α, which can also be adapted in the case of Jacobi coordinates, where only the indirect part differs. We have included here all the necessary material that allows one to write explicitly the secular Hamiltonian at order α¹⁰ for the planar case, α⁵ for the spatial case expressed in terms of the mutual inclination, and α³ for the spatial case in a fixed reference frame. With the additional computation of the required Hansen coefficients, the expressions of the Tisserand functions F_n can also be used for a straightforward computation of the expression of a non-secular inequality kM + k′M′.

As the algorithms that are presented here are very simple, we have not added to this paper any tables in electronic form. Indeed, the reader who needs to use expressions of high order that do not appear in the paper, will have no problems in exploiting the algorithms given here to derive the required expressions. The written results of the paper can then be used to check his computations for the lowest orders. For example, the computation of F₁₀₀, in the spatial case takes less than 8 min on an average 8 core desktop computer in exact rational arithmetics using TRIP, but has 2 343 926 terms, while F₅₀ needs only 7.5 s with already 164 151 terms. As the algorithms require only a few lines of code (fewer than 20 in our case), one understands that it is preferable to compute the terms when they are needed than to store the values in electronic form.

Appendix A: Computation of the Hansen coefficients

The Hansen coefficients are defined as the Fourier coefficients

(A.1)

For arbitrary values of k ∈ Z, the Hansen coefficients can be expressed in an explicit manner as an infinite series involving Bessel functions and hypergeometric functions (Hansen 1855; Hill 1875; Tisserand 1899), or generalized Laplace coefficients (Laskar 2005), but for k = 0, reduces to a simple polynomial in e, 1/e, , and . The literature on the computation of the Hansen coefficients has been huge since the original work of Hansen (1855), but we do not review it here. We concentrate on the obtention of the coefficients and , for n ≥ 0 and 0 ≤ m ≤ n that are required to compute the averaged planetary or lunar Hamiltonian. Using

(A.2)

where v is the true anomaly and E the eccentric anomaly, we have for n ≥ 2

(A.3)

We still need to consider the case n = 1 for which the expansion in true anomaly v is not suitable. We have immediately, using (A.6)

(A.4)

It is important to note that from (A.3), we have

(A.5)

The Hansen coefficients are given in Table A.1 for n ≥ 2 and m ≤ n − 2. On the other hand, the computation of for n ≥ 0 is not as straightforward, as a direct expansion of

(A.6)

in eccentric anomaly is far more complicated than the previous expression (A.3). One can still perform some formal expansion, but this would not be very useful, as it would contain expressions with many summations that cannot be easily reduced to a single summation as in (A.3). In fact, the only explicit computation of available in the literature, was obtained through a complex process in (Hansen 1855; Hill 1875; Tisserand 1899), using the auxiliary variable . In this case, it can be shown that is expressed as a finite hypergeometric series in β². Using a relation among hypergeometric series due to Gauss, and the change of variable e = 2β/(1 + β²), Tisserand (1899) provides a finite expression of in term of hypergeometric series of e². We present here a more direct demonstration for the obtention of an explicit expression of which is a recurrence using the relation among Hansen coefficients

(A.7)

This relation appears in (Hughes 1981), but it was shown by Laskar (2005) that it is equivalent to a recurrence relation obtained on Laplace coefficients by Laplace (1785). The computation of for n ≥ 0 using (A.6) is straightforward and gives

(A.8)

The computation of can also be performed using (A.6). Changing E to  − E, it is immediate to see that the part in sinE cancel, and we are left with

(A.9)

We can now prove by recurrence using the relation (A.7) that the general form of for n ≥ 0 and 0 ≤ m ≤ n is

(A.10)

The computation is delicate but straightforward. One can first show that the two first elements (for l = 0) of the polynomial expressions of and in (A.7) cancel. We then change the index from l to l + 1 and show that the general term of the sum in (A.7) gives the general term of , and that the last term of will give the last term of . The expression (A.10) is equivalent to the expression of (Tisserand 1899).The expressions of the Hansen coefficients are given in Table A.2 for 0 ≤ n ≤ 10 and 0 ≤ m ≤ n.

Appendix B: Computation of Tisserand polynomials and Tisserand functions

Expansions valid for all inclinations were already given by Hansen (1855), and Tisserand (1885) in his researches on asteroidal motions, but these derivations are rather complex. Although in (B.3) we present some of Tisserand’s computations, we make first, as in the planar case some direct expansions and reordering of terms to obtain expressions that are as compact as possible.

B.1. Tisserand functions with respect to mutual inclination

With x = u − u′, y = u + u′, μ = cos²(J/2), ν = sin²(J/2), we have

(B.1)

and with a straightforward expansion in complex notation, we obtain

(B.2)

which we can reorder as

(B.3)

with

(B.4)

With the expansion (8), we have then

(B.5)

with

(B.6)

We can then reorder (B.5) with s = k + r as

(B.7)

exchange the sum in p and k

(B.8)

and set q = p + s which gives

(B.9)

That is, with x = u − u′ and y = u + u′,

(B.10)

with

(B.11)

that is

(B.12)

B.2. Tisserand functions with respect to a fixed reference frame

In the case of a fixed reference frame (Sect. 4.2), the expressions differ slightly, but are more complicated. We have (35)

(B.13)

that is, with μ _*  = a + ia′, ν _*  = b + ib′,

(B.14)

This is similar to the previous expression (B.1), and the computation of the Tisserand functions F_n will follow the same lines and gives the more general expression

(B.15)

with

(B.16)

The coefficients (B.12) and (B.16) in the Tisserand functions are very similar. Indeed, using μ = 1 − ν and (34), we have

(B.17)

with and

(B.18)

B.3. Tisserand simplification

The above expression (B.18) is a double sum. In the mutual inclination case, Tisserand (1885) could reduce it to a single sum using hypergeometric functions. He first considers the differential equation satisfied by Legendre polynomials

(B.19)

As F_n = P_n(z) (Eq. (8)) with z = cosx + ν(cosy − cosx) (Eq. (B.1)), we have from (B.19),

(B.20)

Using expressions (B.10, B.17), one then replaces F_n in (B.20) by

(B.21)

which leads to

(B.22)

where

(B.23)

The equality (B.22) must be satisfied for all x and y, thus all the coefficients (B.23) are equal to 0. The solutions of (B.23) are

(B.24)

or equivalently

(B.25)

where is an unknown constant and F an hypergeometric function (e.g. Whittaker & Watson 1927). Let us assume that the quantity 2n − 2q − 2s + 1 is positive. If it is not the case, one can make the change of variable (s,q) → (n − s′,n − q′). Then 2n − 2q′ − 2s′ + 1 is positive and since F_n is real, (B.10). From (B.25), one thus has

(B.26)

Tisserand (1885) needs then some lengthy computations to determine from the coefficient of νⁿ in F_n. Here, we use the expression (B.12) with (B.17) and (B.26). With μ = 1 − ν, we get

(B.27)

Calculating the coefficient of the term x^s in (1 + x)^n − q(1 + x)^q = (1 + x)ⁿ, one finds

(B.28)

Thus, from (B.27), we obtain

(B.29)

Finally, the most general case, where inclinations are defined with respect to a fixed reference plane, can be as well derived from (B.17) and (B.29). We have

(B.30)

where

(B.31)

and

(B.32)

In (B.31), we used the fact that is the complex conjugate of since F_n is real.

References

All Tables