Numerical studies to detect chaotic motion in the full planar averaged three-body problem

In this paper, the author deals with a well-known problem of Celestial Mechanics, namely the three-body problem. A numerical analysis has been done in order to prove existence of chaotic motions of the full-averaged problem in particular conﬁgurations. Full because all the three bodies have non-negligible masses and averaged because the Hamiltonian describing the system has been averaged with respect to a fast angle. A reduction of degrees of freedom and of the phase-space is performed in order to apply the notion of covering relations and symbolic dynamics


Introduction
The three-body problem is one of the most studied issues in Celestial Mechanics. It deals with the motion of three bodies interacting by only the gravitational force. The problem is, generally, non-integrable, namely, no analytical closed solutions exist to describe the dynamics of the model. Many analytical and numerical studies have been done in order to acquire information about the motions of the associated dynamical system, proving existence of both stable, unstable, regular and chaotic dynamics (see, for example, [3] for a complete treatment of the problem).
The purpose of the current paper is to collect and to link two recent works (we refer to papers [8,9]) where the onset of chaos is numerically proved in two different configurations of the three-body problem. In both the works, we aim to write the Hamiltonian describing the model as the sum of a Keplerian part ruling the motion of two of three bodies plus a part depending on all the three bodies. In principle, we do not have dominant interactions between the three bodies but we choose a region of the phase-space, where the remaining part of the Hamiltonian becomes a perturbation, so we deal with a "perturbed Keplerian problem". It can be possible by introducing new coordinates associated to the three bodies (see [15] for a B Sara Di Ruzza sara.diruzza@unipa.it 1 DMI, University of Palermo, Via Archirafi 34, 90123 Palermo, Italy detailed study of the new coordinates) and through suitable rescalings of variables and time. In works [8,9], two different configurations and phase-spaces are considered and in both cases we use symbolic dynamics to prove, numerically, the existence of chaotic motions. In Sect. 1.1, we aim to write the Hamiltonian of the two different configurations in the same way, in order to have a generic treatment of the problem. Indeed, the goal of the current paper is to show a method to apply to the perturbed Keplerian problem when some conditions are satisfied. Afterwards, the two different cases of papers [8,9] are described.

