Frequency analysis and representation of slowly diffusing planetary solutions

Context . Over short time-intervals, planetary ephemerides have traditionally been represented in analytical form as ﬁnite sums of periodic terms or sums of Poisson terms that are periodic terms with polynomial amplitudes. This representation is not well adapted for the evolution of planetary orbits in the solar system over million of years which present drifts in their main frequencies as a result of the chaotic nature of their dynamics. Aims . We aim to develop a numerical algorithm for slowly di ﬀ using solutions of a perturbed integrable Hamiltonian system that will apply for the representation of chaotic planetary motions with varying frequencies. Methods . By simple analytical considerations, we ﬁrst argue that it is possible to exactly recover a single varying frequency. Then, a function basis involving time-dependent fundamental frequencies is formulated in a semi-analytical way. Finally, starting from a numerical solution, a recursive algorithm is used to numerically decompose the solution into the signiﬁcant elements of the function basis. Results . Simple examples show that this algorithm can be used to give compact representations of di ﬀ erent types of slowly di ﬀ using solutions. As a test example, we show that this algorithm can be successfully applied to obtain a very compact approximation of the La2004 solution of the orbital motion of the Earth over 40Myr ([ − 35Myr, 5Myr]). This example was chosen because this solution is widely used in the reconstruction of the past climates.


