Closed-form perturbation theory in the restricted three-body problem without relegation

We propose a closed-form normalization method suitable for the study of the secular dynamics of small bodies in heliocentric orbits perturbed by the tidal potential of a planet with orbit external to the orbit of the small body. The method makes no use of relegation, thus circumventing all convergence issues related to that technique. The method is based on a convenient use of a book-keeping parameter keeping simultaneously track of all the small quantities in the problem. The book-keeping affects both the Lie series and the Poisson structure employed in successive perturbative steps. In particular, it affects the definition of the normal form remainder at every normalization step. We show the results obtained by assuming Jupiter as perturbing planet, and we discuss the validity and limits of the method.


Introduction
The present paper is about the development of a method to compute a secular normal form in the framework of the restricted three-body problem (R3BP). We are interested in the heliocentric dynamics of a massless particle perturbed by the tidal potential of an external planet, i.e. a planet whose orbit is exterior to the particle's one. The objective is to define a transformation leading to a Hamiltonian function suitable to study the particle's secular dynamics, i.e. we search for a normal form not depending on the fast angles characterizing the problem. Using modified Delaunay variables, the latter are the mean longitudes of the particle and of the planet.
The Hamiltonian for the problem of interest is equal to the sum of two components, a leading term Z 0 not depending on the fast angles and a disturbing function R: The leading term is where n P is the planet's mean motion and I P is a dummy action variable canonically conjugated to the planet's mean longitude λ P such that dλ P dt = ∂H ∂ I P .
At the jth iteration, we perform a Lie transformation leading to a new Hamiltonian H ( j) given by where Z i , i = 1, . . . , j − 1 are normal form terms arising at previous normalization steps. We have ) s = R s , ∀s ∈ [1, +∞). Following the above iterative procedure, the Hamiltonian H (r ) after r normalization steps is given by the sum between the normal form Z = r s=0 Z s , and a remainder with a (hopefully) small size. Each generating function χ j , j = 1 . . . r , is determined by solving the so-called homological equation with {Z 0 , χ j } = − n ∂ ∂λ + n P ∂ ∂λ P χ j , n = (G M ) 2 3 .
In our problem, solving (2) can be complex, since the disturbing function R is not directly available as a sum of trigonometric terms over the problem's canonical angles. In fact, R depends on the mean longitudes of the planet and of the particle through geometric angles, i.e. the true anomalies or the eccentric anomalies. This implies that we have to solve Kepler's equation in series form in order to obtain the required trigonometric expansions in the angles λ, λ P . Typically, to overcome this difficulty the original Hamiltonian is approximated by means of a Taylor expansion in some small parameter truncated at an adequate order to make explicit the dependence on the fast angles [(see Tisserand (1889); Brouwer and Clemence (1961); Kaula (1966); Vinh (1970); Brumberg and Fukushima (1994)]. As an example, consider an expansion over the eccentricities of both the particle and the planet. The book-keeping parameter in (1) will depend on e and e P and each term R ( j−1) s will have the form s,k ( , , ; e P , a P ) are the so-called Laplace coefficients (Brouwer and Clemence 1961). The generating function χ j satisfying (2) is χ j = k R ( j−1) j,k nk 1 + n P k 2 sin(k 1 λ + k 2 λ P + k 3 γ + k 4 θ).
Some examples of application of the above method of expansion can be found in Metris et al. (1993), Wnuk (1988) and Lara et al. (2011). However, an important drawback of this technique is that it can be applied only for lowly eccentric orbits. To address this issue, an alternative, introduced in Palacián (1992) and formalized in Deprit et al. (2001) and Segerman and Coffey (2000), is the so-called relegation method; it consists in neglecting one of the two components of the leading term in the homological equation, so that this can be solved in closed form. A multipole expansion of the initial disturbing function R is performed so that each R ( j−1) s has the form R ( j−1) s = k f k ( , , ; e P , a P ) r ( , , λ) cos k 1 u( , , λ) + k 2 f P (λ P ) + k 3 γ + k 4 θ , where r is the heliocentric radius of the particle equal to r = a 1−e cos u . Since the planet's trajectory is external to the particle's trajectory, we have n P < n. Hence, treating n P as a small quantity (of 'book-keeping order 1', see below), instead of the homological equation (2), we work with a homological equation involving only the particle's Keplerian terms in the kernel, namely [see Palacián (2002)]: so that whereR (4), now, can be solved in 'closedform' (see Palacián (2002)), i.e. without expanding the function R ( j−1) j in the eccentricities e, e P . However, it is easy to see that, by applying the Lie transformation, the new Hamiltonian contains terms due to the contribute {n P I P , χ j } of the type R ( j,0) j,1 = δ f k ( , , ) cos(k 1 u + k 2 f P + k 3 γ + k 4 θ), δ = k 2 n P k 1 n , whose size may not be sufficiently small, i.e. comparable to the size of the next term to normalize, namely R ( j) j+1 . As a consequence, the iterative normalization process must be adjusted: additional iterations must be added to handle first the terms coming from {n P I P , χ j }, before dealing with R ( j) j+1 . This leads also to remainder terms depending on increasing powers of δ, whose size decreases, provided that δ < 1. Generating functions χ k j , k = 1, . . . , m, satisfying are iteratively computed to normalize the contributions R ( j,k−1) j,k depending on δ k . After m steps, the final remainder R ( j,m) j,m+1 will depend on δ m+1 and will have a size comparable to R ( j) j+1 . While the relegation technique successfully remedies the issue of the poor convergence of series depending on powers of the orbital eccentricities (see San-Juan et al. (2004); Ceccaroni and Biggs ( 2013); Feng et al. (2015); Palacián et al. (2006) for applications), the practical applicability of the technique is severely limited by the requirement of smallness of the ratio n P /n. To understand this, let us use the example of a two degree-of-freedom harmonic oscillator system in action-angle variables as proposed in Segerman and Coffey (2000): with J 1 , J 2 the actions and φ 1 , φ 2 the angles. The homological equation to solve is where is a formal parameter to assess the size of the terms. The classical solution is However, in the case ω 1 /ω 2 << 1 we can apply the relegation technique: we can neglect ω 1 J 1 in the leading term and determine a generating function χ R satisfying Because of the term ω 1 J 1 in the Hamiltonian, the Lie transformation gives If k 1 ω 1 /k 2 ω 2 < 1 the remainder size is lower than the size of the normalized term. However, it may be higher than the targeted size, so that the relegation process must be further iterated. As discussed in Segerman and Coffey (2000), the iterations produce the generating function It is, now, trivial to see that the generating function χ R , obtained by relegation, corresponds to the series expansion of the usual generating function χ C , obtained without relegation, in powers of the ratio k 1 ω 1 /k 2 ω 2 . However, it is obvious that, even if ω 1 << ω 2 , the method may not converge if the coefficients k 1 , k 2 are such that k 1 ω 1 /k 2 ω 2 ≥ 1. We refer to Sansottera and Ceccaroni (2017) for more details about the convergence of the relegation algorithm. Some methods alternative to relegation have been proposed in the literature to solve the homological equation in closed form. In Mahajan et al. (2018), a technique based on the method of characteristics is developed; its application is shown in Mahajan and Alfriend (2019). In Lara et al. (2013), the homological equation is solved in closed form for orbits with low eccentricity by accepting a remainder of small size depending on e.
All the above techniques were applied, so far, in the so-called satellite problem, i.e. the motion of a test body in the multiple expansion of a planet's gravitational potential (e.g. with the J 2 and C 22 terms). In the present paper, we examine, instead, the applicability of a closed-form normalization method in the framework of the R3BP suitable for orbits with relatively high eccentricities and not using relegation. Our method is similar in spirit to the one introduced in Lara et al. (2013) for satellite motions in the geopotential. In particular, after a multipole expansion of the initial disturbing function, we introduce a book-keeping symbol (with numerical value equal to 1), and write the initial Hamiltonian as where we have The exponent of the book-keeping parameter in each perturbing term R (0) s keeps track of the order of smallness of the term, which, in turn, may depend on one or more of the following three small quantities: e, e P and the ratio between the planet and the Sun's masses. As in Lara et al. (2013), to overcome the difficulty of solving the homological equation in closed form the main idea is to accept a remainder coming from the homological equation itself; at each jth iteration, j = 1 . . . r , we determine a generating function χ ( j) where Z s 0 + j−1 does not depend on λ and λ P . The new Hamiltonian is s contains also the remainder contributions coming from the homological equation. The structure of the paper is as follows. The method will be detailed in Sect. 2. In Sect. 3, we apply the method and give numerical results for the simplest case of the planar circular restricted three-body problem (PCR3BP); an analysis of the results is performed to assess the validity of the method. In the present study, we consider Jupiter as the perturbing planet and a main belt asteroid as the test particle. In Sect. 4, we report the outcomes obtained by applying, instead, the method to some orbits in the more general elliptic R3BP.

Normalization method
In this section, we describe the formal steps required to apply the proposed closed-form normalization method. They include the preparation of the initial Hamiltonian, the choice of the book-keeping scheme, the definitions related to the used Poisson structure as well as the normalization process through the composition of Lie series.

Hamiltonian preparation
Let us consider a heliocentric inertial reference frame with the x axis pointing toward the planet's perihelion and the z axis parallel to the planet's orbital angular momentum. The Hamiltonian of the R3BP is where r is the particle's heliocentric position vector, r = |r|, and p is the conjugated canonical momenta vector, with p = |p|. In Eq. (8), R is the perturbing planet's tidal potential equal to with r P the position vector of the planet, r P = |r P |.

Multipolar expansion
We are interested in analysing the motion of small bodies orbiting the Sun for which we always have r < r P . Then, the function R can be approximated with its truncated multipole expansion: where cos α = r · r P rr P and P j (·) are Legendre polynomials. The time-dependent term 1/r P is omitted in (9) since it does not contribute to the particle's equations of motion.

Extended Hamiltonian
The Hamiltonian (8) can be expressed as a function of orbital elements, using the relations and f is the true anomaly. To avoid trigonometric functions at the denominator in R, it turns convenient to introduce the eccentric anomaly u in place of f through the relations The planet orbit is assumed Keplerian, so that only the true anomaly f P varies in time. The variable f P depends on time through the orbit's mean longitude λ P . However, the Hamiltonian can be formally extended to an autonomous one by adding a term depending on a dummy action I P conjugated to the angle λ P . The extended Hamiltonian is where n P is the planet's mean motion. The dependence of H on the modified Delaunay variables ( , , , λ, γ, θ) is implicit, through the orbital elements, and the dependence on λ P is also implicit, through f P .

Expansion of the semi-major axis
A key element of our proposed method is the following: for algorithmic convenience purposes, it turns out quite useful to have constant frequencies appearing at the kernel of the homological equation to be solved at successive normalization steps. This can be achieved in the following way: recalling that we first choose a reference value a * for the semi-major axis, around which the trajectories to be studied by the normal form oscillate. Setting, now we readily find Then, the Keplerian term in the Hamiltonian can be expanded in powers of the small quantity δ : δ 2 a * 2 + · · · . Introducing the above expansion, the Hamiltonian takes the form (apart from a constant) H = n P I P + n * δ − 3 2 δ 2 a * 2 + · · · + R(a * + 2 n * a * δ + · · · , e, i, ω, , u, f P ; a P , e P ), with The angle u depends on the canonical variables, u = u(δ , , λ, γ ), through Kepler's equation where M = λ + γ is the mean anomaly. We note that this expansion of the Hamiltonian in powers of δ is equivalent to the canonical transformation , , , I P , λ, γ, θ, λ P → δ , , , I P , λ, γ, θ, λ P .
The trigonometric reduction (13) yields a sum of trigonometric monomials cos(k 1 u +k 2 f P + k 3 ω + k 4 ); moreover, after RM-reduction all terms in H appear divided by r except for the terms n P I P and n * δ .

Book keeping
A book-keeping symbol , with numerical value = 1, is used in order to keep track of the relative size of the various terms in the Hamiltonian. There are four different small parameters to consider in the problem: μ, δ and the two eccentricities e and e P . We adopt the following 'book-keeping rules' to assign a unique power of the symbol (reflecting the order of smallness) to each term in the Hamiltonian: • all terms depending on powers of the eccentricities e j e k P , with j, k ∈ Z, are multiplied by the book-keeping factor ( j+k) ; • all terms depending on (1 + η) j and (1 − η) k , with j, k ∈ N, are multiplied by 0 and 2k , respectively; • all terms depending on μ j δ k , with j, k ∈ N, are multiplied by ( j+k)s 0 with s 0 given in (7); • all terms depending on δ k , with k ∈ N, coming from the Keplerian contribute in the Hamiltonian, are multiplied by (k−1)s 0 ; • all terms depending on φ k , with k ∈ N, are multiplied by k .
The quantity φ = u − M is called 'equation of the center'. By Kepler's equation, we have φ = e sin u. After the assignment of the above book-keeping factors, the Hamiltonian is split into two main components, i.e. a leading term Z 0 and the disturbing function R, where To perform the above operation, and in particular to specify the value of the lowest bookkeeping order s 0 in the perturbation, we must have an estimate of the size of e along any individual trajectory: in the numerical examples below we use the initial value e(t 0 ) for this purpose.
In terms of the above book-keeping, the goal of the normalization becomes, now, to define a Lie series transformation leading to a final Hamiltonian normalized up to a pre-selected order s m in the book-keeping parameter . In particular, after s m − s 0 + 1 normalization steps, the Hamiltonian will have the form: where Z is in normal form. All terms with book-keeping order higher than s m are considered negligible in the initial Hamiltonian. Then, the starting Hamiltonian for computing the normal form is set as: The order s m is called the maximum truncation order of the expansion. Let us remark that if we target a remainder with a size of order m P M 2 , we must impose .
For a remainder of order m P M , with > 2, we have, instead, s m = s 0 − 1.

Poisson structure
All along the normalization in closed form, we need to compute Poisson brackets of the form {A 1 , A 2 }, where A 1 and A 2 are functions of (δ , , , I P , λ, γ, θ, λ P ) whose explicit expressions are given in terms of the orbital elements (e, i, ω, , u, f P ) and of the variables r , φ, η: To compute {A 1 , A 2 } we use the formula The partial derivatives in the Poisson bracket (19) are computed using the formulas in "Appendix A", which allow for various simplifications during the normalization process. An important remark regarding these formulas is that the book-keeping parameter appears explicitly in all the partial derivatives which depend on the variables e, e P and δ . This is an essential element of the method: supposing that A 1 =Ã 1 j and A 2 =Ã 2 k , j, k ≥ s 0 , the result of {A 1 , A 2 } will not be exclusively of order j + k in , but will contain also several terms with different powers of . In particular, we have the following Proposition 1 Given two functions A 1 =Ã 1 j and A 2 =Ã 2 k withÃ 1 andÃ 2 in the form (18), the Poisson bracket {A 1 , A 2 } generates terms whose minimum order in is equal to The proof of Proposition 1 is given in "Appendix B".

Homological equation
As mentioned in the introduction, at the j-th iteration of the normalization process, we must determine a generating function χ ( j) s 0 + j−1 satisfying a homological equation of the form We now give the precise form of the homological equation.
By applying formula (19) and using the expressions in "Appendix A", the Poisson bracket We then define χ ( j) s 0 + j−1 by solving the equation: The solution of (22) is found as follows: the function Z s 0 + j−1 contains all the terms of R ( j−1) s 0 + j−1 not depending on λ and λ P . Beside these terms, the function R ( j−1) s 0 + j−1 contains four more different types of terms: Depending on the type of encountered term to be normalized, the generating function χ ( j) s 0 + j−1 must acquire a corresponding term equal to: Then, the outcome of the operation {Z 0 , χ ( j) s 0 + j−1 yields terms in the normal form having the form as follows: • for each normalized term of type 1: f (e, i, η, ω, ), • for each normalized term of type 2: 0, • for each normalized term of type 3:f (e,i,η,ω, ) • for each normalized term of type 4: 0.
We note that the residual of the normalization is equal to zero only for the terms of type 1. Another important remark regards the average value < χ ( j) s 0 + j−1 > of the generating function χ ( j) s 0 + j−1 with respect to the angles λ,λ P . We have that Thus, the average < χ ( j) s 0 + j−1 > is different from zero. This generates no problem for the iterative application of the method. However, it is customary to subtract from χ ( j) s 0 + j−1 > in order that the elements found in the normal form properly correspond to mean elements (see Lara et al. (2014)). We collect in "Appendix C" all the formulas required for the computation of the average < χ ( j) From Proposition 1, we have that the terms generated by the Lie transformation are of order higher than the term normalized at each step when s 0 > 1. In case s 0 = 1, this no longer holds true. We, then, have two distinct algorithms to perform the normalization depending on whether s 0 > 1 or s 0 = 1.

Normalization process for s 0 > 1
The normalization process consists of determining a succession of Lie transformations leading to the targeted normal form. If s m is the targeted order of the final normal form, s m −s 0 +1 steps must be performed. At each step j, the goal is to normalize the Hamiltonian H ( j−1) obtained at the previous step. For this purpose, the homological equation (Eqs. (20) and (22)) is solved to determine the generating function χ ( j) s 0 + j−1 ; the new Hamiltonian is (16)).
The remainder terms R ( j) (s) , with s ≥ s 0 + j, contain three parts: (20); (iii) the terms generated by the the Lie transformation, i.e. coming from Concerning the last part, from Proposition 1 we have Then, the smallest order of the terms coming from the Lie transformation is equal to 2s 0 , for j = 1, or equal to 2s 0 + j − 2, for j > 1. Since s 0 > 1, we have that the remainder is always of order higher than s 0 + j − 1, i.e. the order of the normalized term in the Hamiltonian H ( j−1) . A detailed example of the normalization process for s 0 > 1 is given in "Appendix D". We note that the case s 0 > 1 is rather generic, in the sense that it applies to all trajectories except for those with e = O(m P /M ).

Normalization process for s 0 = 1
Size of the remainder If the particle's orbital eccentricity e is very small (e ∼ O(m P /M )), we obtain from (7) s 0 = 1. In this case, at the generic jth iteration of the normalization algorithm it is easy to see that the operator H ( j) = exp(L χ j )H ( j−1) produces remainder terms of the same book-keeping order as those normalized. Consider the Poisson bracket From Proposition 1, we have However, (ii) The normal form term Z 2 satisfies the relation The proof of Proposition 2 is given in "Appendix E". From it, it follows that only at the second step of the normalization process the Lie transformation will generate terms with the same order as the normalized term. We show now how to deal with this problem by performing just one more additional normalization step.

Adjustment of the normalization process
The normalization process must be modified as follows: • The first step is as in the case s 0 > 1.
• The second step consists of two sub-steps; in the first sub-step, the generating function χ (2) 2 is determined leading to the new Hamiltonian In the second sub-step, the generating function χ (2,bis) 2 is computed as described above and the new Hamiltonian is H (2,bis) = exp(L χ (2,bis) 2 )H (2) .
• In the third step, the Hamiltonian H (2,bis) is normalized up to the third order in ; the Lie transformation leads to the new Hamiltonian 3 )H (2,bis) .
• Successive iterations beyond the order 3 are performed as in the case s 0 > 1.