The general setting
Let us consider 3 bodies P 0 , P 1 , P 2 with masses m 0 , m 1 , m 2 interacting by only the gravitational force. By fixing an orthonormal reference frame (i, j, k) in the Euclidean space identified with R 3 , the Hamiltonian H 3b describing the system can be written as where x 0 , x 1 , x 2 ∈ R 3 and y 0 , y 1 , y 2 ∈ R 3 are, respectively, the positions and the impulses of the three bodies; | · | stands for the Euclidean distance and the gravity constant has been fixed to one. At the beginning, we do not impose any restrictions on the masses of the bodies and for this reason we can refer to the Hamiltonian (1) as full three-body problem, where full is opposed, in this case, to restricted where one of the three bodies has negligible mass and it does not affect the motion of the other two bodies. Moreover, in our study, we consider the planar problem, i.e. the motion is restricted to a plane by choosing positions and impulses x 0 , x 1 , x 2 , y 0 , y 1 , y 2 ∈ R 2 . After a suitable change of the Cartesian coordinates and a rescaling of coordinates and Hamiltonian itself, we aim to write the Hamiltonian H 3b in the following way: H 3b (y , y, x , x) = |y| 2 2 − 1 |x| +Ĥ (y , y, x , x; η) (2) where (x, y), (x , y ) ∈ R 2 × R 2 are Cartesian positions and impulses of two of the bodies whose origin is set in a suitable point 1 and η = η(m 0 , m 1 , m 2 ) is a mass parameter. The first two terms |y| 2 2 − 1 |x| are called Keplerian part and, assuming it takes negative values, it generates motions for the body of coordinates (x, y) on an ellipse. The termĤ is the remaining part of the Hamiltonian. In principle, we do not ask thatĤ is necessarily a perturbation since the coordinates are not centered at the most massive body, although we will work in regions of the phase-space such thatĤ is small. An analysis of this point will be addressed in the following sections.
The Hamiltonian (2) has four degrees of freedom and in the following steps we aim to reduce it to a two degrees of freedom system. First of all, the Hamiltonian H 3b is invariant under the group of transformations SO(2), namely, orthogonal rotations. Introducing a set of canonical coordinates in a suitable rotating system we reduce this symmetry obtaining a three degrees of freedom system. To this aim, we introduce the total angular momentum, which is a constant of motion, as The other body or the barycenter, for example.
where k = i × j is the unit vector orthogonal to the plane where motions take place, and then, we define a six-dimensional "rotation-reduced phase space" introducing 6 new canonical coordinates ( , R, G, , r , g).
To define this set of coordinates we denote as "first body" the body P 1 of Cartesian coordinates (x, y) and as "second body" the body P 2 of Cartesian coordinates (x , y ). We denote by E the ellipse generated by the Keplerian term of Eq. (2) for given values of (x, y). Finally, we define two pairs of Delaunay variables ( , G, , g) ∈ R 2 × T 2 for the first body P 1 as follows: where a is the semi-major axis of E, -G = x × y · k is the angular momentum of the first body, -is the mean anomaly of x from the pericenter P of E, g is the angle detecting the pericenter of E with respect the direction of x , and a radial-polar pair of coordinates (R, r ) ∈ R 2 for the second body P 2 as namely, R is the radial velocity of x and r is the Euclidean length of x . We underline that the coordinates (R, G, , r , g, ) refer to a reference frame which is not inertial with respect to the "second body" and they have been introduced to reduce the symmetry of rotations (see [15] for a complete and detailed discussion about these variables). The Hamiltonian (2) written in the new coordinates appears as where the first term −1/(2 2 ) is the Keplerian term and the total angular momentum C is considered as parameter. The next step consists in averaging Eq. (3) with respect to the fast angle, the mean anomaly , in order to reduce again the degrees of freedom from three to two. The new averaged Hamiltonian takes the form: The three terms in parentheses-in which, for the moment, we omit the dependence on the mass parameter η-will be named as kinetic term, disturbing term and Newtonian potential term, respectively. Under two explicit assumptions about the terms appearing in Hamiltonian (4), namely, and we claim that, at first approximation, the motion of ( (t), R(t), G(t), (t), r (t), g(t)) under the Hamiltonian flow (4) is described by the following three rules: (a) (t) remains almost constant and (t) ≈ t moves fast; (b) the motion (R(t), r (t)) is ruled by K C ; (c) the motion (G(t), g(t)) is ruled by the non-autonomous Hamiltonian U (r (t), ·, ·).
Precise computations and assumptions on each term will be addressed in the following sections where two explicit examples will be numerically analyzed.

Euler integral function
Let us recall some important properties about the function U (as it has been proved in [15]). Its general formula is given by: and we remark that it is integrable; in particular, there exists a function F of two arguments, such that where The function E(r , G, g) is called Euler integral, because it appears in the integration of the two-fixed centers problem. By Eq. (7), the level sets of E are also level sets of U and the phase portrait of E can be easily studied (see [17]): in fact, the level sets of E are given by the curves where r is fixed and E takes different values; by periodicity of g, we consider as domain for the coordinates (g, G), the rectangle [0, 2π) × (−1, 1). Then, according to the fixed value of the variable r , three different phase portraits appear: a first case if 0 < r < 1, the second case if 1 < r < 2 and the third if r > 2. Let us recall that, the variable r is the norm of the vector x , thus, the three different cases stand for different configurations of the third body P 2 with respect to the first two bodies P 0 and P 1 , whose motion is ruled by the Keplerian part of the Hamiltonian. Let us now show in detail the phase space of the Euler integral function. In Fig. 1, we can distinguish three completely different phase portraits. In panel (a) the case 0 < r < 1 is represented. The point (0, 0) is a minimum, then there are two symmetric (with respect to G = 0) maxima in π, ± 1 − r 2 4 and finally the point (π, 0) is a saddle. Two separatrices split the phase space in different regions: the first separatrix is the curve S r (r ) := (G, g) : G 2 − r √ 1 − G 2 cos g = r passing through the saddle point (π, 0) and the second is the curve S 1 (r ) := (G, g) : G 2 − r √ 1 − G 2 cos g = 1 passing through the points ( π 2 , ±1). The curve S r (r ) delimits a region of librational motions around the minimum (0, 0); between S r (r ) and S 1 (r ) a region of rotational motions exists. Finally, S 1 (r ) delimits a region of librational motions around the two maxima. In panel (b) the case 1 < r < 2 is shown. The minimum (0, 0) persists as well as the two symmetric maxima, the saddle point and the two separatrices S r (r ) and S 1 (r ). Rotational motions disappear, S r (r ) delimits libration areas around the two maxima, other librational motions surround the separatrix S r (r ) and thus, the maxima and the saddle point, up to the separatrix S 1 (r ) which Fig. 1 Contour lines of the Euler integral function for fixed r : a 0 < r < 1; b 1 < r < 2; c r > 2 delimits the area with librational motion around the minimum. The structure completely changes in panel (c) representing the case r > 2. The saddle point and the separatrix S r (r ) disappear; (π, 0) turns to be the only maximum and (0, 0) remains a minimum. The separatrix S 1 (r ) still exists and it delimits two regions of librational motions around the minimum and the maximum, respectively.
The analysis of the dynamics of the three-body problem ruled by Hamiltonian (4) has been deeply developed in the last years. In particular, from an analytical point of view it has been studied by using normal form theory in order to prove stability estimates. For a complete treatment of the subject we refer the reader to the articles [4,6,16,17]. At the same time, the problem has been studied from a numerical point of view in order to prove existence of chaotic motions under particular initial configurations. We refer, in particular to papers [8,9].
The two numerical studies will be presented and compared in this paper as follow: in Sect. 2 we present two different configurations and write how the respective Hamiltonian functions are reduced; in Sect. 3, we explicitly define the Poincaré maps and study the hyperbolic structures of the two models; in Sect. 4 we provide some definitions and properties of covering relations and symbolic dynamics; in Sect. 5, we prove, numerically, the existence of symbolic dynamics and the onset of chaos in the two studied models. Finally, in Sect. 6, we provide conclusions and ongoing works.

