Issue 
A&A
Volume 522, November 2010



Article Number  A60  
Number of page(s)  11  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201014496  
Published online  03 November 2010 
Explicit expansion of the threebody disturbing function for arbitrary eccentricities and inclinations
ASD, IMCCECNRS UMR8028, Observatoire de Paris,
UPMC, 77 avenue DenfertRochereau,
75014
Paris,
France
email: laskar@imcce.fr
Received:
24
March
2010
Accepted:
27
June
2010
Since the original work of Hansen and Tisserand in the XIXth century, there have been many variations in the analytical expansion of the threebody disturbing function in series of the semimajor axis ratio. With the increasing number of planetary systems of large eccentricity, these expansions are even more interesting as they allow us to obtain for the secular systems finite expressions that are valid for all eccentricities and inclinations. We revisited the derivation of the disturbing function in Legendre polynomial, with a special focus on the secular system. We provide here expressions of the disturbing function for the planar and spatial case at any order with respect to the ratio of the semimajor axes. Moreover, for orders in the ratio of semimajor axis up to ten in the planar case and five in the spatial case, we provide explicit expansions of the secular system, and simple algorithms with minimal computation to extend this to higher order, as well as the algorithms for the computation of non secular terms.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability / planetary systems / methods: analytical
© ESO, 2010
1. Introduction
In the threebody 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; AbuElAta & Chapront 1975; Laskar & Robutel 1995), but this approach, which is well suited to the study of the Solar System, has some limitations for some extrasolar 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 semimajor 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; AbuElAta & 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^{1}. 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 threebody system was expanded up to the octupolar order (α^{3}). 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 selfcontained 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 Nbody system, we consider here a threebody 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))
where
is the Keplerian interaction with elliptical elements a,e,i,M,ω, and Ω that denote, respectively, the semimajor 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_{1} part of the perturbation and indirect part T_{1} are then
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_{1}, 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
where P_{n}(x) are the Legendre polynomials that can be written as
with
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
and
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 [ ,
and, with n = p + q, and after changing q to n − q,
Thus
with
for 0 ≤ q ≤ n. If we write x = v − v′ + ω − ω′, ℱ_{n} can now be expressed in the form
with
The quantities ℱ_{n,q} can then be expressed in term of Hansen coefficients defined for n,m ∈ Z as
For convenience, we denote and . We have thus
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
we have
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)
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 α
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
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
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
with the same notations as Tisserand (1885)
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
where the are polynomials in μ,ν of degree n that are called the Tisserand polynomials^{2} 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_{20} 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}/γ′^{n + 1})F_{n} in term of Hansen coefficients. We thus have
with
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
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
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
Tisserand functions for the spatial case.
Tisserand functions for the spatial case 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 threebody 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; AbuElAta & Chapront 1975; Laskar & Robutel 1995)
with as before x = u − u′, y = u + u′, and
with c = cos(i/2), s = sin(i/2), and the same for the primes. We write
with
With these notations, we have
and
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
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_{1} (3) is more straightforward, as it compares with the computation of r·r′. Indeed, we have (Laskar & Robutel 1995, Eq. (24)),
where are the velocities in the corresponding Keplerian problem. With classical computation, we obtain
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
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)),
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^{2}(J/2), ν_{ * } = sin^{2}(J/2). In the planar case (Sect. 3), μ_{ * } = 1, ν_{ * } = 0.
We considered here heliocentric coordinates. One could also use Jacobi coordinates for the threebody 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 nonsecular inequalities.
6. Conclusion
We have presented a selfcontained exposition of the expansion of the threebody Hamiltonian in canonical heliocentric coordinates in term of Legendre polynomials at any order in the ratio of semimajor 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 α^{10} for the planar case, α^{5} for the spatial case expressed in terms of the mutual inclination, and α^{3} 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 nonsecular 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_{100}, 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_{50} 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.
One should consult (Aksenov 1986) for a detailed discussion of the relation between the Tisserand polynomials and the inclination functions of Kaula (1962).
Appendix A: Computation of the Hansen coefficients
is the polynomial part of the Hansen coefficients .
Hansen coefficients for (0 ≤ n ≤ 10,0 ≤ m ≤ n).
The Hansen coefficients are defined as the Fourier coefficients
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
where v is the true anomaly and E the eccentric anomaly, we have for n ≥ 2
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)
It is important to note that from (A.3), we have
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
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 β^{2}. Using a relation among hypergeometric series due to Gauss, and the change of variable e = 2β/(1 + β^{2}), Tisserand (1899) provides a finite expression of in term of hypergeometric series of e^{2}. 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
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
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
We can now prove by recurrence using the relation (A.7) that the general form of for n ≥ 0 and 0 ≤ m ≤ n is
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^{2}(J/2), ν = sin^{2}(J/2), we have
and with a straightforward expansion in complex notation, we obtain
which we can reorder as
with
With the expansion (8), we have then
with
We can then reorder (B.5) with s = k + r as
exchange the sum in p and k
and set q = p + s which gives
That is, with x = u − u′ and y = u + u′,
with
that is
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)
that is, with μ_{ * } = a + ia′, ν_{ * } = b + ib′,
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
with
The coefficients (B.12) and (B.16) in the Tisserand functions are very similar. Indeed, using μ = 1 − ν and (34), we have
with and
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
As F_{n} = P_{n}(z) (Eq. (8)) with z = cosx + ν(cosy − cosx) (Eq. (B.1)), we have from (B.19),
Using expressions (B.10, B.17), one then replaces F_{n} in (B.20) by
which leads to
where
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
or equivalently
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
Tisserand (1885) needs then some lengthy computations to determine from the coefficient of ν^{n} in F_{n}. Here, we use the expression (B.12) with (B.17) and (B.26). With μ = 1 − ν, we get
Calculating the coefficient of the term x^{s} in (1 + x)^{n − q}(1 + x)^{q} = (1 + x)^{n}, one finds
Thus, from (B.27), we obtain
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
where
and
In (B.31), we used the fact that is the complex conjugate of since F_{n} is real.
References
 AbuElAta, N., & Chapront, J. 1975, A&A, 38, 57 [NASA ADS] [Google Scholar]
 Aksenov, E. P. 1986, SvA, 30, 221 [NASA ADS] [Google Scholar]
 Blaes, O., Lee, M. H., & Socrates, A. 2002, ApJ, 578, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Brumberg, V. A. 1967, Bull. Inst. Theor. Astron., XI, 73 [Google Scholar]
 Brumberg, V. A., Evdokimova, L. S., & Kochina, N. G. 1971, Celest. Mech., 3, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Gastineau, M., & Laskar, J. 2009, TRIP 1.0, TRIP Reference manual, IMCCE, Paris Observatory, http://www.imcce.fr/Equipes/ASD/trip/trip.html [Google Scholar]
 Giacaglia, G. E. O. 1974, Celest. Mech., 9, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, P. 1855, Abhandl. d. K. S. Ges. d. Wissensch, IV, 182 [Google Scholar]
 Hill, G. W. 1875, The Analyst, II, 161 [CrossRef] [Google Scholar]
 Hughes, S. 1981, Celest. Mech., 25, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 1962, AJ, 67, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1959, AJ, 64, 367 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Krymolowski, Y., & Mazeh, T. 1999, MNRAS, 304, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Laplace, P. 1785, in Oeuvres Complètes, T11 (Paris: GauthierVillars) [Google Scholar]
 Laskar, J. 1990, in Les Méthodes Modernes de la Mécanique Céleste, ed. C. F. D. Benest [Google Scholar]
 Laskar, J. 2005, Celest. Mech. Dynam. Astron., 91, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Laskar, J., & Robutel, P. 1995, Celest Mech. Dynam. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Lidov, M. L., & Ziglin, S. L. 1976, Celest. Mech., 13, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Migaszewski, C., & Goździewski, K. 2008, MNRAS, 388, 789 [NASA ADS] [CrossRef] [Google Scholar]
 Tisserand, F. 1885, Annales de l’Observatoire de Paris, 18, 1 [Google Scholar]
 Tisserand, F. 1899, Traité de Mécanique Céleste (Paris: GauthierVillars) [Google Scholar]
 Whittaker, E., & Watson, G. 1927, A course of modern analysis (Cambridge University Press) [Google Scholar]
All Tables
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.