Numerical application in the PCR3BP
We applied the method described in Sect. 2 in the case of the PCR3BP considering Jupiter as the perturbing planet. The orbital planes of the body and planet coincide and the planet orbit is assumed circular (e P = 0). This implies that f P = λ P , i = 0 and = 0 so that the Hamiltonian does not depend on the Delaunay variables and θ . We perform two tests to assess the applicability and precision of the method. As a first test, we estimate the size of reminder H R = +∞ s≥s m +1 R (s m −s 0 +1) s of the normal form and compare it to the size of the initial disturbing function R = H (0) − Z 0 (see (16)). We perform a multipolar expansion of degree 10 and a normalization up to a certain order s m in book-keeping set as s m = min(2s 0 − 1, s 0 + 10) with s 0 given in (7). Such a choice is empirically found to yield a good compromise between computational load and requirements for precision.
To obtain estimates of the remainder size, we consider a truncation of the remainder up to terms of book-keeping order s m + 3: the size of H R can be estimated through the norm The norm of R was computed with the same definition. Figure 1 shows log 10 ||H R ||/||R|| in color scale in a grid of values for the initial semi-major axis a(0) = a * and eccentricity . The quantity log 10 ||H R ||/||R|| gives an estimate of the relative size of the remainder with respect to the initial perturbation, which estimates, in turn, the relative error of the semi-analytically computed trajectory with respect to the true trajectory. Denoting e 0 = e(0) and a 0 = a(0), in Fig. 1 the red line corresponds to the set of points (a 0 , e 0 ) such that a 0 (1 + e 0 ) = a P , i.e. the radius at the apocenter coincides with the Sun-Jupiter distance (for the PR3BP r P = a P ). Since our method is applicable to particles with trajectories lying entirely inside the trajectory of the planet, the red line represents an upper boundary of the region in the (a 0 , e 0 ) plane in which the method can be applied. The black line represents the upper boundary of the values (a 0 , e 0 ) for which we have Hill-stable orbits. An orbit is defined as Hill-stable when its Jacobi constant C is larger than the Jacobi constant C L 1 at the Lagrangian point L 1 . The boundary drawn was determined as described in Ramos et al. (2015).
From Fig. 1, we can observe that the relative error is lower than 10 −2 for each value of the eccentricity up to an initial semi-major axis lower than ∼ 0.3a P ; up to ∼ 0.45a P it raises above 10 −1 only for the higher values of e(0). On the other hand, getting closer to the Hillunstable region, the error becomes higher and keeps having acceptable values (of few percent) only in regions with low eccentricity. In the Hill-unstable region, the error is everywhere high. In the figure, we can notice also several vertical strips along which the error is always higher than in their neighbourhood. These strips correspond to mean motion resonances, in which the method fails due to small divisors appearing along the normalization process; note that in the Hill-stable region, the error is high at those domains where the concentration of these strips becomes more conspicuous.
As a second test, we compare the semi-analytical computation of the evolution of the orbital semi-major axis and eccentricity using the normal form with the results obtained through the numerical propagation of the particle's trajectory.
After s m − s 0 + 1 steps, the normalization process transforms the original canonical variables, (δ (0) , (0) , I with analogous formulas holding for all the other variables. Since the initial conditions of any trajectory are given in the original variables, to compute the initial conditions in the new variables, the inverse transformation must be used; for example, we have Having computed their evolution, we can obtain also the evolution of the semi-major axis and the eccentricity as n * a * 2 + δ (0) . The disturbing function was computed as in the previous test; to save computational time, only 4 steps were carried out in the normalization process. The numerical propagation was performed with MATLAB using the function ode45. Figure 2 shows the maximum relative errors obtained in function of the initial values of the semi-major axis and eccentricity (a 0 , e 0 ). To interpret the results, we also computed the stability map shown in Fig. 3. It was obtained by computing the fast Lyapunov indicators (FLI) (Froeschlé et al. 1997) for orbits with the same initial conditions presented above, using a propagation time equal to 20 orbital periods. In Figs. 2a, b and 3, the red and black lines are the same as described above. We observe that in the domain left to the black line, the error is generally small, except along the vertical strips corresponding to mean motion resonances and their neighbourhood, similar as in Fig. 1. The stability map (Fig. 3) confirms these features due to mean motion resonances. Figure 2 shows that the error of the error increases, in general, as a and e increase. For higher values of a, the error is mostly dominated by the truncation level of the multipolar expansion. Fixing a, the error as e increases is regulated, instead, by the choice of maximum normalization order. From the tests performed, we can conclude that with the adopted truncation and normalization orders the method produces accurate results up to an initial semi-major axis a 0 ≤ 0.6a P ; for higher values of the initial a, either we accept an higher error or we must increase the order of the multipolar expansion which implies a substantial increase in computational time. On the other hand, the method gives accurate results for still high values of e 0 , up to almost 0.7: again more accurate results can be obtained for still higher values of e 0 , by performing a large number of normalization steps at higher computational cost.

Note that both a (s
Let us finally remark that, while in the case of the PCR3BP the normalization method described above can be used directly for the computation of the proper semi-major axis and proper eccentricity, the method lends itself conveniently as the starting point for the computation of proper elements also in more complex cases, e.g. when the perturbation effects of external planets are considered or in the more general spatial elliptic R3BP.

Numerical application in the spatial elliptic R3BP
We reproduced the second test described in Sect. 3 in the more general of the elliptic R3BP selecting few initial conditions for the particle's trajectory. In particular, we considered the following cases In all cases also impose δ (0) = 0, ω(0) = 90 • , M(0) = 90 • , i(0) = 20 • , (0) = 0 • , λ P (0) = 0 • . The number of terms in the initial disturbing function is higher than in the case of the PCR3BP; to keep the number of operations relatively low we perform a multipolar expansion of order 5. We carried out 4 steps during the normalization process as in the numerical examples for the PC3BP. Figure 4 shows the outcomes obtained for the semimajor. The method works well in the first three cases yielding a maximum relative error ∼ 10 −4.3 for case 1 and ∼ 10 −3.7 in the other two cases. Instead, in case 4 the method does not work properly; indeed the maximum relative error we get is 10 −1 . Similarly, for the eccentricity (Fig. 5) the maximum relative error is ∼ 10 −3.9 for case 1, ∼ 10 −3.7 for case 2, ∼ 10 −2.6 for case 3 and ∼ 10 −0.12 for case 4.
These results generally confirm the conclusions obtained for the PCR3BP. The method is able to produce accurate outcomes also for high eccentricity if the initial semi-major axis a 0 is sufficiently lower than a P . For high values of a 0 the relative error depends also on the maximum order of the multipolar expansion of the original Hamiltonian. Increasing the multiple order produces a lower error, but also causes a significant increase of the computational time. Moreover, for a fixed normalization order the error increases with the orbital eccentricity: to rectify this trend it is necessary to perform a larger number of steps during the normalization process as the initial value of the eccentricity grows. For example, repeating the test for case 3, but performing 7 steps of the normalization process, the resulting maximum relative error decreases: it reduces to ∼ 10 −4.3 for the semi-major axis and to ∼ 10 −2.8 for the eccentricity (see Fig. 6).

Conclusions
In the present paper, we discussed a closed-form normalization method which can be applied to the study of the secular dynamics of the massless particles whose heliocentric orbit is affected by a planet external to the orbit. The main novelties of the hereby presented method are: (i) the lack of any relegation, and (ii) the use of a book-keeping scheme applied not only to the Hamiltonian, but also to the Poisson structure used at all perturbative steps. In particular: • Section 2 explains in detail all the algorithmic steps for the implementation of the method.
These steps include: i) the preparation of the initial Hamiltonian (Sect. 2.1); ii) the implementation of the book-keeping scheme (Sects. 2.2, 2.3); and iii) the strategy to solve the homological equation in closed-form (Sect. 2.4). The consistency of the method is rigorously demonstrated with the help of Propositions 1 and 2. We also discuss variants of the method which have to be used in the case of orbits with small eccentricities (e ∼ O(m P /M )). • We have illustrated a numerical application of the method in the framework of the planar restricted three-body problem. Due to the relatively small number of series terms in that case, all expansions can be pushed to relatively high orders in our chosen book-keeping parameter, which practically corresponds to computing a O (m P /M ) 2 secular Hamiltonian, along with the complete formulas for the corresponding normalizing canonical transformation in closed-form. A map of the precision level of the method is given in Fig. 1 (relative error of the remainder) and through numerical estimates in Fig. 2a, b. Overall, in the Sun-Jupiter case, we obtained relative errors in the reproduction of the time series of the osculating elements of numerically propagated trajectories of the order of 10 −6 for NEA objects and 10 −3 for a main-belt asteroid with eccentricities up to 0.7. The exact structure of the map of the precision in the (a, e) plane depends mostly on the existence of mean motion resonances; this is expected since very small divisors are introduced near any such resonance. In fact, the treatment of mean motion resonances requires a different method, albeit still possible to implement in closed-form. This subject is proposed for further studies. • We have also illustrated a numerical application of the method in the framework of the more general spatial elliptic restricted three-body problem. In that case, the number of series terms to be normalized quickly grows and poses the most severe limits to the implementation of the method in the computer. Preliminary results suggest an overall error of the order of ∼ 10 −3 , even for highly eccentric orbits (e ∼ 0.7), but in areas with a < 0.6a P in the Sun-Jupiter case.
As a general conclusion, the method allows to produce quite accurate results up to high values of the eccentricity, but with a drawback relative to the computational cost; indeed, the number of terms to be normalized quickly reaches the limits of the personal computer for trajectories with eccentricities e > 0.7 and semi-major axis a < 0.6a P . This is normal, since every term in closed-form encodes the information of many terms in the associated expansion in powers of the orbital eccentricity. On the other hand, it is well known that closed-form methods can benefit from algebraic simplifications (see Lara (2021)) that can be suitably implemented in special manipulators. In the present paper, we do not discuss any potential optimization of the algebraic structure of the method, since our main focus is on demonstrating the possibility to avoid relegation, which would otherwise be the major limitation to any closed-form approach in the restricted three-body problem. We, thus, leave open the possibility of such optimization. We finally remark that a method analogous to the one presented here can be implemented also to deal with the perturbation on a small body by an internal planet, as discussed in Rossi and Efthymiopoulos (2022).
Acknowledgements I.C. has been supported by the H2020 MSCA ETN Stardust-Reloaded, Grant Agreement n. 813644. C.E. acknowledges the support of MIUR-PRIN 20178CJA2B 'New frontiers of Celestial Mechanics: theory and applications'.
Funding Open access funding provided by Università di Pisa within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Poisson structure: formulas
We report the formulas used to compute the Poisson bracket (21). For any function for any A = A 1,2 of the form (18): The partial derivatives in the formulas above are: Whenever needed, the higher order terms O( s 0 δ ) in the above formulas are easy to obtain by computer algebra.
If implemented in the manipulator in the exact form given above, the previous expressions allow for various simplifications during the normalization process.

