Relegation-free closed-form perturbation theory and the domain of secular motions in the restricted three-body problem

We propose a closed-form (i.e., without expansion in the orbital eccentricities) scheme for computations in perturbation theory in the restricted three-body problem (R3BP) when the massless particle is in an orbit exterior to the one of the primary perturber. Starting with a multipole expansion of the barycentric (Jacobi-reduced) Hamiltonian, we carry out a sequence of normalizations in Delaunay variables by Lie series, leading to a secular Hamiltonian model without use of relegation. To this end, we introduce a book-keeping analogous to the one proposed in Cavallari and Efthymiopoulos (Celest Mech Dyn Astron 134(2):1–36, 2022) for test particle orbits interior to the one of the primary perturber, but here adapted, instead, to the case of exterior orbits. We give numerical examples of the performance of the method in both the planar circular and the spatial elliptic restricted three-body problem, for parameters pertinent to the Sun-Jupiter system. In particular, we demonstrate the method’s accuracy in terms of reproducibility of the orbital elements’ variations far from mean-motion resonances. As a basic outcome of the method, we show how, using as criterion the size of the series’ remainder, we reach to obtain an accurate semi-analytical estimate of the boundary (in the space of orbital elements) where the secular Hamiltonian model arrived at after eliminating the particle’s fast degree of freedom provides a valid approximation of the true dynamics.