Two different configurations
We want to present two different cases related to two different configurations of the threebody problem where chaos appears. The two cases have been studied, separately, in papers [9] and [8], respectively.

First case: 0 < r < 1
Let us show the first model. We fix three bodies P 0 , P 1 , P 2 with masses, respectively, m 0 = μ, m 1 = κ, m 2 = 1, with μ, κ < 1. Through a translation of coordinates, we center the reference system in the body P 0 with mass μ and write the Hamiltonian (1) where x and x are, respectively, the positions of bodies P 1 and P 2 with masses κ and 1. The Hamiltonian takes the form: We underline that this Hamiltonian is not centered in the more massive body: the reference frame has the origin in the body P 0 with mass μ, thus, the Hamiltonian describes the motions of the bodies P 1 , P 2 with masses κ, 1, respectively. Figure 2 shows the model. A suitable rescaling of variables (y , y, x , x) and of Hamiltonian through a rescaled time provides the new Hamiltonian in the form given by Eq. (2) as with From now on, in our computations, we choose α and β as independent mass parameters. Hamiltonian (8) has four degrees of freedom and, as it has been discussed in Sect. 1.1, we reduce one degree of freedom by exploiting the rotational invariance and using the conservation of the total angular momentum C. We introduce the Delaunay variables ( , G, , g) relative to the Keplerian motion of P 1 with respect to P 0 , and the couple (R, r ) of radial coordinates relative to the body P 2 . In the new coordinates, Hamiltonian (8) can be written as Eq. (3), with η = η(α, β). To reduce one more degree of freedom, we average the Hamiltonian with respect the fast angle , obtaining Eq. (4). Without loss of generality, we choose = 1 (which implies that the semimajor axis a is equal to 1) and consider the Keplerian term −1/(2 2 ) as an additive constant which does not change the equations of motion; thus, we will study the motion ruled by the averaged Hamiltonian: with e = √ 1 − G 2 being the eccentricity, ν the true anomaly and ξ the eccentric anomaly related with the mean anomaly by means of the Kepler equation Let us recall some properties of the function U (r , G, g) (the reader can refer to papers [4,9,16,17] for more details): -the function U is singular if and only if 0 < r < 2 and (g, G) ∈ S 1 (r ); -the rate of divergence of U is logarithmic with respect to the distance from S 1 (r ).
In this section we are dealing with the case 0 < r < 1 and, moreover, we require that r be less than the radius of the Keplerian ellipse, meaning that P 2 is moving inside the Keplerian ellipse which is defined by the motion of P 1 (consider again Fig. 2). Being r small, the potential U can be expanded in power series of r as where N is a suitable order of truncation.
Let us now choose values of parameters 2 fulfilling assumptions (5) and (6): we fix an energy level H fix and, on such level, we find initial conditions e : such that the orbit starting at (R e , G e , r e , g e ) and evolving under the Hamiltonian H is periodic in the four-dimensional phase space (R, G, r , g). The projection of such orbit on the three-dimensional phase space (G, r , g) is still periodic. Now, we are ready to study the hyperbolic structure of this problem by introducing a suitable Poincaré map. Before continuing, we introduce a new model and then, we proceed with the definition of first return maps on Sect. 3.

