Symbolic dynamics in a binary asteroid system

We highlight the existence of a topological horseshoe arising from a a--priori stable model of the binary asteroid dynamics. The inspection is numerical and uses correctly aligned windows, as described in a recent paper by A. Gierzkiewicz and P. Zgliczy\'nski, combined with a recent analysis of an associated secular problem.


Purpose of the paper
This paper aims to highlight chaos in the secular motions of a binary asteroid system interacting with a planet whose orbit is external to the orbits of the asteroids. These chaotic motions turn to bifurcate from an a-priori stable configuration, in the sense of Chierchia and Gallavotti (1994). We shall not provide rigorous proofs, besides the heuristic arguments that we are going to present in this introduction. In fact, our study will be purely numerical. Moreover, we shall not implement any algorithm to control machine errors. We are however convinced that our computations are correct thanks to a-posteriori checks that we shall describe in the course of the paper. Let us describe the physical setting. Three point masses constrained on a plane undergo Newtonian attraction. Two of them (the asteroids) have comparable (in fact, equal) mass and, approximately, orbit their common barycentre. The orbit of a much more massive body (the planet) keeps external to the couple, for a sufficiently long time. We do not assume 1 any prescribed trajectory for any of the bodies, but just Newton law as a mutual interaction. We fix a reference frame centred with one of the asteroids and we look at the motions of the other one and the planet. As no Newtonian interaction can This research is funded by the ERC project 677793 StableChaoticPlanetM. 1 See, e.g., Paez and Locatelli (2015) for a study based on a restricted model. be regarded as dominant -as, for example, in the cases investigated in Arnold (1963); Féjoz (2004); Laskar and Robutel (1995); Pinzari (2009); Chierchia and Pinzari (2011) and Giorgilli et al. (2017); Volpi et al. (2018) -in order to simplify the analysis, we look at a certain secular system, obtained, roughly, averaging out the proper time of the reference asteroid. This means that we are assuming that the time scale of the movements of the planet is much longer. Beware that our secular problem has nothing to do with the one usually considered in the literature, where the average is performed with respect to two proper times (e.g., Féjoz and Guardia (2016)). Let us look, for a moment, to the case where the planet is constrained on a circular trajectory. In such case, the only observables are the eccentricity and the pericentre of the instantaneous ellipse of the asteroid. Quantitatively, this system may be described by only two conjugate Hamiltonian coordinates: the angular momentum G (related to the eccentricity) and the pericentre coordinate g of the asteroidal ellipse. There is a limiting situation, which roughly corresponds to the planet being at infinite distance, where, exploiting results from Pinzari (2019Pinzari ( , 2020a, the phase portrait of the system in the plane pg, Gq reveals only librational periodic motions. Physically, such motions correspond to the perihelion direction of the asteroidal ellipse affording small oscillations about one equilibrium position, with the ellipse highly eccentric and periodically squeezing to a segment. The movements are accompanied by a change of sense of motion every half-period. The purpose of this paper is to highlight the onset of chaos in the full secular problem, when the planet is far and moves almost circularly.
The Hamiltonian governing the motions of three point masses undergoing Newtonian attraction is, as well known, Here, x 0 , x 1 , x 2 and y 0 , y 1 , y 2 are, respectively, positions and impulses of the three particles relatively to a prefixed orthonormal frame pi, j, kq Ă R 3 ; m 0 , m 1 " µm 0 , m 2 " κm 0 , with y i " m i 9 x i , are their respective gravitational masses; |¨| denotes the Euclidean distance and the gravity constant has been taken equal to one, by a proper choice of the unit system. In the sequel, in accordance to our problem, we shall take x i , y i P R 2ˆt 0u » R 2 and µ " 1 ! κ, so that x 0 , x 1 correspond to the position coordinates of the asteroids; x 2 is the planet. The Hamiltonian H is translation invariant, so we rapidly switch to a translation-free Hamiltonian by applying the well known Jacobi reduction. We recall that this reduction consists of using, as position coordinates, the centre of mass r 0 of the system (which moves linearly in time); the relative distance x of two of the three particles; the distance x 1 of the third particle with respect to the centre of mass of the former two. Namely, Note that, under the choice of the masses specified above, we are choosing the asteroidal coordinate x 0 as the "starting point" of the reduction. This reverses a bit the usual practice, as x 0 is most often chosen as the coordinate of the most massive body; see Figure 1. At this point, the procedure is classical: the new impulses pp 0 , y, y 1 ) are uniquely defined by the constraint of symplecticity, with p 0 ("total linear momentum") being proportional to the velocity of the barycentre. Choosing (as it is possible to do) a reference frame centred at, and moving with, r 0 , so to have r 0 " 0 " p 0 , after a suitable rescaling, one obtains (see Appendix A) The choice κ " µ " 1 gives β "β " 1 and simplifies H to Schematic representation of the model we are dealing with. The model is composed by three bodies P 0 , P 1 , P 2 , where the first two have equal masses m 0 and the third body has the largest mass κm 0 , κ ą 1. The point b 1 is the barycentre of P 0 and P 1 , whilst b 2 is the barycentre of all the three points (close but different from P 2 ).
From now on, we regard β as mass parameter, with β " κ and σ " β 2 . By choosing a region of the phase-space where |x 1 | ą |βx|, (1.6) we ensure the denominators of the two last terms in (1.5) to be different from zero. The Hamiltonian (1.5) with x, x 1 , y, y 1 P R 2 , has 4-degrees-of-freedom (DoF, from now on), but is SO(2)-invariant. We choose a system of canonical coordinates which reduces this symmetry and hence lowers the number of DoF to 3. If k " iˆj is normal to the plane of the orbits, we denote as the total angular momentum, which is a constant of the motion. Then we take a 3-DoF system of coordinates, which we name pΛ, G, R, , g, rq P R 3ˆT2ˆR , where pΛ, G, , gq are "Delaunay coordinates for the asteroid, relatively to x 1 ", while pR, rq are "radial coordinates for the planet". More precisely, they are defined as where, considering the instantaneous ellipse generated by the first two terms in the Hamiltonian (1.5), a is the semi-major axis (see again Figure 1), S and S tot are the area of the ellipse spanned from the perihelion P and the total area and α x 1 ,P is the angle between the direction of x 1 and P relatively to the positive direction established by xˆy. With these notations, represents the mean anomaly, G is the projection of the angular momentum of the asteroid on the direction of the unit vector k and g is the anomaly of the perihelion P with respect to the direction of x 1 . Using the coordinates (1.7), condition (1.6) becomes ε ă 1 2 , ε :" βa r (1.8) as a body moving on an ellipse does not go further than twice the semi-axis from the focus of the ellipse. The canonical character of the coordinates (1.7) has been discussed, in a more general setting, in Pinzari (2019). In terms of the coordinates (1.7), the Hamiltonian (1.5) reads where, for short, we have let " pΛ, G, q " 1´e cos ξpλ, G, q , p " ppΛ, G, , gq " pcos ξ´eq cos g´G Λ sin ξ sin g. Here, e " epΛ, Gq " c 1´G 2 Λ 2 is the eccentricity, and ξ " ξpΛ, G, q denotes the eccentric anomaly, defined as the solution of Kepler's equation ξ´epΛ, Gq sin ξ " . The next step is to switch to the 2-DoF -averaged (hereafter, secular) Hamiltonian, which we write asH pG, R, g, rq " 1 2π 2Λ 2`σ KpR, r, Gq`σUpG, g, rq , (1.10) with KpR, r, Gq :" UpG, g, rq :" U`pG, g, rq`U´pG, g, rq`2 m 2 0 r (1.11) where U˘pG, g, rq :"´m In (1.10) we have omitted to write Λ and C among the arguments ofH, as they now play the rôle of parameters. Observe that the function U is π-periodic in g, as changing g with g`π corresponds to swap U`and U´, as one readily sees from (1.9)-(1.12). We do not provide rigorous bounds ensuring that the secular problem may be regarded as a good model for the full problem. Heuristically, we expect that this is true as soon as (1.8) is strengthened requiring, also, r " β 3{2 a.
(1.13) Indeed, extracting r from the denominators of the two latter functions in (1.9) and expanding the resulting functions in powers of βa r , one sees that the lowest order terms depending on have size m 2 0 σβa r 2 " m 2 0 β 3 a r 2 (recall that σ " β 2 ). So, such terms are negligible compared to the size m 2 0 a of the Keplerian term, provided that (1.13) is verified. Neglecting the constant term´m 5 0 2Λ 2 and, after a further change of time, the common factor σ in the remaining terms, the secular Hamiltonian (1.10) reduces tô HpG, R, g, rq " KpR, r, Gq`UpG, g, rq . (1.14) We now specify the range of parameters C, Λ and β and the region of the phase space for the coordinates pG, R, g, rq that we consider in this paper. In particular, we look for values of parameters and coordinates where the Hamiltonian (1.14) is weakly coupled, and describe the motions we expect to find in such region. As above, our discussion will be extremely informal. First of all, we take Λ and C verifying This condition implies that also |G| ! C (as |G| ă Λ) and hence K affords the natural splitting We consider a region of phase-space where r and R take values These are the values where K 0 attains its minimum, and correspond to circular motions of the planet, with r 0 being the radius of the circle. In the region of phase space defined by (1.16), the relative sizes of K 1 and U to K 0 are where c i are independent of m 0 , β, Λ and C. Even though (by (1.15) and (1.13)) K 1 and U are small compared to K 0 , however, they cannot be neglected, as their sum governs the slow motions of the coordinates G and g, which do not appear in K 0 . Remark that K 1 and U are coupled with K 0 , since they depend on r. It is however reasonable to expect that, as long as the minimum of K 0 cages r to be close to the value r 0 , the coupling is weak and the dynamics of G and g is, at a first approximation, governed by the 1 DoF Hamiltonian FpG, gq :" pK 1`U q| r"r0 .
( 1.18) To understand the global phase portrait of F in the plane pg, Gq, we need to recall some results from Pinzari (2020b). We go back to the functions U˘in (1.12), which enter in the definition of U. In (Pinzari, 2020b, Section 3), it is proved that, under the assumption (1.8), the following identity holds with ε as in (1.8) and t˘pG, g, εq :" The equality (1.19) has two main consequences. The former is that, even though the transformation (1.7) looses its meaning when G " 0, however, U˘keep their regularity, provided that (1.8) holds. Indeed, the functions t˘are regular at G " 0 and, being bounded below by´1 and above by 1, the denominator of the function under the integral never vanishes, under (1.8), as it is immediate to verify. Secondly, the phase portrait of the functions U`p¨,¨, rq, U´p¨,¨, rq coincides, a part for a rescaling, with the one of t`p¨,¨, εq, t´p¨,¨, εq, respectively. In particular, U`p¨,¨, rq and U´p¨,¨, rq have elliptic equilibria at pG, gq " p0, 0q and pG, gq " p0, πq, because this is true for t˘, as it is immediate to check. The phase portrait of t`p¨,¨, εq for ε ă 1 2 is shown in Figure 2 (left); the one of t´p¨,¨, εq is specular, interchanging the equilibria. We now merge these informations, in order to build the phase portrait of the function F in (1.18). By the Implicit Function Theorem, one can argue that, for an open set of values of the parameters, due to the linear term in G in K 1 , the equilibria of U`and Uá re shifted along the G-axis, but are not destroyed. Quantifying this shift is not easy, as U has an involved dependence on t`, t´. Based on the ε-expansion of U, with m 0 " 1 , C " 75 , Λ " ? a " 3 , β " 40 (1.20) (which comply with (1.8), (1.13), (1.15)) we obtain the phase portrait of F as in Figure 2 (right). We observe that, at contrast with the figure at left-hand side, where the motions are purely of elliptic kind, the phase portrait at right-hand side also includes rotational motions. The linear term of K 1 is responsible of this fact, breaking the symmetry G Ñ´G. We underline at this respect that the present framework is in a sense complementary to the one studied in Pinzari (2020b), where the phase portrait of F has, in fact, only elliptic motions: in that case, the linear term of K 1 does not exist, as C is fixed to 0. Remark also that the vanishing of C in Pinzari (2020b) affects condition (1.15) (which is not satisfied) and the motions generated by K 0 (which are collisional, rather than circular). The purpose of this paper is to show that, if the parameters are chosen about (1.20) and the energy is fixed to the level of a suitable initial datum pG ‹ , R ‹ , g ‹ , r ‹ q satisfying (1.16) (see Appendix B.1 for the exact values), then, in the system (1.14) a topological horseshoe wakes up in the plane pg, G{Λq. The analysis will be purely numerical, based on techniques developed in Gierzkiewicz and Zgliczyński (2019), Zgliczynski and Gidea (2004). More details on the methodological strategy are given along the following sections.

Poincaré mapping
From now on, we neglect to write the "hat" in (1.14). Moreover, for the purposes of the computation, we replace the function U with a finite sum where q ν pG, g, rq are the Taylor coefficients in the expansion of U with ν " 1, . . ., k. Using the parity of U as a function of r, these coefficients have the form In our numerical implementation, we use the truncation in (2.1) with k " k max " 10, so as to balance accuracy and number of produced terms. We still denote as H the resulting Hamiltonian: The study of the secular 2-DoF Hamiltonian in the continuous time t can be reduced to the study of a discrete mapping through the introduction of ad-hoc Poincaré's section Meiss (1992). The advantage consists in reducing further the dimensionality of the phase-space, and, in the case of n " 2, to sharpen the visualisation of the dynamical system. In fact, for a 2-DoF system, the phase-space has dimension 4 and, due to the conservation of the energy (the Hamiltonian H itself), orbits evolve on a three-dimensional manifold M . By choosing an appropriate surface Σ transverse to the flow, one can look at the intersections of the orbits on the intersection of M X Σ, i.e., a two-dimensional surface. The surface Σ chosen is a plan passing through a given point pG ‹ , g ‹ , r ‹ q and normal to the associated orbit, i.e., to the velocity vector pv ‹ Let us now formally introduce the Poincaré map. We start by defining two operators l and π consisting in "lifting" the initial two-dimensional seed z " pG, gq to the four-dimensional space pG, R, g, rq and "projecting" it back to plan after the action of the flow-map Φ t H during the first return time τ . The lift operator reconstructs the four-dimensional state vector from a seed on DˆT{2, where the domain D of the variable G is a compact subset of the form r´Λ, Λs. For a suitable pA, Aq Ă R 2ˆR2 , its definition reads wherez " pG, g, R, rq satisfies the two following conditions: (1) Planarity condition. The triplet pG, g, rq belongs to the plane Σ, i.e., r solves the algebraic The Hamiltonian is separable in R, so this condition amounts to solve a quadratic equation. If R 2 ě 0, then we choose the root associated to the "positive" branch`?R 2 . If R 2 ă 0, then we are led to the notion of inadmissible seed. The set of admissible seeds, noted by A, for the chosen section Σ is portrayed in Figure 3.
The Poincaré mapping is therefore defined and constructed as P : DˆT{2 Ñ DˆT{2 z Þ Ñ z 1 " P pzq "`π˝Φ τ pzq H˝l˘p zq . The mapping is nothing else than a "snapshots" of the whole flow at specific return time τ . It should be noted that the successive (first) return time is in general function of the current seed (initial condition or current state), i.e., τ " τ pzq, formally defined (if it exists) as where`Gptq, gptq, rptq˘is obtained though Φ t H pG, R, g, rq. The Poincaré return map we described has been constructed numerically based on the numerical integration of the Hamiltonian equation of motions (the details regarding our numerical settings are presented in the Appendix B.) This mapping being now explicit, we are able to unveil the phase-space structures through successive iterations of P . Figure 4 presents the successive coordinates of tP n pzqu where the initial seeds z cover a discretisation of DˆT{2 domain (mesh) and n " 10 3 . The phase-space structures can be roughly categorised in three distinct zones. In the lower part, say for G ă´2, we can distinguish one "pic" centred around g " π{2. One elliptic zone is immersed inside this structure, surrounded by "scattered dots", indicative of chaos. There is a large region of the phase-space foliated by circulational tori. The last upper region is a large zone where almost all regular structures have disappeared. The panel provided by Figure 4 presents some magnifications of phase-space structures. The obtained phase-space structures have been confirmed using a finite time dynamical chaos indicator, the Fast Lyapunov Indicator (FLI) computed with the whole flow on an iso-energetic section (see Appendix C for more details). The FLIs computation relies on monitoring the growth over time of the tangent vector under the action of the tangent flow-map (variational dynamics). The final FLIs values are colour coded according to their values and projected onto the section to provide a stability chart. Stable orbits correspond to dark regions, orbits possessing the sensitivity to initial conditions appear in reddish/yellow color. As shown in Figure 4, the FLIs confirm nicely the global structures depicted via the mapping. Moreover, numeric suggests that the lift of P on the variables pG, g, rq (i.e., the map obtained from Φ τ H by projection on pG, g, rq) is generically twist. 2.1. Hyperbolic structures and heteroclinic intersections. Equilibrium points of the mapping P (i.e., periodic orbits of the Hamiltonian system (2.2)), have been found using a Newton algorithm with initial guesses distributed on a resolved grid of initial conditions in DˆT{2 (again, see Appendix B for more details regarding the numerical setup). We found more than 20 fixed points x ‹ whose coordinates have been reported in Appendix B.4. The eigensystems associated to the fixed points have been computed to determine the local stability properties. The point x ‹ is hyperbolic when one of its real eigenvalues has modulus greater than one, the other less than one (expanding and contracting directions, respectively). In the case of complex eigenvalues, the point is elliptical. The result of the analysis is displayed on Figure 5 along with the following convention: hyperbolic fixed points appear as red crosses, elliptical points are marked with blue circles. As intuitively expected, the hyperbolic points are embedded within the chaotic sea. On the contrary, the stable islands host the elliptic points. Note that even the fixed-point in the small stability island has been recovered with the Newton scheme. In the vicinity of the unstable fixed-points, the dynamics is dominated by Figure 5. Phase-space of P together with its fixed-points. Hyperbolic points appear with red crosses, elliptical points appear with blue circles. the stable and unstable manifolds who have the eigenvectors of DP px ‹ q asymptotically tangents near x ‹ . The local stable manifold associated to an hyperbolic point x ‹ , can be grown by computing the images of a fundamental domain I Ă E s px ‹ q, E s px ‹ q being the stable eigenspace associated to the saddle point x ‹ . We considered the simplest parametrisation of I, namely a normalised version of the eigenvector associated to the saddle point x ‹ . This allowed us to compute a piece of W s loc. px ‹ q under the action of the flow-map Simó (1990); Krauskopf et al. (2006). To compute the unstable manifold, the same computations are performed by reversing the time and changing E s by E u . Finite pieces of those manifolds are presented in Figure 6 for two saddle points. Following the well established conventions of the cardiovascular system (as reported in Meiss (2008)), the stable manifolds are displayed with blue tones, unstable manifolds appear in red tones. As we can observe, those curves intersect transversally forming the sets of heteroclinic points, trademark of the heteroclinic tangle and chaos Morbidelli (2002). We now have at hands all the necessary ingredients and tools to prove the existence of symbolic dynamics using covering relationships and their images under P .

Symbolic dynamics via covering relations
In this section we prove the existence of symbolic dynamics for the considered model. The tools rely on ad-hoc covering relations that we present briefly following Gierzkiewicz and Zgliczyński (2019), in particular for the case n " 2.
3.1. Covering relations and topological horseshoe. Let us introduce some notations. Let N be a compact set contained in R 2 and upN q " spN q " 1 being, respectively, the exit and entry dimension (two real numbers such that their sum is equal to the dimension of the space containing N ); let c N : R 2 Ñ R 2 be an homeomorphism such that c N pN q " r´1, 1s 2 ; let N c " r´1, 1s 2 , Figure 6. Finite pieces of manifolds of two hyperbolic fixed points q 1 , q 2 . Their stable and unstable manifolds intersect transversally in (more than one) heteroclinic points. Stable manifolds are in blue while unstable manifolds are in red.
Nć " t´1, 1uˆr´1, 1s, Nc " r´1, 1sˆt´1, 1u; then, the two set N´" c´1 N pNć q and N`" c´1 N pNc q are, respectively, the exit set and the entry set. In the case of dimension 2, they are topologically a sum of two disjoint intervals. The quadruple pN, upN q, spN q, c N q is called a h-set and N is called support of the h-set. Finally, let SpN q l c " p´8,´1qˆR, SpN q r c " p1, 8qˆR, and SpN q l " c´1 N pSpN q l c q, SpN q r " c´1 N pSpN q r c q be, respectively, the left and the right side of N . The general definition of covering relation can be found in Gierzkiewicz and Zgliczyński (2019). Here we provide a simplified notion, suited to the case that N is two-dimensional, based on 2 (Zgliczynski and Gidea, 2004, Theorem 16).
Definition 3.1. Let f : R 2 Ñ R 2 be a continuous map and N and M the supports of two h-sets.
We say that M f -covers N and we denote it by M f ùñ N if: Conditions (2) and (3) are called, respectively, exit and entry condition.
The case of self-covering is not excluded. The Figure 7 shows two schematic examples of covering relation between two different sets N, M and a self-covering relation of N . The notions of covering relationships are useful in defining topological horseshoe (confer Gierzkiewicz and Zgliczyński (2019); 2 More precisely, Definition 3.1 is based on the proof of (Zgliczynski and Gidea, 2004, Theorem 16). Indeed, (Zgliczynski and Gidea, 2004, Theorem 16) asserts that under conditions (1), (3) and one of the inclusions in (Zgliczynski and Gidea, 2004, (78) or (79)), one has M f ùñ N in the sense of Gierzkiewicz and Zgliczyński (2019). However, during the proof of (Zgliczynski and Gidea, 2004, Theorem 16), inclusions (Zgliczynski and Gidea, 2004, (78) or (79)) are only used to check the validity of (2).   (2004)).

Zgliczynski and Gidea
Definition 3.2. Let N 1 and N 2 be the supports of two disjoint h-sets in R 2 . A continuous map f : R 2 Ñ R 2 is said to be a topological horseshoe for N 1 and N 2 if Topological horseshoes are associated to symbolic dynamics as presented in Theorem 2 in Gierzkiewicz and Zgliczyński (2019) and Theorem 18 in Zgliczynski and Gidea (2004), where the authors show that the existence of a horseshoe for a map f provides a semi-conjugacy between f and a shift map t0, 1u Z , meaning that for any sequence of symbols 0 and 1 there exists an orbit generated by f passing through the sets N 1 and N 2 in the order given by the sequence, guaranteeing the existence of "any kind of orbit" (periodic orbits, chaotic orbits, etc.).
From the Definition 3.1, the covering relation N 1 P ùñ N 2 is verified if the three following conditions are satisfied: (1) the image P pN 1 q of N 1 lies in the strip between the top and the bottom edges of N 2 , (2) the image of the left part of N1 lies on the left of N 2 , (3) the image of the right part of N1 lies on the right of N 2 ; the conditions can be easily checked in Figure 8 and, then, in Figure 9.
3.2. Existence of a topological horseshoe. In this section we describe how we construct explicitly a topological horseshoe for the Poincaré map of the Hamiltonian (1.10). We start by considering one hyperbolic fixed point q for the Poincaré map P and we denote by v s and v u , respectively, the stable and the unstable eigenvectors related to DP pqq. We construct a parallelogram N containing q whose edges are parallel to v s and v u and thus we define N as where A and B are suitable chosen closed real intervals. If the intervals A and B are sufficiently small, under the action of the map P , the parallelogram N will be contracted in the stable direction and expanded in the unstable direction. We denote by P pN q the image of N through the map P . In practice, we choose two hyperbolic fixed points q 1 and q 2 having the important property of transversal intersection of their stable and unstable manifolds as shown in Figure 6. This property is a good indication of the existence of a topological horseshoe. Based on this couple of fixed points whose coordinates read # q 1 " pg 1 , G 1 q " p0.203945459, 2.06302430q, q 2 " pg 2 , G 2 q " p0.278077917, 2.21418596q, we define two sets N 1 , N 2 Ă R 2 which are supports of two h-sets as follows: and v s 1 , v u 1 , v s 2 , v u 2 are the stable and the unstable eigenvectors related to q 1 , q 2 , respectively. Then the following covering relations hold proving the existence of a topological horseshoe for P , i.e., existence of symbolic dynamics for P . The obtained horseshoe associated to q 1 and q 2 with the aforementioned parameters is illustrated in Figure 9, providing the existence of symbolic dynamics.

Conclusions and open problems
This work originates from Pinzari (2019), where it has been pointed out that the average U`(1.12) of the Newtonian potential with respect to one of the two mean anomalies is an integrable function which in turn may be written as a function of another function t`, whose dynamics is completely known. The functional dependence (1.19) between these two functions, holding in the case of the planar problem, has been pointed out in (Pinzari, 2020b, Section 3). The identity (1.19) raises the very natural question whether and at which extent such relation has a consequence on the dynamics of the three-body problem. Giving an answer to such question is in fact demanding, as it requires to understand whether it is possible to find a region of phase space where the three-body Hamiltonian is well represented by its simple average (here "simple average" is used as opposite to "double average", most often encountered in the literature, e.g., Laskar and Robutel (1995)) and, simultaneously, the kinetic term K in (1.11) does not interfere with U too much. In Pinzari (2020b) it has been proved that if the total angular momentum C of the system vanishes, by symmetry reasons, and using a well-suited perturbation theory, the librational motions of t`reported in Figure 2 (left) have a continuation in the averaged three-body problem. In this paper we investigated the case C ‰ 0. With purely pioneering spirit, in order to simplify the analysis, we focused on the very peculiar situation where the two minor bodies have equal mass and we fixed an energy level once forever. We believe that both such choices can be removed without affecting the results too much, because, as informally discussed in Figure 9. Horseshoe connecting the points q 1 and q 2 proving symbolic dynamics for the map P . Red represents the entry sets and their images and blue the exit sets and their images.
the introduction, what really matters is the relative weight of K and U. Figures 2 and 4 not only show that, in our simplified model, this continuation is numerically evident, but also exhibit the onset of chaos in certain zones, clearly highlighted along the paper using techniques of Gierzkiewicz and Zgliczyński (2019). Even though the results are encouraging, many questions are still pending (some of them have been pointed out in Pinzari (2020b)), and we aim to face them in the future: pQ 1 q If C ‰ 0, is there a choice of parameters and phase space where the phase portrait of F includes only librational motions? pQ 2 q In the case that the orbit of the planet is inner to the one of the asteroids, the phase portrait of U`includes a saddle and a separatrix through it (see (Pinzari, 2020b, Figures 1, 2 and 3)). How does this affect the three-body problem motions? pQ 3 q By Pinzari (2019), relation (1.19) has a generalisation to the spatial problem. What are the consequences on the spatial three-body problem? pQ 4 q Is the onset of chaos in the averaged problem present also in the full (non-secular) system? pQ 5 q What can we prove analytically? pQ 6 q What can we prove with computer-assisted techniques?

B. Numerical setups and results
B.1. Choice of the parameters. The analysis we have done is related to the choice of parameters and initial data we started with. The Hamiltonian (1.10) is composed by three parts where the first one is the unperturbed and constant part depending on Λ, the second one represents the kinetic part and the third is the perturbing part. To ensure the non-resonant terms of P to be small with respect to H 0 we choose, as mentioned in the introduction, The initial datum is taken to be Note that R ‹ and r ‹ verify (1.16) but are not exactly centred at 0 and r 0 because the r-component of the Hamiltonian vector-field vanishes for R " 0, while it needs to be different from zero in order that the Poincaré map is well defined. The values of G ‹ and g ‹ have been empirically chosen such that the orbit from from pG ‹ , R ‹ , g ‹ , r ‹ q is approximately periodic and hence the Poincaré map is well defined.
B.2. Flow. The Hamiltonian equations of motion have been numerically propagated using a fixed time-step RK4 method Press et al. (1992). Even though the step has been kept fixed, no numerical issues have been encountered and the integration times were reasonable for the whole numeric exploration. Under the choice of our time-step δ, the flow-map preserves the Hamiltonian itself, a conserved quantity (first integral ), with a relative error of about 10´1 4 for stable orbits and 10´1 2 for chaotic orbits on a arc length of about τ " 10 2 orbital revolutions. Besides the first integral being numerically well preserved, the quality of the integration has been assessed further using a forwards/backwards strategy. The method consists in propagating forwards in time (say on r0, τ s) the Cauchy problem # 9 x " v H pxq, xp0q " x 0 , and then to back-propagate (from τ to 0) the new Cauchy problem # 9 x " v H pxq, where the initial seed x τ is obtained from the forward numerical flow-map, x τ " Φ τ px 0 q. Then the relative error is estimated. On a selection of orbits, we found ∆ to be of the order of 10´1 2 for regular orbits, 10´8 for chaotic orbits on timescale of about 10 2 orbital revolutions.
B.3. Poincaré mapping P . The construction of the Poincaré map P is based on the time evolution of the whole flow and a bisection procedure. Given an initial point z, to find its next state z 1 " P pzq we compute xptq " Φ t px 0 q, x " lpzq, until following conditions are met: (1) Section condition: X " px 1 ptq, x 2 ptq, x 4 ptqq P Σ up to a numerical tolerance ε Σ " 10´1 0 . This step relies on a bisection method halving the length of the numerical step δ until we drop under the tolerance ε Σ . (2) Orientation condition: The scalar product 9 Xp0q¨9 Xptq is positive, meaning that the orbit is intersecting the plan Σ in the same direction as the starting point.

C. The Fast Lyapunov Indicator & dynamical timescales
The Fast Lyapunov Indicator (FLI) is an easily implementable tool suited to detect phase-space structures and local divergence of nearby orbits. It has a long-lasting tradition with problems motivated by Celestial Mechanics Froeschlé et al. (1997). The indicator can be used in the context of deterministic ODEs, mappings, and is able to detect manifolds and global phase-space structures Froeschlé et al. (2000); Guzzo and Lega (2014); Lega et al. (2016). A large literature exists with the FLI tested on idealised systems (e.g., low dimensional quasi-integrable Hamiltonian system Froeschlé et al. (2000); Guzzo and Lega (2013), drift in volume-preserving mappings Guillery and Meiss (2017)) but also on many applied gravitational problems across a variety of scales, ranging from the near-Earth space environment Daquin et al. (2018) to exoplanetary systems Páez and Efthymiopoulos (2015). For simplicity, let us present the tool in the case of ODEs. Let us assume we are dealing with a n-dimensional autonomous ODE system. If the system is non-autonomous, we classically extend the dimension of the phase-space by 1 dimension. The FLI indicator is based on the variational dynamics in R 2n , # 9 x " f pxq, 9 w " Df pxq¨w, w P T x M , and is defined at time t as FLIpx 0 , w 0 , tq " sup 0ďτ ďt log wpτ q . (C.1) The FLI is able to distinguish quickly the nature of the orbit emanating from x 0 . Orbits containing the germ of hyperbolicity will have their final FLI values larger than regular orbit (for the same horizon time τ ). More precisely, chaotic orbits will display a linear growth (with respect to time) of their FLIs, whilst regular orbits have their FLIs growing logarithmically. In order to reduce the parametric dependence of the FLIs upon the choice of the initial tangent vector, the FLIs are computed over an orthonormal basis of the tangent space (i.e., we compute Eq. (C.1) 4 times with a different initial w 0 ) and averaged Guzzo and Lega (2013). As a rule of thumb, the FLI is computed over a few Lyapunov Figure 10. On the left: calibration of the finite time chaos indicators (FLI) for three distinct orbits. After a transient time of about t " 5, 000 (i.e., " 10 orbital revolutions) the discrimination of the nature of the orbits is sharp enough. On the right: time evolution of the maximal Lyapunov exponents χ. For chaotic orbits, they define Lyapunov timescales of about 0.76 orbital revolutions.
times τ L , but in practice, the choice of the final τ requires a calibration procedure by testing few orbits. By computing FLIs on discretised domains of initial conditions, the color coding of the FLIs (using a divergent color palette) reveals the global topology of the phase-space (e.g., web of resonances and preferred routes of transport, see Guzzo and Lega (2013)) furnishing a so-called stability map. The Lyapunov time τ L is obtained as the inverse of the maximal Lyapunov characteristic exponent (we refer to Skokos (2010) for computational aspects related to characteristic exponents), where χ denotes the maximal Lyapunov characteristic exponent χpx 0 , w 0 q " lim tÑ`8 1 t log wptq .
Stable orbits do satisfy χ Ñ 0 and hence τ L tends to be large. On the contrary, chaotic orbits are characterised by χ Ñ r P R ‹ and therefore τ L converges to a finite value. The panel shown in Figure 10 presents the calibration procedure based on three orbits. The stable orbit displayed in black (logarithmic growth of the FLI) admits for initial condition pG, gq " p´2, πq. The two others orbits are chaotic but one (red) is less hyperbolic than the other (blue). The respective initial conditions read pG, gq " p´2, 1.6q and pG, gq " p2, πq. As it is observed, after a transient time of about t " 5, 000 (i.e., 10 orbital revolutions), safe conclusions can be formulated regarding the stability of the orbits (left panel). The respective maximal Lyapunov characteristic exponents are presented in the right panel of Figure 10. The inverse, the Lyapunov time, defines timescales of " 380 for the most chaotic one (which is about 0.76 revolutions) and " 1, 100 for the second chaotic one (2.2 orbital revolutions).