Introduction
As opposed to the usual (Laplace-Lagrange) theory, closed-form perturbation theory [11] provides a framework for series calculations in perturbed Keplerian problems without expansions in powers of the bodies' orbital eccentricities. This is mainly motivated by the necessity to construct secular models for sufficiently eccentric orbits, like those of many asteroids, in our solar system, or the planets in extrasolar planetary systems.
The efficiency of the usual series methods of expansion in the orbital eccentricities is limited by the fact that the inversion of Kepler's equation in powers of the eccentricity converges only up to the so-called Laplace limit e L ≈ 0.66274 [6]. Generally, such convergence slows down way before this value (around e ∼ 0.3 − 0.4 in many applications). In order to address this issue, closed form perturbation theory aims at solving in 'closed-form' the homological equation by which the Lie generating function is computed at every perturbative step (see for example [3,5]). The process is far from being priceless: a major obstruction appears when the kernel of the homological equations contains addenda beyond the Keplerian terms. The most common such addendum ( [11]) is the centrifugal term −νH, where ν is the angular frequency in a frame co-rotating with the primary perturber, and H is the Delaunay action equal to the particle's angular momentum in the direction of the axis of rotation. In the case of a planet's orbiter, ν is equal to the planet's rotation frequency, and the problem appears for all non-axisymmetric terms (tesseral harmonics) of the planet's multipole potential. In the R3BP, instead, ν represents the mean motion of the primary perturber (e.g. Jupiter in the Sun-Jupiter system), while the problem appears in a similar way after introducing a multipole expansion of the disturbing function in the particle's Hamiltonian.
An algorithm to overcome the above issue, called the relegation algorithm, has been proposed in works by Deprit, Palaciań and collaborators [2,4,7,13,15]. Briefly, given a quasi-integrable Hamiltonian H = H 0 + εH 1 , where ε is a small parameter, suppose that H 0 = H 0 + H 0 , where, in a domain in phase space we have that H 0 yields the dominant contribution to the Hamiltonian flow of H 0 versus the H 0 term. In usual perturbation theory, we seek to partly normalize the perturbation H 1 via a sequence of canonical transformations defined by generating functions χ (r) , r = 1, 2, . . . satisfying a homological equation of the form {H 0 , χ (r) } + h identifying H 0 with the Keplerian term (when ν is small) leads to a homological equation that can be solved in closed form (we set, instead, H 0 = −νH when ν is large). However, all Poisson brackets of χ (r) with the part H 0 left out of the kernel lead to terms which need to be 'relegated', i.e., pushed to normalization in subsequent steps. For reasons explained in detail in [14], only a finite number or relegation steps can be performed before reaching a point beyond which the scheme generates divergent sequences of terms (see also [13]). This implies that the process necessarily stops after some steps, leading to a finite, albeit possibly quite small remainder.
Relegation is a technique particularly suitable to the limiting situation of a strongly hierarchical problem, when the integrable part H 0 depends on a frequency vector involving n frequencies ω = (ω 1 , . . . , ω n ) out of which one, say ω i for some i with 1 ≤ i ≤ n is significantly larger in absolute value than the rest. In particular, the harmonics cos(k · ϕ) in the Hamiltonian whose normalization can be 'relegated' should satisfy |k i ω i | |k j ω j |, j = 1, . . . , n, j = i, for every integer k i , k j ∈ Z \ {0} (assuming also the non-resonant condition k · ω = 0, k = (k 1 , . . . , k n )). For example, as explained in [14] in the simple case with n = 2 and ω 2 ω 1 , the generating function χ (N ) produced after N relegation steps contains terms with coefficients growing as a geometric sequence with ratio k 1 ω 1 /k 2 ω 2 . Thus, relagation is limited to those terms for which the above ratio is smaller than unity. This includes most harmonics of low Fourier order in the Hamiltonian perturbation when ω 2 ω 1 , but only few when the two frequencies become comparable in size. Hence, by construction, relegation has limited applicability in this latter, non-hierarchical, case.
Variants of the relegation technique have been discussed in literature to address perturbed Keplerian problems in which the gravitational potential is due to an extended body expanded in spherical harmonics (e.g. [7,10]). To address the non-hierarchical case, a techique similar to the one of the present paper is discussed in [7], referring to the averaging of the tesseral harmonics in the case of the Earth's artificial satellites. In the case of the R3BP, instead, Cavallari and Efthymiopoulos [1] discuss a relegation-free algorithm for the elimination of short-period terms in the particle's Hamiltonian, when the orbit of the particle (e.g. an asteroid) is totally interior to the orbit of the primary perturber (e.g. Jupiter). We are aware of no relegation-free algorithm proposed in literature which addresses, instead, the case when the particle's orbit is exterior to the orbit of the primary perturber. Providing such an algorithm, discussing some of its important differences with past-proposed algorithms, as well as checking its limits of applicability, constitutes the primary goal of our present paper.
The R3BP is defined by the motion of a body P of negligible mass in the gravitational field of two massive bodies P 0 (the primary or central body) and P 1 (the secondary or primary perturber), which perform a motion r 1 (t) either elliptic in the more general version (ER3BP) or circular (CR3BP). The starting point for our analysis in the sequel is the Hamiltonian of the model, obtained after reduction via Jacobi coordinates (R, P ). 1 . Expressing time through the secondary's mean anomaly M 1 = n 1 t, where n 1 is the mean motion of the secondary, and canonically conjugating M 1 with a dummy action variable J 1 allows to express the Hamiltonian as where G is the gravitational constant and µ = m 1 m 0 + m 1 ∈ (0, 1/2] is the mass parameter; is the elliptic revolution of P 0 − P 1 around their barycenter with eccentricity e 1 and semimajor axis a 1 , in which the dependence of the system's eccentric anomaly E 1 ∈ T = R/(2πZ) on the mean anomaly M 1 ∈ T is given through Kepler's equation according to standard twobody problem setting; (R = (X, Y, Z), P = (P X , P Y , P Z )) ∈ T * (R 3 \ {−µr 1 , (1 − µ)r 1 }) is the position-momentum couple of P and the phase space is endowed with standard symplectic form We make use then of Delaunay elements ( , g, h, L, G, H), defined by where a, e, i, M, Ω, ω stand for the semi-major axis, the eccentricity, the inclination, the mean anomaly, the longitude of the ascending node, the argument of pericenter of the particle.
A key ingredient of the method proposed below is the following: similarly as in [1], we introduce a book-keeping symbol σ with numerical value equal to 1, whose role is to organize the perturbative scheme so as to successively normalize terms of similar order of smallness, treating together all small quantities of the problem, i.e., -the eccentricities e, e 1 (when e 1 = 0), -the mass ratio µ, -the semi-major axis fluctuation δL around the mean L * for a particular particle trajectory. 1 In the R3BP problem the Jacobi transformation is implemented when R > r1 .
The book-keeping symbol acts by assigning powers σ 1 and σ ν 1 , σ ν , σ ν respectively, for non-zero natural numbers ν, ν 1 defined below, to all the terms in the original Hamiltonian as well as in the Hamiltonian produced after every normalization step. Given this baseline, we arrive (in Section 2) to the following result: we demonstrate that, for k µ , k mp ∈ N \ {0} with k µ > 1, the combination of expansions of (1) up to µ kµ and ( r 1 / R ) kmp is canonically conjugate by ν(k µ − 1) near-identity transformations to a secular model, obtained as a normal form with respect to the fast angles , M 1 with The dependencies f = f ( , δL, G) for the true anomaly, e = e(δL, G) and i = i(G, H) are implied in all the above expressions; c l,p , d νkµ,s are real coefficients. A crucial point is the way by which the positive integers ν = ν(e * , µ) ≥ 1, ν 1 = ν 1 (e * , e 1 ) ≥ 1 are chosen. As detailed below, these integers, which regulate the book-keeping scheme, are suitably tuned on the basis of a selected reference value e * ∈ (0, 1): ν = log 10 µ log 10 e * , ν 1 = log 10 e 1 log 10 e * , where · is the ceiling function. The normalizing scheme leading to (4) is local: knowing that the semi-major axis is preserved under the flow of the (secular) normal form, we introduce the splitting L = L * + δL, where L * = √ Gm 0 a * δL, n * = √ Gm 0 a −3/2 * is a targeted reference value for the semi-major axis a * , and expand the Hamiltonian in powers of δL, rendering δL the new action variable canonically conjugated to the particle's mean anomaly.
Given the above, the normalization algorithm provides a sequence of Lie generating functions χ (j) ν+j−1 = O(σ ν+j−1 ), j = 1, . . . , ν(k µ −1), which yields the Lie canonical transformation allowing to recursively normalize all terms depending on the angles f and E 1 in the Hamiltonian. The normalizing trasformations are possible to define for values of the frequencies n * (mean motion of the particle at the semi-major axis a * ) and n 1 far from mean-motion resonances (see Remark 3). Furthermore, the generating functions are computed as solutions of a homological equation of the form where Z 0 = n * δL + n 1 J 1 and R (j−1) ν+j−1,ν+j−1 ∼ σ ν+j−1 collects the trigonometric monomials of O(σ ν+j−1 ) depending on at least one of the two anomalies. The key to obtaining a closedform solution for (8) is, precisely, the appropriate choice of a O(σ ν+j−1 ) remainder left in the second hand of the equation. In words, we do not seek for an exact cancellation of the terms R (j−1) ν+j−1,ν+j−1 , but only for an approximate cancellation, leading to a remainder, which, however, is of higher order in book-keeping, and, hence, possible to reduce at subsequent steps.
As discussed in Section 3, a relevant outcome of the analysis of the behavior of the remainder obtained by the above method stems from an estimation of the optimal number of normalization steps j opt , where the remainder becomes of order ν + j opt − 1 in the book-keeping parameter, with j opt ≤ ν(k µ − 1). The value of j opt is defined as the one where the error bound E (j) (a * , e * ) = ν+j≤l≤νkµ,s |d l,s as in (6) after j normalization steps. As typical in perturbation theory, the value of j opt depends on the chosen reference values (a * , e * ). With the present method one can then obtain a map of the size of the optimal remainder as a function of (a * , e * ) in the semi-plane a > a 1 . Using this information, we compute the limiting locus uniting all points in (a * , e * ) such that the normal form computation yields no improvement with increasing number of normalization steps, i.e., where j opt = 1. Comparing with numerical stability maps obtained with the Fast Lyapunov Indicator (FLI) [9], one sees that, the limiting locus found semi-analytically essentially coincides with the numerical (FLI map) limit where no harmonic in the Hamiltonian associated with one of the exterior mean-motion resonances affects the dynamics. As a consequence, all motions in the sub-domain of the plane (a * , e * ) below the limiting locus are stable in the secular sense, i.e., protected against instabilities caused by short-period resonant effects. For this reason, we identify this locus as the border of the domain of secular motions, and substantiate the fact that its semi-analytical computation (through the normal forms) yields results in precise agreement with those found by the heuristic definition of the same border via the fully numerical (FLI) computation of stability maps.
The paper is structured as follows. Section 2 presents step-by-step the algorithm that gives rise to (5) and (6), supplemented with the formulas for the Poisson algebra in Keplerian elements used in all closed-form computations. Section 3 is devoted to a numerical investigation of the method's accuracy for an asteroid in the Sun-Jupiter system, first in the spatial ER3BP, and then in the planar CR3BP; in the latter case, the computations are short enough to allow for a specification of the optimal normalization order in a grid of values in the (a * , e * ) plane, leading to the semi-analytical determination of the border of the domain of secular motions. Section 4 summarizes the basic conclusions of the present study and gives some relevant comments for future work.
2 The closed-form method for the outermost R3BP