B Proof of Proposition 1
We analyse the powers in of the various terms produced in the Poisson bracket {A 1 , A 2 }. The partial derivatives with respect to θ and do not introduce any order variation, thus Similarly, In fact, the partial derivatives with respect to λ P generate terms of orders j + k or higher. The same holds true for The only part of {A 1 , A 2 } which produces terms of order lower than j + k is We have Then, the quantity generates terms of order j + k − 2, while the quantity generates terms of order j + k − 1. If both j, k > s 0 , the minimum book-keeping order encountered in the above derivatives is j +k −2 and the proof of Proposition 1 is completed. We will show, now, that the minimum book-keeping order becomes j + k − 1 if either j = s 0 , k > s 0 or k = s 0 , j > s 0 . To this end, the following lemma must be used:

Lemma 1 Consider the angular variables
where L E is the eccentric longitude, L T ,P is the planet's true longitude, andω andω P are the longitudes of the pericenter of the particle. The Hamiltonian function contains terms of the form f m (δ , i, η, r ; a * , a P , e P )e p P e q cos(m 1 L E + m 2 L T ,P + m 3ω + m 4 + m 5ωP + m 6 P ) fulfilling the following D'Alembert rules: • m 1 + m 2 + m 3 + m 4 + m 5 + m 6 = 0 ( 2 4 ) • q − |m 3 | is nonnegative and positive (25) • p − |m 5 | is nonnegative and positive.
We recall that in the selected reference frame P = ω P = 0. If, now, j = s 0 , A 1 does not depend explicitly on e and φ (as a consequence of the adopted book-keeping rules). The only small parameter on which it can depend will be either μ or δ 2 . Thus Moreover, from Lemma 1 it follows that A 1 does not containω. This implies that either it does not depend on u and ω or it is of the form A 1 = j f k (η, i, r ; e P , a P , a * ) cos(k 1 u + k 2 f P + k 3 ω + k 4 ), k 3 = k 1 .
the Poisson bracket {A 1 , A 2 } produces terms of order equal or larger than j + k − 1 (since the quantity in (27) is zero).