Introduction
Before the computer age, long-term solutions for planetary orbits were derived by perturbation methods, and were obtained in the form of the sum of periodic terms.The first such solution was obtained by Lagrange (1782) and was later on improved by LeVerrier (1840LeVerrier ( , 1841)), who took Uranus into account as well.These long-term solutions have been found to be of fundamental importance for the understanding of the past climate of the Earth, when it was understood that the changes of the orbit of the Earth also induce some change in its obliquity and in the insolation at the Earth surface (Milankovitch 1941; for a detailed review, see Laskar et al. 2004).With the advent of computers, two different approaches become possible.Computer algebra allowed extending the perturbation methods (e.g., Bretagnon 1974), and as computer speed increased, direct integrations of the full solar system become possible (Quinn et al. 1991;Sussman & Wisdom 1992).Instead, Laskar (1985Laskar ( , 1986Laskar ( , 1988) developed a mixed strategy in which an analytical averaging of the planetary equations by perturbation methods was obtained with dedicated computer algebra, followed by numerical integration of the averaged system.In order to compare the output of the numerical integrations with the quasi-periodic solutions of the perturbative meth-Full Tables A.3 and A.4 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http: //cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/628/A84 ods, Laskar (1988) introduced the frequency analysis method that allows obtaining in a very efficient way a precise approximation of the numerical solution in quasiperiodic form (e.g., Laskar 2003).
One outcome of these computations was also to demonstrate that the solar system motion is chaotic (Laskar 1989(Laskar , 1990)).As a consequence, the solutions are not quasi-periodic, although they can be approximated by a quasi-periodic expression over a limited time of a few million years (Laskar 1988(Laskar , 1990;;Laskar et al. 2004).
We here derive a more adapted strategy for the slowly diffusing trajectories of a dynamical system that is very well suited to the construction of compact forms for the long-time behavior of planetary orbits.We introduce an algorithm that is derived from the frequency analysis to which a slow variation of the frequencies and amplitudes is added.By slowly diffusing solution, we mean a solution that while it experiences significant frequency drifts during the entire considered time interval is nearly quasiperiodic in time subintervals.Similar to the frequency analysis, a key step of our algorithm is to construct a frequency-dependent function basis on which the considered solution is decomposed.
After recalling some fundamental results of the frequency analysis, we consider in Sect. 2 a single-term model and show that it is possible to exactly recover its varying frequency as a function of time.This is important when the details of how a solution diffuses in the frequency space are of interest.To construct a compact solution representation with a reasonably high precision, however, small flickers without significant cumulative effects can be smoothed out from a varying frequency.It is therefore preferable to use a model with a limited number of parameters, for instance, low-order polynomials, to approximate the frequency.To do so, we sample the averaged frequencies over a sliding time-interval.For a slowly diffusing solution, it is assumed that all fundamental frequencies can be sampled in this way using the algorithm called numerical analysis of fundamental frequencies (NAFF; e.g., Laskar 2003).
We design a general algorithm that represents a slowly diffusing solution in Sect.3, where a frequency-dependent function basis is constructed based on the Chebyshev approximations of fundamental frequencies.The basis functions with significant but non-fundamental frequencies are generated according to the assumption that at any instant, a slowly diffusing solution is nearly quasi-periodic with the smoothed fundamental frequencies at that instant: all main frequencies are integral linear combinations of the fundamental frequencies.This effectively avoids the difficulties of determining the frequencies of longperiod terms and/or groups of neighboring-period terms.
Two simple examples are given in Sect. 4. In the first example, a weakly dissipated system is considered, and the represented solution diffuses because of dissipation.In the second example, a Hamiltonian system of degree 1.5 is considered, and a solution starting from an obvious resonance overlap zone is represented.To show the flexibility of our algorithm in representing such a solution, we exclude all possible libration frequencies from the set of fundamental frequencies used in constructing the function basis.
Applications to the representation of ephemerides of the solar system bodies are provided in Sect. 5.As examples, we represent the eccentricity and inclination of the Earth as given by the long-term numerical solution La2004 (Laskar et al. 2004).For such a realistic solution, it is natural to take restrictions on the representation model from previous results into account, such as the important libration frequencies from numerical analysis (e.g., Laskar et al. 2004), so that the resulting representations can be as similar as possible to the physical model.

Frequency analysis and time-dependent frequencies
For a KAM solution (after Kolmogorov, Arnold, Moser) of a dynamical system (Kolmogorov 1954;Arnol'd 1963;Moser 1962), with fundamental frequency vector ν = (ν 1 , . . ., ν N ) ∈ R N , where N is the number of degrees of freedom of the system, the motion of any given degree of freedom variable z(t) can be described in complex form as where k = (k 1 , . . ., k N ) is the frequency index vector, a k complex amplitude, and •, • denotes the usual inner product of two real vectors.
The NAFF algorithm was designed to numerically recover significant terms of Eq. (1) from a numerical sample set of z(t), for example, an ephemeris obtained by numerical integration (Laskar 1988).However, the application of NAFF is not limited to numerically recovering KAM solutions.Actually, most of the relevant works concern frequency drifts.In particular, the variations in fundamental frequencies of the secular solar system over 200 Myr was used to study the chaotic behavior of the solar system (Laskar 1990).
We first briefly recall the theorem about the convergence of the NAFF algorithm (for the complete version with proof, see Laskar 2003).With a set of N appropriate variables, the complex function (1) describing a degree of freedom of an analytic KAM solution can be written as where ν 1 is a fundamental frequency.We have then the following (e.g.Laskar 1999).
Theorem 1.Let ν T 1 be the value of σ ∈ R that maximizes the function φ(σ) = | z(t), e iσt χ T |, where χ = χ(t/T ) > 0 is a weight function, and •, • χ T the inner product of two complex functions of t defined as We then have lim T →∞ ν T 1 = ν 1 .Similar to the fundamental frequency ν 1 , any other main frequency can be recovered by searching its neighborhood or R for the value of σ that maximizes φ(σ), defined with the remaining z(t), that is, the original z(t) minus the recovered terms.
A diffusion solution is then characterized by fundamental frequency variations (Laskar et al. 1992;Laskar 1993).To determine a varying frequency, we first consider the following simplest case, where ν(t) and a(t) > 0 are real integrable functions.Writing where θ(t) = t 0 (ν(τ) − σ(τ))dτ belongs to C 0 , the set of continuous function on R, we have the theorem below.
Proof.The right-hand side of Eq. ( 5) can be written as Using the Cauchy-Bunyakovsky-Schwarz inequality1 , we can easily deduce from this formula that where • denotes the norm induced by the inner product (3).Moreover, the equality is true if and only if √ a(t) and √ a(t)e −iθ(t) are linearly dependent, that is, when θ(t) is a constant.This constant is 0, since θ(0) = 0 by definition.The above arguments imply that φ(σ(t)) has a unique maximum attained at θ(t) ≡ 0. Because ν(t) is continuous and the searched σ(t) is also continuous, θ(t) ≡ 0 is equivalent to σ(t) ≡ ν(t).We therefore conclude for any given T > 0 that φ(σ(t)) has one and only one maximum attained at σ(t) ≡ ν(t) .
Theorem 2 implies that even in the case of a varying frequency, it is still possible to recover the frequency exactly, as a function of time.Now, we consider the following more general complex function: where a 1 (a k ) ∈ S a with S a a subspace of real function space, θ k ∈ R, and {ν n } N n=1 ∈ S ν with S ν a linear subspace of the real integrable function space.As an extension of Eq. ( 2), this function inherits an important time-varying character of Eq. ( 2), that is, all phase increments are described by a single fundamental frequency vector ν(t).We make the heuristic assumption that under suitable conditions, a similar result as for Theorem 1 holds for a more general expression with varying frequencies as in Eq. ( 6).

Assumption. If ν T
1 − ν 1 is sufficiently small, where ν T 1 ∈ S ν and • is an appropriately defined T -dependent norm (e.g., the one induced by the inner product (3)), and

Representation of slowly diffusing solutions
From now on, we restrict ourselves to slowly diffusing solutions of an ordinary differential equation system obtained by slightly perturbing an integrable Hamiltonian system.We denote by (I, θ) = {I n , θ n } N n=1 the action-angle variables of the integrable Hamiltonian system, which is used to express the ephemeris of a diffusing solution {z n (t) ≡ I n (t)e iθ n (t) } N n=1 .

Representation procedure
We start from a sample set of {z n (t)} N n=1 .We assume that the samples are given at grid points from t 0 to t 1 with fixed time step h, which is much smaller than the minimum of the fundamental periods.In accordance with the condition of slow diffusion, we assume that the solution is close to quasi-periodic in the time subintervals (of [t 0 , t 1 ]) described below, with a length far exceeding the maximum of the fundamental periods.These roughly stated preconditions are required because we use NAFF algorithm to estimate the changing fundamental frequencies.
As the first step, we apply the NAFF algorithm to obtain fundamental frequency samples.For each given degree (n), this algorithm is applied to z n (t) samples over evenly spaced timesubintervals . For each of these subintervals, the first recovered frequency is taken as the averaged value of the fundamental frequency ν n over the same subinterval.The N averaged fundamental frequencies obtained in this way are then taken as the instant frequencies at τ λ , resulting in the fundamental frequency samples {τ λ , ν n (τ λ )} Λ λ=1 , n = 1, . . ., N. Second, we fit for each given degree the frequency samples to a Chebyshev expansion valid on [τ 1 , τ Λ ], where and T m (x) is the Chebyshev polynomial2 of degree m.We then numerically construct a frequency-dependent function basis B (see the next subsection), on which the considered solution is decomposed.
The final step, that is, decomposing z n (t) on B, is the same as that of the NAFF algorithm (Laskar 1999), except for the searched function bases for significant terms.The function bases are {e iω k t , ω k ∈ R} for NAFF and B for the procedure described here.

Representation model
By variation of parameters, Eq. ( 1) becomes Now, let be a truncated set of the frequency index vector k.For each k ∈ ẐN , let where a l,k ∈ C.These lead us from Eq. ( 8) to the following representation model of z v (t) : From the Chebyshev approximations of the fundamental frequencies, it is easy to obtain by integration the Chebyshev expansion representing the phase increment associated with the frequency k, ν(τ) , from the value, as is preferred, at the middle of the time interval where ϕ 0k ∈ R, and ϕ k (t) gathers all Chebyshev polynomials of degree larger than 0. With ϕ 0k and ϕ k (t), Eq. ( 11) can be written as where a l,k = a l,k e iϕ 0k ∈ C is simply the coordinate of z v (t) when it is decomposed on the function basis For the obtained representation, it should be noted that although we decompose a solution on the basis B, and correspondingly, express its representation by Eq. ( 13), ϕ k (t) is not used to specify the representation.{ν n (t)} N n=1 is used instead.In other words, the representation is specified by the coordinates a l,k , and the Chebyshev coefficients of {ν n (t)} N n=1 .This choice is made for the following two reasons.One is that the representation could otherwise be unnecessarily cumbersome.The other reason is that from the Chebyshev coefficients of {ν n (t)} N n=1 , those of ϕ k (t) can be easily obtained according to Eq. ( 12).Here, we also recall that the solution representation, because of its dependence on the Chebyshev approximations of ν n (t) and a k (t), is valid on [τ 1 , τ Λ ] rather than [t 0 , t 1 ].
To complete the description of the representation model, we point out that the integers {M n , K n } N n=1 and {L k } k∈ ẐN , introduced in Eqs. ( 7), ( 9), and (10), respectively, are necessary parameters for defining a particular representation procedure.Of course, the number of these parameters can be reduced by requiring that some or all of these parameters take the same sufficiently high value, for instance, L k = L for all k ∈ ẐN , at the price of unnecessarily increasing the basis dimension.Another important point to mention is that to terminate the procedure, the maximum number of representation terms (J) and/or the required precision has to be set beforehand.The required precision is specified by using either the absolute truncation error δ or the relative truncation error δ r .Correspondingly, the procedure is terminated if where the module • is induced by the inner product (3) with prescribed χ ≡ 1.In the following, these parameters, as summarized in Table 1, are referred to as procedure parameters.

Examples
In the following two subsections, we illustrate our representation procedure with two slowly diffusing solutions, of which one diffuses as a result of a dissipative perturbation, and the other as a result of its chaotic nature.In order to illustrate the convergence property of this procedure, it is convenient to write the resulting representation with a single term index, that is, where the terms are arranged in the same order as they are obtained in the procedure, which roughly corresponds to the order of decreasing |a j |.

Example of a dissipated solution
Consider the following weakly dissipated system:  18)-( 20).Also shown is the phase orbit of the unperturbed pendulum system, namely (Eq.( 18)) with ε = 0, which passes through the initial phase point of the dissipated solution.The dashed line depicts the separatrix of the unperturbed pendulum.where If we nullify the dissipative term −ε dθ dt in Eq. ( 18), the system is a simple pendulum with Hamiltonian H 0 (I) = I 2 2 + cos(θ), where I = θ.Take as an example the solution z(t) ≡ I(t)e i θ(t) of Eq. ( 18) with the following initial conditions: Because of the dissipative perturbation, the phase orbit starting from (θ 0 , I 0 ) gradually decays away from the unperturbed orbit passing through the same phase point.This dissipation effect is significant in the long run, as is shown in Fig. 1.As a solution of a system of 1 degree of freedom with a dissipative perturbation, the decaying z(t) has a changing fundamental frequency ν, which inherits the characteristic frequency of the unperturbed system and varies as the solution decays.25), a representation of the dissipated solution z(t) specified by Eqs. ( 18)-( 20), which is obtained with {M, K, L, δ r } = {9, 10, 9, 10 −5 }.

Practical implementation
By numerical integration, an ephemeris of z(t) is obtained at {t i = ih : i = 0, . . ., 65535, h = 0.15}.We then apply NAFF to 129 evenly spaced time subintervals with length d = 2048 h and midpoints {τ λ = d 2 + 496(λ − 1)h} 129 λ=1 , respectively, to obtain a sample set of the changing frequency, {τ λ , ν(τ λ )} 129 λ=1 .This frequency sample set is shown in Fig. 2. We also show a Chebyshev expansion approximating ν(t).This expansion is of degree 9 and is expressed as where the Chebyshev polynomials T m (x) are listed in The values of the coefficients {c m } 9 m=0 are given in Table 2.We write the Chebyshev expansion approximating the phase increment, from the value at x = 0 or t = τ 1 +τ 129 2 , associated with ν in two parts, where ϕ 0 is a real constant and gathers all Chebyshev polynomials of degree larger than 0. In Eq. ( 23), τ 129 −τ 1 2 is the time-scaling factor, and 3 , with This completes the determination of ϕ(x), which is used in the following.Although the constant ϕ 0 is not used in the following, it can be easily determined as ϕ 0 = −ϕ(0) because the phase increment ( 22) is 0 at t = τ 1 +τ 129 2 or x = 0.With procedure parameters {M, K, L, δ r } = {9, 10, 9, 10 −5 }, we obtain a 37-term representation, 3 For the general case, see Eq. (A.7).26) and ( 27).In the chaotic zone formed by resonance overlap lies the phase point (θ 0 , I 0 ) = (0, 1.535), which will be chosen as the initial phase point of the considered chaotic solution.26)-( 28).The samples of this frequency are computed by applying NAFF algorithm over a sliding time-interval, and the errors are estimated as their respective differences from the frequencies of the quasi-periodic approximation of the solution over the same sliding time interval.Also shown is a Chebyshev appoximation of ν 3 (t).
of the solution z(t), where T l( j) (x) is the Chebyshev polynomial of degree l( j), and ϕ k( j) (x) = k( j)ϕ(x) the phase increment associated with the frequency k( j)ν.The data needed for specifying the first ten representation terms are given in Table 3.We now discuss how the procedure parameters affect the precision of the resulting representation.For this, we first consider a 150-term representation obtained by resetting M = 1.This resetting introduces only a small discrepancy in ν(t) because |c m /c 0 | < 6 × 10 −4 for m > 1.The precision of this representation is several orders less precise than the previous one.As seen by comparing the two panels of Fig. 3, this is because the whole considered time interval is long, and the small discrepancy in ν(t) can accordingly result in significant phase errors.On the other hand, however, the fact that there is no term with either |k| > 8 or l > 4 in the 37-term representation indicates that  26)-( 28).This representation is obtained with  26)-( 28).The residual of representation is plotted against the number J of representation terms.The other procedure parameters are fixed as {M 1 , M 2 , M 3 , K 1 , K 2 , K 3 , L k } = {0, 0, 9, 10, 10, 10, 9}.
the representation model practically converges with respect to K and L. Our representation procedure also converges rapidly with respect to J.This is illustrated in Fig. 4, for which the procedure parameters excluding J are fixed as {M, K, L} = {9, 10, 9}.