Second case: r > 2
Let us now present a new model: recalling again Hamiltonian (1), we choose m 1 = μm 0 and m 2 = κm 0 with μ = 1 κ so that P 0 , P 1 are bodies with equal mass and much smaller than the body P 2 . The Hamiltonian (1) is translational invariant, so we rapidly switch to a translation-free Hamiltonian by applying the well-known Jacobi reduction (see, for example, [12]). We recall that this reduction consists of using, as position coordinates, the centre of mass of the system (which moves linearly in time), the relative position vector x of two of the three bodies and the position vector x of the third body with respect to the centre of mass of the former two. Without loss of generality, we choose m 0 = 1 and rescale the coordinates and Hamiltonian as described in [8] to transform Hamiltonian (2) in the following explicit expression 3 : The choice κ μ = 1 gives β =β From now on, we regard β as mass parameter, with β ∼ κ and σ ∼ β 2 . By choosing a region of the phase-space where we ensure the denominators of the last two terms in (11) to be different from zero. All the assumptions made until now lead to a configuration where two bodies P 0 , P 1 with the same mass (for example two asteroids) move around their barycenter b and a third much more massive body P 2 (a planet) moves far from the first two. A schematic representation can be seen in Fig. 3.
As it has been described in Sect. 1.1, we introduce the couples of Delaunay variables ( , G, , g) for the asteroid relatively to x (the angle g is the angle of pericenter between the line passing through b, P and x direction) and the radial-polar coordinates (R, r ) for the planet. More precisely, the variables are defined as follows: where, considering the instantaneous ellipse generated by the first two terms in Hamiltonian (11), a is the semimajor axis, S and S tot are, respectively, the area of the ellipse spanned from the pericenter P and the total area. With this notation, G is the projection of the angular momentum of P 0 , P 1 on the direction orthogonal to the plane where motions take place, is the mean anomaly of P 1 , g is the argument of pericenter, R is the radial velocity of P 2 and r is its distance from the center of mass of P 0 and P 1 .
In the new coordinates the Hamiltonian reads: where is the eccentricity; ξ = ξ( , G, ) denotes the eccentric anomaly, defined as the solution of Kepler's equation (9). In the new coordinates, the assumption (12) takes the form: Thus, the last two terms of Eq. (13) can be expanded in powers of the small parameter β a/r ; then, the Hamiltonian can be averaged with respect to the fast angle in order to reduce the system from three to two degrees of freedom and to study the secular (averaged) problem. Let us write Equation (4) for this case: with where After expanding U ± in power series of β a/r and averaging with respect to , we replace the series of the Newtonian potential with a finite sum given explicitly by: where q ν (G, g, r ) 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 truncate Eq. (16) at k = k max = 10, so as to balance accuracy and number of produced terms. In this work we do not provide rigorous estimates to prove that the averaged problem is a good model for the full problem, but we fulfill some requirements which guarantee that the averaged problem well describes the full one. In a work in progress, we are better comparing the averaged and non-averaged problems. Neglecting the Keplerian term in Eq. (14) (being an additive constant) and rescaling time to get the parameter σ = 1, Hamiltonian (14) becomes with K C , K C , U as in (15). From now on, we will deal with the secular Hamiltonian (18). The dependence on β, , C does not appear in an explicit way being they considered as parameters.
Let us now show how the initial conditions for (G, R, g, r ) and the parameters C, , β are chosen. The aim is to find a region of the phase space and a range of parameters such that the terms in Eq. (18) become weakly coupled. The choice β 1 has been already made, as well as β a r < 1 2 . To ensure that the averaged problem well describes the full one, we also, heuristically, require that r β 3/2 a (see [8] for further details). Then, we take and C such that C, which implies also |G| C (being |G| < ). Moreover, we consider a region of the phase space where r = r 0 ∼ C 2 /2 and R ∼ 0, which allow the term K C to reach its minimum, so that the coupling between the term K C and the two terms K C , U is weaker and at a first approximation the motion of (G, g) is ruled by ( K C + U )| r =r 0 . In Fig. 4, we plot the phase portrait of the function ( K C + U )| r =r 0 on the plane (g, G), with U = U (r , G, g) truncated up to the order 1/r 3 and parameters as in Eq. (19). Such plot approximates the motion of the variables (g, G) replacing panel (c) of Fig. 1.
According to assumptions (5) and (6), we choose the following values for the parameters 4 : and the following initial conditions: We do not choose r = r 0 and R = 0 but values close to them, in order to well define, in the following, a section plane transversal to the orbit (as it will be described in Sect. 3). Moreover, G and g are chosen such that the orbit starting from (Ĝ,R,ĝ,r ) under the flow of H is periodic in the four-dimensional phase space (G, R, g, r ).