Remark 4
To automatically obtain all the terms with the correct book-keeping order, all terms generated by expressions of the form ∂ A 1 ∂γ are automatically adjusted to appear with the same exponent of r in the denominator. Then, in the numerator of all the resulting terms the variable r is substituted with its expansion r = a * (1 − e cos u) + O( s 0 δ ).
Remark 5 If we target a normal form of order s m , as defined in (17), all the contributions O( s 0 δ ) in the partial derivatives used to compute the Poisson bracket can be neglected.

Remark 6
Once having computed the Poisson bracket by applying formula (19), we substitute φ with e sin u in all produced terms depending on the equation of the center.
Proof of Lemma 1 Consider the expression of the initial Hamiltonian (Eq. (13)) before the introduction of the book-keeping parameter. The terms obtained after the substitution a = a * + δa and the expansion in δ (see Sect. 2.1) come from two parts: i) the initial Keplerian term, and ii) R determined through the multipole expansion of the planet's tidal potential (see (9)). Moreover, the Hamiltonian contains the term n P I P .
From Eq. (9), R contains terms of the type We have where By performing the transformation (23), we obtain T 1 = − 1 2 ae cos i + 1 cos(−ω + L T ,P ) + 1 2 ae cos i − 1 cos(ω + L T ,P − 2 ), Then, cos α fulfils the D'Alembert rules (24), (25), (26). Moreover, since r = a(1 − e cos u) = a(1 − e cos(L E −ω)), 1 r P = 1 + e P cos f P a P η 2 P = 1 + e P cos(L T ,P −ω P ) a P η 2 P , it follows that r and 1/r P also fulfil the D'Alembert rules. The product between terms fulfilling the D'Alembert rules fulfils them as well. It follows that the terms coming from R in H fulfil the D'Alembert rules. Now, all the terms of R are multiplied by Q defined in (14) (RM-reduction). Q results from the expansion in δ of Thus, Q fulfils the D'Alembert rules, implying that the product RQ fulfils the D'Alembert rules as well. Finally, the Lie transformation preserves the d'Alembert rules. Thus, all the Hamiltonians computed throughout the normalization process fulfil the D'Alembert rules.