Multipole expansion of the perturbation
Referring to section 1, let H be given in barycentric Cartesian coordinates as in (1): Figure 1: Representation of the R3BP in the barycentric frame (or equivalently in Jacobi variables) with R > r 1 .
indicates the generalized binomial coefficient (equal to 1 for l = 0).

Remark 1.
For l = 1 in Eq.(10) the coefficients of the dipole term (r 1 · R)/ R 3 in the two sums in the r.h.s. of the equation cancel each other exactly. Thus, no dipole term appears in the disturbing function. This is a consequence of the choice of Jacobi coordinates.

Canonical form of the Hamiltonian
Performing an extra series expansion in powers of µ < 1 yields the standard nearly-integrable form where the Keplerian part reads and the disturbing function becomes We now move to Delaunay action-angle variables (3) by replacing into (11) the relationships as well as (2) for the vector r 1 . We get Remark 2. Only the square of the norm r 1 2 = r 1 · r 1 is required in Eq. (13), while the norm R appears only in the denominator of the above equation, in powers equal to or higher than quadratic. Then equations (15) and (2), respectively dependent on f and E 1 , lead to a representation of the disturbing function as a sum of trigonometric polynomials depending on harmonics of the form cos(s 1 f + s 2 g + s 3 h + s 4 E 1 ). This is a key ingredient of the closed-form method, i.e., working with the angles f and E 1 , instead of the mean anomalies M, M 1 , no series reversion of Kepler's equation is used throughout the whole perturbative scheme.
In order to avoid relegation, our method discussed below works locally, by constructing a model for the secular Hamiltonian valid for a particle's semi-major axis varying as a = a * +δa(t), i.e., by a small quantity δL around some reference value a * . By standard secular theory, we have the estimate δa = O(µ) far from mean-motion resonances. Formally, introducing the new canonical variable δL as and expanding the Hamiltonian in powers of the quantity δL around L * , we obtain where a constant term −G 2 m 2 0 /(2L 2 * ) was dropped from the expansion. The constant n * = G 2 m 2 0 /L 3 * is equal to the particle's mean motion under Keplerian orbit at the semi-major axis a * . Remark 3. The choice of the reference value a * determines the kind of divisors appearing in the normalization procedure. In the present paper, we deal only with the 'non-resonant' case, in which the frequencies n * and n 1 satisfy no-commensurability condition. For example, to be far from any resonance we may require that n * and n 1 satisfy a diophantine condition with |k| = |k * | + |k 1 | and some suitable γ > 0, τ > 1. However, the algorithm presented below can be readily extended to cases of mean-motion resonance. We leave the details for a future work, noting only that in resonant cases we have the estimate δL = O(µ 1/2 ), instead of O(µ). The effect of approaching close to a mean-motion resonance with the present series is seen, instead, as a rise in the value of the series' remainder, caused by (non-zero) small divisors in the series (as visible, for example, in Fig. 7 discussed in section 3 below).

Poisson structure and book-keeping
implemented to the closed-form version of the functions F 1 , F 2 . The closed-form version of a function F is defined as: The derivatives in the canonical variables of a function F as in Eq.(21) are computed by the chain rule formulas where A sketch of the derivation of the above formulas can be found in Appendix A. They are strictly valid with e ∈ (0, 1), i ∈ (0, π). However, several cancellations lead to no singular behavior of the Poisson bracket formulas arising throughout the various perturbative steps also when e = 0 or i = 0.