Poincaré map and hyperbolic structure
In this section we define a first return map in order to reduce the dimension of the phase space of the models introduced in Sect. 2. The simple idea is to start from the four-dimensional phase space (R, G, r , g), fix an energy level H fix to derive R as a function of (H fix , G, r , g) and consider the three-dimensional space (g, G, r ). In such a space, we consider the projection of propagated orbits under the flow of H , which are unidimensional curves in the threedimensional space (g, G, r ). Then, we fix a two-dimensional plane section orthogonal to a given orbit in a given point and study the first return map of orbits on such a plane. We introduce a generalized first return map as follows: consider a compact domain D ⊂ R 2 ; fix a two-dimensional plane in the three-dimensional space (g, G, r ) and define a two-dimensional map P as where z = (z 1 , z 2 ) ∈ D is an unordered couple of combinations of variables (g, G, r ), G), (g, r ), (G, r )}; we call z 3 the remaining variable such that (z 1 , z 2 , z 3 ) ∈ ; then we consider the initial condition (z 1 , z 2 , z 3 , R(H fix , G, r , g)) and propagate it under the flow t H obtaining the orbit t namely, τ is the time of the first return on the plane .
The aim is to find a suitable Poincaré map of the flow t H and look for fixed points (corresponding to periodic orbits of the flow). If some hyperbolic fixed points exist, it is possible to construct stable and unstable manifolds of such points finding homoclinic or heteroclinic intersections which are the base to apply symbolic dynamics as described in Sect. 4.
In the following subsections we define suitable Poincaré maps for the two cases described in Sect. 2 (studied, respectively, in papers [9] and [8]).