Example of chaotic solutions of Hamiltonian system
We consider the following Hamiltonian system: where A solution of this system has three fundamental frequencies.
The first two are the forcing frequencies ν 1 and ν 2 , which are A84, page 6 of 13 time (Myr) time (Myr) frequency (''/yr) Fig. 9. Variations of major fundamental frequencies of the solar system over the time interval from −35 Myr to 5 Myr with origin at J2000.These frequencies (in arcsec yr −1 ) are computed by applying the NAFF algorithm to the proper modes of the secular solar system associated with major planets.The errors in the resulting frequency samples are estimated as the difference between the frequencies computed respectively from the original ephemeris and its quasi-periodic approximation.Also shown in this figure are the Chebyshev approximations of these varing fundamental frequencies (red curves), which are specified in Table A.2.
non-commensurable.The other, denoted as ν 3 , can be approximated in a similar way as in the previous subsection.
From the analysis of the frequency map (e.g., Laskar 1999), defined as I 0 → ν 3 with θ 0 = 0 and shown in Fig. 5, we know that the phase point lies in a chaotic zone formed by resonance overlap, and the solution z(t) starting from this point is chaotic.
An ephemeris of z(t) is obtained by the symplectic integrator SBABc4 (Laskar & Robutel 2001) at {t i = ih : i = 0, . . ., 4095, h = 0.186058}.We then apply NAFF to 129 evenly spaced time intervals with length d = 512 h and midpoints {τ λ = d 2 + 28(λ − 1)h} 129 λ=1 .This gives the sample set {τ λ , ν 3,λ } 129 λ=1 , partly shown in Fig. 6 together with a Chebyshev approximation.It should be noted that there are two intrinsically different error sources in this way of approximating a changing frequency.One error source is related to the NAFF process that gives the samples of the frequency, while the other is related to the fitting process that leads to a Chebyshev approximation of the frequency.Accordingly, we consider the Chebyshev approximation as sufficiently good if it deviates from the frequency sample set by much less than the sample uncertainties.As shown in Fig. 6, this requirement can be met with a Chebyshev polynomial of degree M 3 = 9.
Test calculations show that our representation procedure cannot lead to an acceptable representation of the solution on the whole sampling time-interval.There are two possible reasons for this.
The most intrinsic reason would be that our solution experiences passages into or out of resonance zones.Such a passage is associated with a shift between circulation and libration of the corresponding resonance angle, and accordingly, with occurrence or disappearance of certain terms.If some of these terms are significant enough in the whole considered time interval (τ 1 , τ 129 ), then there would be no way to obtain any acceptable non-piecewise representation.Therefore, we restrict ourselves to the shorter time interval (τ 1 , τ 50 ).
Another possible reason is that there are one or more significant libration frequencies, which are not taken into consideration when we generate the function basis B. While it is easy to make a necessary extension of B in order to include known libration frequencies (see Sect. 5), it is not that straightforward to identify and sample these frequencies (Laskar 1990).To show the flexibility of our algorithm, we do not search for any libration frequency and make the corresponding extension of B. The flexibility arises because when we choose reasonably high values of our procedure parameters K n , then the resulting set of frequencies would be dense enough over a large frequency interval, in the sense that every important libration frequency is not far from at least one element of the frequency set.
Setting the procedure parameters as 9,10,10,10,9, 100}, we obtain a 100-term representation of our solution.Figure 7 shows that the errors of this representation are typically smaller than 10 −4 I 0 .Figure 8 also illustrates, in the same way as in the previous subsection, the convergence property of the representation procedure.