C Computation of the average of the disturbing function or the generating functions
We report some useful formulas to apply for the computation of the average of the disturbing function and of any generating function with respect to λ. Since dλ = dM, the average of any trigonometric quantity over λ coincides with the average over the mean anomaly M.
From Kelly (1989) where k ∈ Z and ν is a generic angle. From Kelly (1989) and Kozai (1962), we have the following formulas applicable whenever f is used in place of u: D Example of the normalization algorithm for s 0 > 1 We give below a detailed example of the proposed normalization algorithm in the generic case s 0 > 1 (i.e. e >> m P /M ). We consider a toy Hamiltonian in which the initial disturbing function R is given by the quadrupolar expansion (see (9)): with Note that in R s 0 +1 has only terms of type 2.

The new Hamiltonian is
s contains the following contributions : (i) R (s) ; (ii) the remainder of the homological equation; (iii) the terms generated by the Lie transformation. From the previous considerations, we have that these last ones have book-keeping order equal to or larger than 2s 0 .
At the second step, the procedure is repeated with the goal of determining the generating function χ (2) s 0 +1 to normalize H (1) up to the order s 0 + 1 in . The term to normalize is (34) andR (1) s 0 +1 equal to the order s 0 + 1 (in ) term of the remainder computed at the previous step. The generating function χ (2) s 0 +1 then is computed as 2n P − n * + eη + (8 + η + )C 4 sin(2 f P + u + 2g − 2h) 2n P + n * + 2 n P n P − n * e P η 2 while we also have Subsequent steps can be computed using analogous formulas.