A case of homoclinic intersection
Let us continue introducing first return maps for the model described in Sect. 2.1. The initial condition (10), corresponds to a periodic orbit in the four dimensional space. We consider the projection of such orbit in the three-dimensional space (g, G, r ) and define a plane e , in such a space, orthogonal to the orbit at the point (g e , G e , r e ). We construct the two-dimensional map: where (g , G ) is the first return value on the plane e . We take a grid of initial conditions (g, G) on the domain D = [0; 2 π] × [−1; 1] and complete it so that the initial point (g, G, r ) ∈ e and (R, G, r , g) has a fixed energy level; then, we iterate the first return map many times (in the order of 10 3 ) to get a Poincaré section as shown in Fig. 7; we implement a Newton method in order to find fixed points of the map P H , e . We found two fixed points: (g e , G e ) and (g h , G h ). Studying the linear part of P H , e and computing its eigenvalues, we see that the former is an elliptic point and the latter is a hyperbolic point (having real eigenvalues, one greater and one smaller than 1). The two points correspond to periodic orbits in the three-dimensional space (g, G, r ). The first one, e , called elliptic, has a rotational behavior in the variable g while in the second, h , called hyperbolic, the variable g librates around the value π (see the propagated variables in Figs. 5 and 6, respectively). In Fig. 7, we provide a visualization of the Poincaré section: on the upper panel, the two dimensional map is shown, and the elliptic and hyperbolic fixed points are represented in blue and red, respectively: we can see elliptic islands surrounding the elliptic point, rotational curves corresponding to regular tori and some chaos in the bottom part of the plot; on the bottom panel, we represent a three-dimensional visualization where the plane section e and the two periodic orbits are shown (the elliptic one in blue and the hyperbolic one in red).
We want to focus on the hyperbolic orbit to prove, numerically, existence of chaos. Let us recall Fig. 1a: level curves of the Euler integral function are shown for a fixed value of the variable r ∈ (0, 1). We define two-dimensional manifolds in the three-dimensional space (r , G, g) as: which, as the E changes, are essentially surfaces where the Euler integral E(r , G, g) remains equal to the constant value E: We are also interested in the manifold corresponding to the separatrix S r (r ), namely, M(r ) = M r = {(r , G, g) : E(r , G, g) = r }: We analyze the relation between the hyperbolic orbit and the manifold M r to show the presence of chaos close to this manifold as we shown in Fig. 8: let us start by comparing the variation of the variable r and the variation of the Euler integral versus time along the hyperbolic orbit: the value of the Euler integral is always greater than the value of the variable r , and in particular the minimum value E min of the Euler integral is greater than the maximum This implies that the hyperbolic orbit touches the surface M(E min ) (in particular it lies on the surface in g = π being tangent on that point) and never touches the surface M r being the orbit always above this surface. In Fig. 9, it is possible to see a three-dimensional visualization of the hyperbolic orbit and of the two surfaces M(E min ) and M r .
Observing the orbit numerically, we see that the maximum value r max = 0.274496 and the minimum value r min = 0.022 of the variable r , are both reached when g = π. The section plane e intersects the orbit h in an intermediate point being the associated quadruplet (R, G, r , g) given by: We study, initially, the stable and unstable 5 manifolds of the hyperbolic fixed point (g h , G h ) of the Poincaré map P H , e and we claim (see Numerical Evidence 5.1 in [9]) that a Transverse 5 We recall that local stable and unstable manifolds associated to the hyperbolic fixed point x * of a map P are defined as Homoclinic Intersection between the stable and unstable manifolds for P H , e exists. 6 In fact, the two manifolds have a transverse intersection at the point (g h , G h ). Then, we consider different planes i h orthogonal to h at different points (g i , G i , r i ) on the curve. We consider first return maps on i h . We denote as * the plane orthogonal to h at (g * , G * , r * ) = (1.27 π, 0.346, 0.270) .
This point of h has been chosen for being "close" to M r . With this choice, we detect a homoclinic tangency 7 and absence of splitting 8 and claim (see Numerical Evidence 5.2 in [9]) that as soon as r i is chosen closer and closer to r max the stable tori zone becomes smaller and smaller. For i h = * , stable motions are not numerically detected, and the unstable, stable manifolds have a homoclinic tangency at (g * , G * ). In other words a splitting of such is depicted in purple; the elliptic and hyperbolic fixed points are represented with blue and red dots, respectively. Bottom: spatial visualization of the section plane e with the section map in purple, the periodic orbits e , h in blue and in red, respectively (figures already used in paper [9]) (color figure online) manifolds (which have the shape of S 0 (r * )) is not numerically detected. The idea is that the "near-homoclinic tangency" in (g * , G * ) becomes a "homoclinic tangency" in (π, G * * ) (where G * * is such that the point (π, G * * , r max ) ∈ h ). We consider sections up to the point (22) and not up to (π, G * * , r max ) because this point is a cuspid for h so that it is not possible to construct an orthogonal plane there. In Fig. 10, four first return maps (21) are shown where we represent the following objects: -the elliptic (dark-blue dot) and the hyperbolic (dark-red dot) fixed points; -rotational tori (purple curves); -chaotic motions (dotted purple); -the transverse/near-tangent homoclinic intersection between the stable (blue curve) and unstable (red curve) manifolds from (g i , G i ). The last step we propose in this section in order to prove numerically the detection of chaos is to define a new first return map on the curve h at the point (g * , G * , r * ). The new section plane is taken to be vertical and in particular we consider the plane * = {g = g * }.
The two-dimensional first return map is depicted in Fig. 11. The aspect of the stable and unstable manifolds changes drastically, but homoclinic intersections are present.