Book-keeping: Hamiltonian
We introduce in the series a book-keeping symbol σ (see [5] for an introduction to the bookkeeping technique), with numerical value σ = 1, whose role is to provide a grouping of all the various terms in the series according to their 'order of smallness'. Hence, a group of terms with common factor σ l , l ∈ Z, indicates a term considered as of the 'l-th order of smallness'. Since in our series there are several small quantities, we introduce a book-keeping scheme allowing to simultaneously deal with all small quantities while maintaining the closed-form character of the series. To this end, we make the following substitutions, called 'book-keeping rules', within the initial Hamiltonian: • BK-Rule 1: e σ 1 e = σe (not applicable to the quantity e 2 within η = √ 1 − e 2 ), • BK-Rule 2: η σ 0 η = η, • BK-Rule 3: µ σ ν µ, with ν as in Eq. (7), • BK-Rule 4: e 1 σ ν 1 e 1 , with ν 1 as in Eq. (7) (not applicable to the quantity e 2 1 within • BK-Rule 5: 1 Since σ = 1, the above substitutions affect the structure of the series only at the formal level, and can be substituted directly into the original Hamiltomian, whereby they propagate at subsequent normalization steps once these steps are organized in successive powers σ, σ 2 , etc., of the book-keeping symbol. The BK-Rules 1 to 7 above are justified on physical ground as well as on motives of algorithmic convenience. In particular: -BK-Rule 1 implies that, despite the use of closed-form formulas, the basic small quantity in powers of which the series are organized is the eccentricity of the test particle. -BK-Rule 3 implies that a factor µ in front of a series term should be treated as of comparable order of smallness as a term of order e ν , with ν given by Eq. (7). Similarly, BK-Rule 4 implies that a term containing a factor e 1 raised to some power should be treated as of comparable order of smallness with a term e ν 1 raised to the same power. Note that the eccentricity e is a quantity variable in time, so that to compute the exponents ν, ν 1 we need to use, for any examined trajectory, a reference value e * yielding an estimate of the overall level of eccentricity all along the orbital evolution for that trajectory. Note that, by standard secular theory we have e(t) = e * + O(µ) if e * is close to the mean eccentricity (see also discussion at the introduction). Note finally that we obtain exponents ν, ν 1 ≥ 1 in the typical case in which e > µ and e ≥ e 1 . These inequalities arise naturally in the case of small bodies in highly eccentric orbits perturbed by some planet of, say, our solar system, which are the cases of main interest in applying the present method (see, nevertheless, Remark 4 on the treatment of cases where the above conditions are not met).
-BK-Rule 7 stems from the estimate δL = O(δa) = O(µ) holding for the oscillations in semimajor axis of trajectories far from mean-motion resonances (as already pointed outin the latter case, instead, we have in general δL = O(δa) = O(µ 1/2 ) and the corresponding rule has to be adapted accordingly). The lowering of the book-keeping power by one for within H 0 is introduced for reasons of algorithmic convenience, i.e., in order to maintain n * δL in the kernel of the homological equation.
-BK-Rules 5 and 6 imply just a partition of the unity aiming at keeping the perturbative scheme in closed-form while splitting the corresponding expressions (involving η and η 1 respectively) in two parts, of orders O(1) and O(e 2 ), or O(e 2 1 ).

Book-keeping: Poisson structure
Some of the formulas in Subsection 2.3.1 imply differentiation with respect to e through the corresponding partial derivatives in (27), (28), thus yielding a lowering of the power of the eccentricity in some terms arising through Poisson brackets at consecutive steps of perturbation theory. To account for this fact, similarly as in [1] we introduce the use of the book-keeping symbol σ in the formulas of the Poisson algebra as follows: first, we re-write the derivatives with respect to the angles , g, h, M 1 as and with respect to the actions δL, G as Note that in (47) use was made of the identity φ 1 = e 1 sin E 1 (Kepler's equation). Finally, we revise formulas (31), (32), (34)-(43), attributing a book-keeping to all factors involving the eccentricity function η as Remark 4. The small eccentricity problem consists of the fact that the above-proposed bookkeeping rules are not applicable in the case 0 < e * µ < e 1 , since, by (7), the exponents ν, ν 1 would be smaller than unity. The simple solution of rounding these exponents to 1, while maintaining the same book-keeping rules as above, fails, since, at any given normalization order r, the presence of σ −1 , σ −ν 1 terms in the formulas of the Poisson algebra leads to the generation of terms of order lower than r in the normal form's remainder. Notwithstanding our focus on a method dealing with large eccentricity orbits (for which the problem does not appear), we discuss below a variant of the main algorithm that deals with trajectories in the case ν = 1, i.e., when e * µ.

Preliminary step: Hamiltonian preparation
After implementing BK-Rules 1 to 7 the Hamiltonian (19) resumes the form: where σ s ∈ {σ ν , σ ν+1 , . . .} and, by D'Alembert rules, only cosines and real coefficients q s appear (invariance under simultaneous change of sign of all angles). Setting Z 0 = n * δL + n 1 J 1 , for obtaining a closed-form normalization algorithm it turns convenient to re-express the Hamiltonian according to The Hamiltonian (63) resumes the form: where ν the remainder at the zero-th normalization step (i.e. in the original Hamiltonian). The terms R (0) ν,l contain terms of book-keeping order σ l , with l ≥ ν.

Step 1: normalization of the σ ν -terms
For a suitable generating function χ (1) ν to be determined in a while, we introduce the Lie series operator as where C ω (T 4 × D) denotes the set of real analytic functions in the phase space and is the time derivative along the Hamiltonian vector field generated by χ in which, with the usual abuse of notation, we still indicate with , g, h, M 1 , δL, G, H, J 1 the new canonical variables given by the inverse transformation Our scope will be to define the Lie generating function χ (1) ν in such a way that, after implementing the transformation (68), H (1) contains no terms depending on the angles f and E 1 at order σ ν . The required generating function χ (1) ν is computed as an outcome of the following: Then, it holds that where Furthermore, the function H (1) as computed by Eq.(68) takes the form where the remainder Proof. Setting and recalling the chain rules (44), (47) and (50), (51), (33), we find Requiring that no trigonometric terms depending on f, E 1 be present at order σ ν then leads tô which implies Eq.(70). At order σ ν we then obtain immediately the formula We now consider the function H (1) computed by replacing (70) into (68). The function H (1) can be decomposed as in Eq.(73). We shall demonstrate that the remainder R (1) contains no terms of order lower than σ ν+1 . To this end, it suffices to show that ν contains terms of order equal to or larger than σ ν , while χ (1) ν contains only terms of order σ ν . Thus, except for the Poisson bracket {Z 0 , χ (1) ν }, which only contributes to the secular terms Z (1) ν due to Eq.(71), the first Poisson bracket in (74) contains prefactors of order σ 2ν or higher, while the second contains prefactors σ nν or higher. However, the exponent of σ in these brackets can be lowered due to the negative powers introduced in the book-keeping formulas in the following three classes of factors: (i) partial derivatives with respect to the eccentricity in (48), (49) (carrying σ −1 ) multiplied by corresponding formulae (53), (56) (another σ −1 ), hence a total of σ −2 ; (ii) differentiations (52), (55) involving f (weighting σ −1 ) again in (48), (49), thus a pre-factor σ −1 ; (iii) partial derivatives with respect to φ 1 in (47) (σ −ν 1 , ν 1 ≥ 1), thus a prefactor at least σ −1 .
As regards (i), we first note that χ (1) ν has no explicit dependence on e, but only an implicit dependence through η, which in the closed-form context is treated as an independent symbol. This follows from the fact that χ (1) stems from balancing the coefficients of R (0) ν,ν . The latter term contains a pre-factor µ, which is already O(σ ν ), thus it cannot contain any further factors produced by any explicit power of e. In view of the above, setting ∂χ (1) /∂e = 0, we find that for any F ∈ C ∞ (T 4 × D) the expression in {F, χ (1) ν } pertaining (i) can be factored out as We now have the following lemma: Lemma 1. For every term in the Hamiltonian (63) of the form i.e., explicitly independent on e, we have s 1 = s 2 .
Proof. This is a consequence of D'Alembert rules. Using modified Delaunay angular elements as well as the formulas f = + 2e sin + O(e 2 ), eη(e) −2λ = e + λe 3 + O(e 5 ), λ ∈ N, we find that, after expanding in the eccentricity e, (76) should give the terms However, according to the D'Alembert rules, in a generic trigonometric monomial of the form appearing after expanding H in the eccentricities e, e 1 , we necessarily have that l − |w 2 | must be non-negative and even. Since for any closed-form term in the Hamiltonian, explicitly independent of e, the lowermost term in e produced after the expansion satisfies l = 0, we necessarily have In view, now, of (70), the relation s 1 = s 2 implies ∂χ ν /∂g. Therefore, making use of (50), (53) and (56), Eq.(75) translates into It follows that for any of the functions F = R ν }, . . ., terms produced by derivatives of the type (i) in (68) are subject to a lowering of the exponent of σ per Poisson bracket only by a factor σ −1 , instead of σ −2 . In particular, in the case F = R (0) ν,ν (as well as for any other closed-form function explicitly independent on the eccentricity) we have that (75) is identically vanishing.
As regards (ii), we find that for any F 1 , F 2 ∈ C ∞ (T 4 × D), the derivative ∂f /∂δL (Eq.(52)) participates in the Poisson bracket {F 1 , F 2 } only through the combination On the other hand, the derivative ∂f /∂G (Eq.(55)) participates in the same Poisson bracket through the combination ∂f ∂G which, by Lemma 1, is also equal to zero for F 1 = R (0) ν,ν (or any other term O(σ ν+1 ) in H not depending explicitly on e), and F 2 = χ (1) ν .
In conclusion, returning to (74), and taking all the above deductions into account, we arrive at the expressions and similarly, ν } satisfies Lemma 1. We then have {Z 0 , χ which concludes the proof of the proposition.
These last algebraic operations conclude the first normalization step.

Loop: normalization of the σ ν+j−1 -terms
The procedure followed in the first step can be repeated iteratively in order to normalize consecutively terms of order σ ν+j−1 , with each time an O(σ ν+j ) remainder, for ν, j > 1. As anticipated in Remark 4, the iterative procedure described below fails in the case ν = 1 at step j = 2, so an adjustment (involving one more iteration) is required, as discussed in Subsection 2.4.4 below.
The j-th normalization step is carried out as follows from the next proposition.

The case ν = 1
Coming to ν = 1, one realizes that (91) produces same order σ 2 non-normalized terms via {R 2 }, χ 2 }, . . . , χ 2 }, namely the resulting remainder is R 2 , so the scheme in Proposition 2 is not directly applicable beyond j = 1. Despite this, it is worth noticing that if we manage to get rid of these spurious terms, by performing, for instance, an extra normalization II, such that the new outcome returns R (II) = R (II) 3 , then the algorithm (87) will work for j ≥ 3 upon restarting the recursion from iteration II in place of 2. This is precisely the claim we are about to show to complete the treatment. Let us write (91) as H (2) 2 . Introduce the extra second normalization II based on Proposition 2 targeted to R Moreover the loop composed by (87)-(90) in Proposition 2 holds true for any j ≥ 4 under the modifications Proof. We begin with a necessary generalization of Lemma 1.
Lemma 2. Given F 1 , F 2 ∈ C ω (T × D) trigonometric monomials of the form (76), or equivalently in terms of the sine, fulfilling the property of Lemma 1, addenda of the same type in the Lie series transformation applied to F 1 with respect to F 2 preserve such property.
Proof. Since exp (L F 2 ) F 1 involves the computation of Poisson brackets of functions explicitly independent on e, we have that (92), with F 1 , F 2 in place of F, χ ν+1 , is identically null, as well as (81) because ∂F 1 /∂f = ∂F 1 /∂g, ∂F 2 /∂f = ∂F 2 /∂g by assumption. Thus, the bracket {F 1 , F 2 } in the Lie series either does not introduce any eccentricity dependence at all, or only at numerator through (50) multiplied by cos f or cos 2 f ; therefore its derivatives contain products of cosines (sines) whose coefficients are independent on e like The arguments are now either summed or subtracted, hence they clearly satisfy the property concerned. By cascade reasoning for further nested brackets we conclude. We consider step II: The analysis of the contributions reports these deductions, by which (94) follows.
is independent on f, g, e.
• {Z 2,2 depends on e at most linearly by book-keeping rules, so it does χ  Proof. By Proposition 1, Lemma 1 and 2, the substitution φ 1 = e 1 sin E 1 and the formulas listed in Subsection 2.3.1, Subsection 2.3.3, we take out of (68) the order σ 2 remainder and it is not restrictive to assume ν 1 = 1 in order to include also the e 1 cos E 1 dependent term in (71): where We employ now all D'Alembert rules to show that only the harmonics of interest can exist. Following the same argument as in Lemma 1, let us write the cosine input of (99) using modified Delaunay angles (77) also for P 1 in relation to corresponding orbital elements (3) (subscript '1'): in whichp 1 =q 1 = 0. For the elimination of the apparent singularity at e = 0, we must have 1 − |s 1 + s 2 | ≥ 0 and even, hence s 2 = s 1 ± 1. Then, since R 1,2e is independent on e 1 by book-keeping setting, analogously we must end up with s 4 = s 5 . Regarding instead the regularity at i 1 = 0, because of the absence of i 1 we must conclude that 0 − |s 5 − s 6 | ∈ 2N, namely s 5 = s 6 . At this stage, we invoke the invariance under rotation around the Z axis, which prescribes s 1 − s 1 + s 2 − s 2 + s 3 + s 4 − s 4 + s 5 − s 5 + s 6 = s 3 + s 6 = 0 , and summing up this implies s 3 = −s 4 . Ultimately, concerning the inclination, we must ensure that l − |s 2 − s 3 | ∈ 2N, with l even as well again being i 1 not involved, thus s 2 = s 3 ± 2n, n ≤ l/2 natural number. Putting all together we arrive at which always depends on at least one among f, E 1 since the coefficients s 1 , ±2n ∓ 1 − s 1 never vanish simultaneously. By means of an identical reasoning and given the preservation of D'Alembert rules under exp L χ (1) 1 , we achieve the same outcome for the remaining part of (98) after replacing s 1 → 1 ± s 1 , indeed we find and no solutions to 1 ± s 1 = 0, ±2n ∓ 1 − 1 ∓ s 1 = 0.
Given that the order 2 normal form is sourced from the part of R (1) 2,2 explicitly independent on fast angles, it turns out that it is free of e. Finally, R In order to conclude, we just need to check that the next step gives rise to an O(σ 4 ) perturbation and the cycle of normalizations can restart for j ≥ 4 in light of the bounds on σ from (93) at the end of the proof of Proposition 2. Upon repeating the usual argument, it is easy to see that the only bracket worth investigating is {Z 2,2 independent on e.
Serving as an example, a detailed demonstration of the normalization procedure exposed in the present section for a simple model, containing just few terms of the disturbing function, is presented in Appendix B.

Computer-algebraic implementation of the normalization algorithm
Implementing the above normalization procedure, e.g. by use of a Computer Algebra System (CAS), requires working with a finite truncation of the initial Hamiltonian model (11). To this end, the disturbing function (13) multiplied by µ can be re-arranged as whereh κ 1 ,κ 2 ,κ 3 are real coefficients derived from the coefficients of (13). A convenient truncation of (100) stems from defining two separate truncation orders in powers of µ (truncation order k µ ), and in powers of r 1 / R (multipole truncation order k mp ), through the formula where · is the integer part function. Working with the truncated Hamiltonian H ≤kµ,kmp = H 0 + H ≤kµ,kmp 1 , we then obtain a sequence of secular models Z (j) , j = 1, 2, . . ., where j denotes the normalization step, computed via the formula In particular, we implement the following steps of the CAS algorithm: (i) for a fixed value of µ, choose values for k µ , k mp , perform the corresponding expansions of the Hamiltonian as in (100) and compute the truncated model H ≤kµ,kmp ; (ii) choose the reference values of a * and e * ; (iii) pass to variables (f, g, h, E 1 , δL, e, η, ι c , ι s , J 1 ) and parameters L * , e 1 , a 1 , η 1 on the basis of the selected a * ; (iv) compute ν and ν 1 (Eq. (7)); (v) set the appropriate book-keeping weights following the rules in Subsection 2.3.2 and expand correspondingly the Hamiltonian in δL up to σ νkµ ; (vi) drop constants, perform the identity operation (63), discard book-keeping powers larger than νk µ and introduce n * ; (vii) if ν > 1, compute the generating function (70)  (viii) compute the successive normalizations H (j) , truncated at book-keeping order N bk via the procedure of Subsection 2.4.3, up to a maximum normalization order ν + j max − 1 < N bk , j max ≤ ν(k µ − 1); this allows to obtain truncated Hamiltonian models containing a finite number of normal form terms as well as a finite number of terms provided by the truncated remainder.
In the CAS implementation of the above algorithm we work with numerical coefficients, substituting all constants with their corresponding numerical values. Several types of numerical tests of the precision and overall performance of the method can be carried out as exemplified in the sequel.
Our basic probe of the efficiency of the normalization method in the framework of the ER3BP is given by comparing the short-period oscillations of the orbital elements a(t), e(t), i(t), g(t), h(t), as found by two different methods.
Semi-analytical propagation: following the implementation of the normalization algorithm as described in the previous subsection, the initial osculating element state vector z(0) is transformed  into an initial condition for the corresponding 'mean element' state vector ξ (j) (z(0)), i.e., the element vector corresponding to the new canonical variables conjugated to the original ones after j near-identity normalizing transformations. This is computed by the Lie series composition formula truncated at book-keeping order N bk : using Eq.(69) for the inverse series. We then obtain the evolution of the mean element vector ξ (j) (t) through numerical integration of the secular equations of motioṅ (J standard symplectic unit). This can be back-transformed to yield the evolution of the osculating element vector z(t) using the truncated Lie series composition formula Note that both the direct and inverse transformations (Eqs. (103) and (105)), as well as Hamilton's secular equations (104), can be computed in closed form, using the Poisson algebra rules of Subsection 2.3. We then call semi-analytic the evolution of the element vector z(t) obtained via the formula z(t) = z(ξ (j) (t)) . . In these cases, the distance ratio r 1 / R is small (about 0.1-0.2), a fact implying that the quadrupolar expansion (k mp = 2) suffices to have obtained a relative error of about 0.1% in the representation of the Hamiltonian perturbation H 1 . Going to higher multipoles is straightforward, albeit with a significant computational cost as the number of terms in the Hamiltonian grows significantly. On the other hand, even with low-order truncations of the Hamiltonian we achieve to have an accurate semi-analytical representation of the O(µ) short-period oscillations in all three 'actionlike' elements (semi-major axis, eccentricity, inclination). Most notably, keeping a(0) = 50AU but changing the eccentricity to e(0) = 0.7, i.e., beyond the Laplace value, yields an orbit whose pericenter is at R p = 15AU, implying a distance ratio r 1 / R ≈ 0.3 (Fig. 3). This time, an octupole truncation (k mp = 3) is required to produce an approximation of the Hamiltonian model at the level of a relative error of 0.1%. Still, however, as shown in Fig. 3 the semi-analytical propagation of the orbit is able to track the fully numerical one with an error which does not exceed 0.2% even close to the orbit's pericentric passages.
In the above examples, the maximum number of normalization steps at which the secular Hamiltonian is computed was set equal to j max = 3, j max = 4 and j max = 4 respectively, which corresponds to the best match in all cases. As discussed in the next subsection, an estimate of the minimum possible error in the semi-analytic propagation of the trajectories requires computing first the so-called optimal number of normalizations j opt (or equivalently optimal normalization order ν + j opt − 1) as a function of the reference values (a * , e * ) within a model given by a preset fixed multipole truncation order. Owing to the fact that the same divisors appear in the ER3BP and in the CR3BP, we verify with numerical examples that the error analysis yields essentially identical results in either case. However, the computation of the optimal normalization is easier to perform in the CR3BP, owing to the considerably smaller number of terms produced in the CAS computation of the normal form. Hence, we now turn our attention to this latter computation.

Trajectory propagation: optimal remainder
A considerable reduction of the computational cost occurs in the case of the planar and circular R3BP. This is due, in particular, to the following: • the dependence on M 1 becomes explicit (M 1 = E 1 in (2)), while a 1 = r 1 . As a consequence, φ 1 = 0.
• no terms involving (h, H) appear in the disturbing function, thus ι c , ι s are discarded; • no terms requiring a book-keeping in terms of the exponent ν 1 appear, hence, only ν is defined, as in (7); l,λ,p = 0 for every j, l, λ, p in (90), (86), and consequently p 1 = p 2 ≡ 0 in (88). This is due to the fact that the expression (16) reduces to which always depends on the difference g − M 1 by D'Alembert rules. This implies that, unlike the ER3BP, the action G (and the corresponding eccentricity e) are integrals of the secular Hamiltonian; • as a consequence no lower or equal book-keeping order terms appear in any Poisson bracket of the first normalization step in the case ν = 1. Hence Proposition 3 is redundant.
Owing to the above, in the planar CR3BP we are able to make normal form computations in a grid of points in the plane (a * , e * ) up to a sufficiently high normalization order so that the asymptotic character of the series computed by the algorithm of Section 2 can show up. To this end, we introduce an estimate of the size of the series' remainder after j normalization steps via the upper norm bound where · ∞ denotes the sup norm. Plotting E (j) against the number of normalization steps j allows then to estimate the error committed at any step (size of the remainder). Figure 4 yields an example of such computation. The relevant fact is that there is an optimal number of normalization steps (j = j opt = 6) where the estimate E (j) of the remainder size yields a global minimum.
Although a systematic investigation of the dependence of the optimal number of normalization steps j opt on the parameters (a * , e * ) is beyond our present scope, Figs. 5 and 6 allow to gain some insight into the question. The most relevant remark concerns the dependence of the behavior of the curve E (j) (versus j) on how close to the 'hierarchical' regime the trajectory with reference values (a * , e * ) is. As a measure of the hierarchical character of an orbit we adopt either the ratio of the semi-major axes a 1 /a * , or of the pericentric distances r 1 / R p = a 1 (1 − e 1 )/(a * (1 − e * )) = a 1 /(a * (1 − e * )). Fig. 5 (a * = 30AU, e * = 0.5) implies a pericentric distance ratio r 1 / R p ≈ 0.3 smaller than the one of the example of Fig. 4 ( r 1 / R p ≈ 0.4). We observe that the optimal number of normalization steps in the former case satisfies j opt = 10, i.e., it is larger than in the latter case. Fig. 6 shows, instead, an example of orbit far from the hierarchical limit, satisfying the estimate r 1 / R p ≈ 0.7. In this case a higher order multipole expansion 1 2 3 4 5 6 7 8 (k mp = 5) is required to obtain a precise truncated Hamiltonian model for this orbit. We note, however, that the normalization procedure performs well, producing a decreasing remainder as a function of j up to the point where it is arrested, i.e. j = 6 = ν(k µ − 1). We find numerically that this performance is deteriorated as we gradually approach the condition r 1 / R = 1, beyond which the multipole expansion of the Hamiltonian is no longer convergent.

Semi-analytical determination of the domain of secular motions
The results shown in the two previous subsections refer to isolated examples of orbits treated within various multipole truncation orders as well as different choices of the number of normalization steps, searching each time to arrive at the best approximating secular model given computational restrictions. In the present subsection, we aim to investigate the behavior of the remainder in a closed-form normalization with uniform choice of all truncation orders of the problem, but performed, instead, in a fine grid (100 × 20) of reference values in the plane (a * , e * ). To this end, we set k µ = 2 (second order in the mass parameter), and fix k mp = 3 (octupole approximation). The latter choice, imposed by computational restrictions, yields an initial model whose error with respect to the full Hamiltonian becomes of the order of 1% only for a * > 2a 1 . However, for reasons explained below, a computation within the framework of the octupole approximation becomes relevant to the problem addressed in the sequel also in the range 1.5a 1 < a * < 2a 1 , while higher multipoles are required to address still smaller values of a * . The result of the above computation is summarized in Fig. 7: the left panel shows in logarithmic color scale the size of the remainder, estimated by the value of E (n) (a * , e * ) computed as in (108), corresponding to each point in the plane (a * , e * ), where the number of normalization steps is set as n = min{ν(k µ − 1), 7} = min{ν, 7}. The maximum value n = 7 is, again, imposed by computational restrictions, and it implies that n varies with e * up to about e * = 0.37.
The relevant information in Fig. 7 is provided by the black curve, which corresponds to the isocontour E (n) (a * , e * ) = 10 −2 . Since in the original Hamiltonian we have the estimate E (0) (a, e) := H ≤kµ,kmp 1 = O(10 −2 ), the black curve provides a rough estimate of the limiting border dividing the plane (a * , e * ) in two domains: in the one below the black curve the progressive elimination of the fast angles by the iterative normalization steps leads to a secular model whose remainder decreases with the number of normalization steps j at least up to j = n.
A physical interpretation of the border approximated through the isocontour E (n) (a * , e * ) = 10 −2 can be given through a comparison with a numerical stability map obtained, e.g., as in the right panel of Fig. 7. For each trajectory in a 300 × 900 grid in (a, e), the plot shows in color scale the value of the Fast Lyapunov Indicator (FLI, see [9] for a review) obtained after integrating the variational equations of motion together with the equations of motion of the full Hamiltonian model for a time equal to 50 periods of Jupiter. Thus, deep blue colors indicate the most regular, and light yellow the most chaotic orbits as identified by the value of the FLI. Superposed to the FLI cartography are three curves: For every e = e * , n different normalizations are executed and then evaluated for each a = a * . Right panel: short-period FLI map over a 300 × 900 (a, e) grid of initial data integrated for 50T 1 . As indicated, the three curves represent, respectively, the line of constant pericenter of the particle's trajectory equal to the radius of Jupiter's orbit r 1 = r J (red), Hill's stability criterion (brown) and the isolevel E (n) = 1% (black). Each region enclosed by two consecutive above curves is labeled with the corresponding regime of motion. The main mean-motion resonances are reported below the pictures.
(i) the 'perihelion crossing curve' (red) yields the locus of values satisfying the condition a(1 − e) = r J = a J (in the circular case), that is the points where the pericenter of the test particle's orbit comes at distance equal to the radius of Jupiter's orbit; (ii) the Hill limit [12] (brown) is based on the relationship C Jac (a, e) = C Jac (L 1 ), where C Jac is the particle's Jacobi constant as function of the orbital elements and C Jac (L 1 ) its value at the Lagrangian point L 1 ; (iii) the isocontour E (n) (a, e) = 10 −2 (black, same as in the left panel of Fig. 7).
Of the above three curves, the perihelion crossing curve is analogous, in the R3BP, of the so-called Angular Momentum Deficit criterion (AMD, [8]) used to separate systems protected from perihelia crossings in the case of the full planetary three-body problem. As indicated by the FLI cartography data, Hill's curve gives an overall better approximation separating the domain of strong chaos (yellow) from the domain of regular or weakly-chaotic orbits (all blue nuances). This is expected, since the Hill's curve separates orbits for which Jupiter's gravitational effect becomes (at least temporarily) dominant from those for which it does not. Nevertheless, through the FLI cartography we note the presence of a large domain between the curves (ii) and (iii), where the trajectories, while protected from close encounters, are subject to the long term effects on dynamics produced by resonant multiplets associated with the mean-motion resonances of the problem (the most important of which are marked in the figure). Note that in the octupole approximation, the Hamiltonian contains harmonics including all combinations of the fast angles of the form cos(s 1 f + s 2 (g − M 1 )), with thus including all harmonics associated with the mean-motion resonances detected in the FLI cartography of Fig. 7 for a > 1.5a J . Through the closed-form normalization (Eqs. (70) and (86)) we then obtain small divisors in the series at every value of the semi-major axis a * for which one of the resonant combinations s 1 n * − s 2 n J , n J = n 1 , takes a value near zero. All these incidences lead to Arnold tongue-like spikes pointing downwards in the curve (iii), marking the failure of the approximation of the orbits based on a non-resonant normal form construction. On the other hand, we observe that, for any value of a * there is a threshold value of the eccentricity e * ,s , such that, for e * < e * ,s no visible effects of the harmonics associated with mean-motion resonances are visible in the FLI cartography. This implies that the secular models constructed by eliminating all harmonics involving the fast angles of the problem describe with good precision the dynamics in this domain, called, for this reason, the domain of secular motions. In physical terms, the domain of secular motions corresponds to initial conditions for which the gravitational perturbation of Jupiter is only felt in the 'Laplacian' meaning, i.e., as a mass distributed along a ring coinciding with Jupiter's orbit. The curve (iii) then yields the limit of this domain, which, as found by the FLI cartography, is well distinct from the limit of the Hill domain. The overall situation can therefore be summarized with the identification of four regimes of motion (specified in the FLI chart): • the 'crossing orbit regime' (above curve (i)); • the 'close encounter regime' (between curves (i) and (ii)); • the 'resonant regime' (between curves (ii) and (iii)); • the 'secular regime' (below curve (iii)).

Conclusions
In summary, in the present paper we have proposed a closed-form method for the derivation of secular Hamiltonian models (normal forms) with a small (albeit finite minimum) remainder applicable to the R3BP in the case when the particle's trajectory is exterior to the trajectory of the primary perturber. Also, using this method we were led to the definition of a new heuristic limit separating the motions whose character is 'secular', i.e., not affected by short-period effects, from the rest of motions in the R3BP. In particular: 1. Section 2 develops the formal aspects of the method, which heavily relies on the use of a book-keeping parameter to simultaneously account for all small quantities of the problem as they appear not only in the Hamiltonian and Lie generating functions, but also in the closed-form version of all formulas involved in the Poisson algebra between the Delaunay canonical variables of the problem. A rigorous demonstration of the consistency of the method is then given through Propositions 1, 2 and 3, which also estabilish the explicit formulas for the implementation of one iterative step of the closed-form normalization algorithm.
2. Section 3 gives numerical examples of the implementation and precision of the algorithm in the spatial elliptic, as well as in the planar circular R3BP, examining, also numerically, the method's convergence properties. The effect of choosing different truncation orders (in powers of the mass parameter µ or in the multipole expansion) is discussed, along with several simplifications to the normalization procedure which hold in the circular case. The essentially asymptotic character of the series is established through numerical examples, showing the existence of an optimal number of normalization steps, after which the size of the remainder becomes the minimum possible.

3.
A key aspect of the above presented method lies in the possibility to exploit the behavior of the size of the remainder as a function of the number of normalizing steps in order to obtain a clear separation of two well-distinct domains, as also identified by purely numerical (FLI cartography) means: one, called the domain of secular motions corresponds to the domain where the harmonics in the Hamiltonian associated with resonant combinations of the fast angles (anomalies) of the problem produce no dynamical effect on the orbits visible at the level of the FLI cartography. From the semi-analytical point of view, this turns to be the domain where a non-resonant construction as the one proposed in section 2 produces no (nearly-)resonant divisors up to the optimal normalization step. As a consequence, only the angles associated with the motions of the perihelion and of the line of nodes survive in the final normal form. We show numerically how to use the information on the size of the normal form remainder in order to determine semi-analytically the border of the domain of secular motions in the case of the Sun-Jupiter system. We finally give evidence that this border is well distinct from the border of the domains defined either by the Hill stability or by the perihelion crossing criterion.

B Example of normalization for a µ 2 quadrupolar expansion
Consider the following toy model Hamiltonian with k µ = k mp = ν = 2, ν 1 = 1, according to conventions introduced in §2.4.1: Next, we move on with the second and last iteration j = 2 targeted to R