E Proof of Proposition 2
(i) Consider the jth normalization step in the case s 0 = 1. From Remark 2 it follows that the Poisson bracket {Z s 0 , χ ( j) j } yields terms of order s 0 + j − 1 = j through formula (28), by taking A 1 = Z 1 and A 2 = χ ( j) j . However, since Z 1 is a normal form term, it does not depend on u and r . Moreover, since it has book-keeping order equal to s 0 = 1, it does not depend explicitly on e (as a consequence of the adopted book-keeping rules). By Lemma 1, we conclude that it does not depend on ω. It follows that the term of book-keeping order j coming from {Z s 0 , χ ( j) j } is equal to zero. (ii,iii) To show the second and third points of Proposition 2, we use the following lemma:

Lemma 2 All terms in R
(1) 2 are of the following two types: • terms of type A, not depending on the eccentricity e and φ; • terms of type B, linearly depending on the eccentricity e (and not depending on φ).
All terms of type A are of the form (29). All terms of type B depend on the eccentric anomaly u, the true anomaly f P , or both.
Since at the second step R (1) 2 is the term to be normalized, it follows that all terms in the generating function χ (2) 2 also satisfy Lemma 2. By Lemma 2, the normal form term Z 2 obtained by normalizing R (1) 2 does not depend explicitly on e and φ. Moreover, by definition, any normal form terms cannot depend on u.
To demonstrate point (iii) of Proposition 2, we finally need to prove that {Z s , χ

All terms in R
(1) 2 containing one of the first five factors are of type A; all terms in R (1) 2 containing the factors 6 and 7 are of type B. Considering that φ = e sin u, by substitution we have that also all terms containing factors 8 and 9 are of type B. From Lemma 1, it follows that the terms of type A are of the form (29).
To show that all the terms of type B depend on u, f P , or both, we need to examine the following parts of R (1) 2 : (i) R (0) 2 , stemming from terms of book-keeping order 2 in the original Hamiltonian, (ii) the remainder terms produced by the homological equation and the Lie transformation at the first step of the normalization process.
(1) Analysis of R (0) 2 Let us consider, first, the initial Hamiltonian (13). In view of the expressions for cos α in (31), r in (11) and r P in (10), we obtain that R (see (9) with D k , C k, j ∈ Q, k, j ∈ Z + , k ≥ 3, j ≥ 1. After expanding the semi-major axis (Eq. (12)), the terms of book-keeping order 1 of R have one of the forms Then, considering also the terms from the Keplerian part and performing the product by Q, R  We conclude that R 1 contains only terms of type 1 or type 2 (see Sect. 2.4) which do not depend on the eccentricity and are consequently of the form (29). Moreover, all the terms of R (0) 2 which depend linearly on the eccentricity (Eq. (36)) also necessarily depend on u, f P , or both, as we readily to infer by considering the expressions of T 0 and T 1 in (32).

andR
(1) LT ,2 the parts of the remainder, respectively, coming from the homological equation and the Lie transformation at the first step of the normalization process. It has been already shown that the terms of R (0) 2 linearly depending on e depend also on u, f P , or both. Furthermore, the remainder term R (0) 1 (normalized at the first step) contains only terms of type 1 and 2 of the form a * r f (i, η, ), a * rf k (i, η) cos(k 1 u + k 2 f P + k 3 ω + k 4 ), k 1 = k 3 with either k 1 ≥ 1 or k 2 ≥ 1. It follows that the residual of the homological equation of the first step is produced only by terms of type 2. To normalize these last terms, χ 1 has to acquire terms of the form k 1 n * + k 2 n Pf k (i, η) sin(k 1 u + k 2 f P + k 3 ω + k 4 ).
Applying, now, the homological equation (20) for the function χ (1) 1 , we obtain that the residual of the homological equation yields remainder terms of book-keeping order 2, which are of the form