A case of heteroclinic intersection
Let us consider now the model described in Sect. 2.2 and introduce the first return map as it has been done in the previous section. The initial condition (20) corresponds to a periodic orbit in the four dimensional space. We consider the projection of such orbit in the threedimensional space (g, G, r ) and define a plane , in such a space, orthogonal to the orbit at the point (ĝ,Ĝ,r ). We construct a two-dimensional map P H , as: where (g , G ) is the first return value on the plane . Note that the domain of the variable g is, in this case, [0; π] instead of [0; 2 π] due to the π-periodicity of Eq. (17). We take a grid of initial conditions (g, G) on the domain D = [0; π] × [−1; 1] and complete the quadruplet of the initial conditions such that (g, G, r ) ∈ and (R, G, r , g) has a fixed energy level; then, we iterate the first return map many times (in the order of 10 3 ) to get a Poincaré section and we implement a Newton method in order to find fixed points of the map P H , . This case is completely different from the previous one because, here, many hyperbolic and many elliptic fixed points are found. In Fig. 12 the following objects can be recognized: -elliptic fixed points as blue dots; -hyperbolic fixed points as red crosses; Fig. 11 In the upper panel, the first return map on the plane * = {g = g * } is plotted; the section map is represented in purple, and the stable and unstable manifolds in blue and in red, respectively. In the center panel a three-dimensional visualizations on the space (g, G, r ) is presented, while in the bottom panel we show the projection on the (g, r )-plane. In both panels the elliptic and hyperbolic fixed points are represented, respectively, with dark-blue and dark-red dots; the elliptic and hyperbolic orbits, respectively, with dark-blue and dark-red curves and the stable and unstable manifolds appear in blue and red, respectively. Note that, in bottom panel the stable and unstable manifolds appear as a segment because they belong to a vertical plane parallel to the plane (G, r ) (figures already used in paper [9]) (color figure online) -rotational tori in dark; -librational islands in dark; -chaotic motions as dotted dark.
We computed stable and unstable manifolds of some hyperbolic fixed points. In Fig. 13 left, we consider one fixed point and show finite pieces of stable (in blue) and unstable (in red) manifolds of such a point. In this case, many homoclinic intersections can be seen. In Fig. 12 First return map on the plane . In black the section map is represented; elliptic and hyperbolic fixed points are depicted in blue and red, respectively (figures already used in paper [8]) (color figure online) Fig. 13 Finite pieces of stable and unstable manifolds are depicted in blue and red, respectively. Left: homoclinic intersections between stable and unstable manifolds related to a hyperbolic fixed point. Right: heteroclinic intersections between stable and unstable manifolds of two different hyperbolic fixed points (figures already used in paper [8]) (color figure online) Fig. 13 right, we choose two hyperbolic fixed points q 1 and q 2 and show a small region of their stable (in blue) and unstable (in red) manifolds. In such figure the heteroclinic transversal intersections between the manifolds of the two points are detected and well visible.

Symbolic dynamics
In this section we recall some important definitions of covering relations and we remark the link with symbolic dynamics and existence of chaotic motions. It has been done following papers [11,18,19], and as in [8,9], here we give a simplified version of the subject: in particular the results are valid in R n , while we restrict them to R 2 , being our first return maps defined on such a space.
Definition 1 (h-sets, [11,19]) Let N ⊂ R 2 be a compact set and let The sets are called, respectively, the exit set and the entry set; the sets are called, respectively, the left edge and the right edge of N ; the sets are called, respectively, the left side and the right side of N .
The following definition is fitted to the special case (realised in our study) that the unstable manifold has dimension 1. The simplification compared to the general definition in [11,19] is based on [19,Theorem 16].
In Fig. 14, we represent a schematic example of covering and self-covering relation.
Let us now proceed by defining a topological horseshoe as follows: Topological horseshoes are associated to symbolic dynamics as presented in Theorem 2 in [11] and Theorem 18 in [19], where the authors show that the existence of a horseshoe for a map f provides a semi-conjugacy between f and a shift map {0, 1} 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.). Let us also generalize the notion of symbolic dynamics, giving a weaker definition: we say that f has m-symbolic dynamics if there exist compact subsets with non-empty interior N 0 , N 1 ⊂ D such that, for every n ∈ N and any finite sequence (σ 0 , . . ., σ n ) of symbols σ i ∈ {0, 1} having length n + 1, one can find x 0 ∈ N σ 0 such that the orbit of x 0 under f , namely, x j := f j (x 0 ) is well defined for j = 0, . . ., nm, and x m j ∈ N σ j ∀ j = 0, . . . , n. 1-SYMBOLIC dynamics in N 0 ∪ N 1 is also called horseshoe.
In Fig. 15, we represent a schematic example of a 1-symbolic dynamics and a 3-symbolic dynamics in N 0 ∪ N 1 .
Let us continue by provide the following important results: be a continuous map such that Then, there exists x 0 ∈ N 0 such that We shall use Theorem 1 in the following form: Then f has m-symbolic dynamics in N 0 ∪ N 1 . In addition, an orbit x k corresponding, as in Definition 4, to a given sequence σ 0 , . . ., σ n , can be chosen so that it is well defined for i = 0, . . ., (n + 1)m and, moreover, x (n+1)m = x 0 .