Application to planetary ephemerides
Numerical integration is now an efficient way of constructing ephemerides of the solar system bodies with high precision.
For practical applications, however, it can be useful to represent these discrete solutions analytically, and thus represent them in a continuous way.These representations can be made in the form of segmented Chebyshev expansions, like those that represent the classical planetary ephemerides as INPOP (e.g., Fienga et al. 2008), or other generally applied approximation models without physical basis.A drawback of these representations is that they require large amounts of data.In order to obtain compact representations, Chapront (1995) used a model of Poisson series with fixed main frequencies that were obtained with the NAFF algorithm.Although this model already involves some long-term or long-period-term effects by allowing Poisson terms and therefore works well with planetary ephemerides of five outer planets over a few hundred years, it does not take the frequency drifts into consideration and cannot be used over millions of years.
The frequency drifts are important over a few tens of million years, as shown by Laskar (1990).The algorithm developed here is therefore expected to be more appropriate in representing ephemerides spanning this long time-interval.It is thus interesting to test whether our algorithm can be used to represent the long-term numerical solution of major solar system bodies over this long timescale (e.g., Laskar et al. 2004).For this, we applied our algorithm to the eccentricity and inclination variables of the Earth, that is, To be more precise, e 3 and i 3 are the eccentricity and inclination, respectively, of the instantaneous orbit of the Earth-Moon barycenter, and 3 and Ω 3 are the longitudes of the perihelion and of the node of the same orbit with respect to the fixed J2000.0 equatorial reference system, respectively.The chaotic behavior of the terrestrial orbit certainly limits the time span over which this orbit can be precisely determined, but z 3 and ζ 3 from La2004 are reliable and precise at least over the time span [−35 Myr, +5 Myr], see Laskar et al. (2011).Therefore, we restricted our representations to this time span.Test calculations showed that our algorithm can lead to compact and precise representations for both degrees of freedom.With different procedure parameters, the resulting representations contain very different terms.This confirms the flexibility of our algorithm.
In order to give representations as physical as possible, we resorted to the knowledge we have for the solution.Following Laskar et al. (2004), we computed the fundamental frequencies of the secular solar system by applying the NAFF algorithm over time intervals of length 20 Myr for the proper modes (z  4 .From −35 Myr to 5 Myr with a step 0.1 Myr, we generated the samples of the fundamental frequencies (g 1 , . . ., g 8 , s 1 , . . ., s 8 ) corresponding to these proper modes.The resulting nominal value of s 5 is about 0.00000015 arcsec yr −1 .The other frequency samples are shown in the panels of Fig. 9, where the errors are estimated as the difference between the values of a frequency computed from the associated ephemeris and its quasi-periodic approximation.We also show in this figure the Chebyshev approximations of these fundamental frequencies.All of these Chebyshev approximations, the coefficients of which are listed in Table A.2, were obtained by truncating those of degree 15.The truncation criterion was roughly that the discrepancy in a fundamental frequency should not induce an error in phase larger than 2π/10 4 over several tens of million years.
There are two important libration frequencies, r 1 = 0.251085 arcsec yr −1 of the resonance argument 2( Ω • 2 ).To include them as additional fundamental frequencies, we expressed the main frequency as The d'Alembert characteristic, 8 i=1 (m i + n i ) = 1, was used to exclude non-physical frequency index vectors.
To show to which degree our algorithm is practically useful, we discuss here the following two 100-term representations ( 100 for short), where with f = (g 1 , . . ., g 8 , s 1 , . . ., s 8 , r 1 , r 2 ) and k = (m 1 , . . ., m 8 , n 1 , . . ., n 8 , k 1 , k 2 ), the phase increment ϕ k (t) = t 0 k, f(τ) dτ is associated with the main frequency k, f .The way the residuals decrease with increasing number of representation terms in the two representation procedures of 100 is shown in Fig. 10.The leading 40 terms of these two representations are given in Tables A.3 and A.4, respectively, where the terms are reordered according to their real amplitudes and data are rounded to a convenient number of digits.Comparing these results with those given in Tables 4 and 5 of Laskar (1990), we find that they are coherent with each other in the sense that all the main frequencies explicitly identified previously can be found in our tables.
The full versions of Tables A.3 and A.4 are available at the CDS.They are plain-text tables providing all terms of 100 in the form of Eq. ( 31).Together with the data presented in Table A.2 for computing f , these two tables can be used to compute the eccentricity and inclination variables from 100 .The errors of 100 as solution representation are shown in Fig. 11.Based on the information in this figure, we expect that our algorithm should be efficient in producing compact and precise representations of long-term ephemerides of major solar system bodies.
To conclude, we explicitly illustrated the advantage of representing ephemerides by taking into consideration the frequency drifts by comparing the representations given by a direct use of NAFF and the modified present algorithms.The quasi-periodic representations ( NAFF for short) of the same z 3 -and ζ 3 -ephemeris given by the standard realization of NAFF, which is a built-in tool of the algebraic system TRIP5 , have 50 and 41 terms, respectively.NAFF does not recover more terms because it encounters a frequency that at a given level of precision is already recovered in a previous step.To more precisely show the advantage of taking frequency drifts into consideration, we produced for z 3 and ζ 3 representations ( ) with the same number of terms as the corresponding NAFF representation.The comparison between and NAFF is shown in Fig. 12.This figure clearly shows the improvements brought by introducing frequency drifts into the representation model.Notes.Also listed are the used constant libration frequencies r 1 and r 2 (arcsec yr −1 ).

Fig. 1 .
Fig. 1.Phase trajectory of the dissipated solution specified by Eqs.(18)-(20).Also shown is the phase orbit of the unperturbed pendulum system, namely (Eq.(18)) with ε = 0, which passes through the initial phase point of the dissipated solution.The dashed line depicts the separatrix of the unperturbed pendulum.

Fig. 3 .Fig. 4 .
Fig. 3. Errors of two different representations of the dissipated solution specified by Eqs.(18)-(20).In the case of the upper panel, our representation procedure is terminated with a 37-term representation, because this representation reaches the precision requirement δ r = 10 −5 .In the case of the lower panel, with slightly larger error in the frequency approximation, the representation is about 4 orders less precise, though the number of representation terms is increased to 150.

A84Fig. 5 .
Fig.5.Fundamental frequency map of the system specified by Eqs.(26) and (27).In the chaotic zone formed by resonance overlap lies the phase point (θ 0 , I 0 ) = (0, 1.535), which will be chosen as the initial phase point of the considered chaotic solution.

Fig. 6 .
Fig.6.Change of fundamental frequency ν 3 (t) of the chaotic solution specified by Eqs.(26)-(28).The samples of this frequency are computed by applying NAFF algorithm over a sliding time-interval, and the errors are estimated as their respective differences from the frequencies of the quasi-periodic approximation of the solution over the same sliding time interval.Also shown is a Chebyshev appoximation of ν 3 (t).

Fig. 10 .
Fig. 10.Convergence property of the two representation procedures leading to the 100-term representations of z 3 (t) and ζ 3 (t).The residual of the representation is plotted against the number J of the representation terms.

Fig. 11 .Fig. 12 .
Fig. 11.Comparison between the eccentricity and inclination ephemerides of La2004 and their respective 100-term representations over the time interval from −35 Myr to 5 Myr with origin at J2000.Both representations are in the form of Eq. (31), with the Chebyshev expansions approximating the fundamental frequencies specified in TableA.2, and main frequency index vectors and complex amplitudes of all 100 terms in the full Tables A.3 and A.4 available at the CDS.

Table 1 .
Summary of procedure parameters.

Table A .
2. Coefficients of Chebyshev expansion m c m T m approximating the major fundamental frequencies (g 1 , . . ., g 8 , s 1 , . . ., s 8 ) (arcsec yr −1 ) of the solar system over the time interval from −35 Myr to 5 Myr with origin at J2000.

Table A .
2. continued.Table A.3.Leading 40 terms of the 100-term representation z 3 (t), as expressed in Eq. (31).Abs(a k ) × 10 6 Arg(a k ) (degree) Notes.The full table is available at the CDS.Table A.4.Leading 40 terms of the 100-term representation ζ 3 (t), as expressed in Eq. (31).Abs(b k ) × 10 6 Arg(b k ) (degree) Notes.The full table is available at the CDS.