Application of symbolic dynamics to the studied cases
In this section we apply definitions and theorems provided in Sect. 4. We want to show that the models described in previous sections admit chaotic motions in terms of m-symbolic dynamics. All the following results are numerically proved by an accurate control of propagation errors. Let us present the results in the two following sections.
Moreover, we provide explicit formulas of sets realizing the 3-symbolic dynamics. Let us consider the map P H , * in (23). The stable and unstable eigenvectors related to DP H , * at q * = (G * , r * ) = (0.345986, 0.26987) have directions, respectively, v s = (0.831003, −0.556268), v u = (0.0495113, −0.998774), and the angle between them is α = 0.296467 π. We denote as N 0 the parallelogram through q * with edges parallel to v s and v u , namely: where A 0 , B 0 are the real intervals We define two analogous parallelograms: where Then, we state the following Splitting such relations as and in view of Corollary 1, the Numerical Evidence 1 follows, with N 0 , N 1 , N 2 as in (26), (27). In Fig. 16 the construction of the covering relations is represented and in Fig. 17, it is possible to see that the sets involved are obtained by analyzing the homoclinic intersections of the stable and unstable manifolds through q * and exploiting the property of contraction and expansion, respectively, of the stable and unstable manifolds.

Existence of a horseshoe of the map P H,5
The same result in its stronger version is obtained for the map P H , in Eq. (24). In this section we use the same notations for points, intervals and sets as in the previous section but they are related, in this case, to the map P H , . We state the following Also in this case we provide explicitly the sets which form the associated topological horseshoe. Based on a couple of hyperbolic fixed points of P H , whose coordinates read q 1 = (g 1 , G 1 ) = (0.203945459, 2.06302430), and whose stable and unstable manifolds and their heteroclinic intersections are shown in Fig. 13, 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:

Numerical Evidence 4
The following covering relations hold proving the existence of a topological horseshoe for P H , .
Namely, we prove, numerically, the existence of symbolic dynamics for P H , . The obtained horseshoe associated to q 1 and q 2 with the aforementioned parameters is illustrated in Fig. 18. Looking at Fig. 13 it is possible to see how the two sets N 1 , N 2 are, respectively, contracted and stretched under the map P H , according to the stable and unstable directions.

Conclusions and perspectives
The problem covered in this article arises from the important paper by V.I. Arnold ( [1]), and it has been investigated by many authors in the last decades (see for example [2,5,7,10,13,14]). Related to the current paper are [4,16,17] where the same problem is faced by an analytical point of view in order to prove stability estimates by using normal form theory.
The idea of the numerical investigation and the application of covering relations and symbolic dynamics to this problem comes from the recent paper [11]. In the current paper, we have summarized the results obtained in previous articles (see [8,9]) underlining a method to prove numerically the existence of chaos in the full three-body problem. The authors of these papers are carrying on the study by analyzing how the results change when the nonaveraged problem given by Eq. (3) is considered. We would like also to extend the results to the spatial case and to prove the obtained results about onset of chaotic motions from an analytically point of view.
for the invitation and Ugo Locatelli for interesting discussions. The author also thanks the reviewer for the careful and detailed corrections and for the given suggestions. Finally, the author thanks the INDAM and in particular the GNFM for the support given over the years.
Funding Open access funding provided by Università degli Studi di Palermo within the CRUI-CARE Agreement.
Data Availability All data of this paper are produced by the author with a Fortran Code written by the author and all the figures are produced by the author with Gnuplot. All data could be available on request.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
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/.