Chaotic coexistence of librational and rotational dynamics in the averaged planar three-body problem

Through an appropriate change of reference frame and rescalings of the variables and the parameters introduced, the Hamiltonian of the three-body problem is written as a perturbed Kepler problem. In this system, new Delaunay variables are defined and a suitable configuration of the phase space and the mass parameters is chosen. In such a system, wide regions of librational and rotational motions where orbits are regular and stable are found. Close to the separatrix of these regions, the existence of chaotic motions presenting a double rotational and librational dynamics is proved, numerically, through Poincaré sections and the use of FLI.


Description of the model
We fix an orthogonal reference frame (i, j, k) in the Euclidean space, which we identify with R 3 . In such a space, we consider three bodies P 1 , P 2 , P 3 with masses m 1 , m 2 , m 3 interacting through the gravity only. The Hamiltonian H 3b describing the dynamics of the three bodies can be written as where x 1 , x 2 , x 3 ∈ R 3 are the position coordinates of the three bodies and y 1 , y 2 , y 3 ∈ R 3 are their respective linear momenta; · denotes the Euclidean distance and, finally, the gravity constant has been conventionally fixed to one. In our study, the motion of the three bodies is confined in a plane, so it has been called planar problem, choosing positions and impulses x 1 , x 2 , x 3 , y 1 , y 2 , y 3 ∈ R 2 . We reduce the translation symmetry relating the positions of two (out of three) masses to the position of the third one, as described in ( § 5 Giorgilli 2008).
Contrarily to the usual practice, we choose P 2 as reference body with mass m 2 = μ (usually, the unit mass is chosen) as depicted in Fig. 1. With such choice, the Hamiltonian governing the motions of the P 1 , P 3 with masses m 1 = 1 and m 3 = κ, respectively, is H 3b (x, x , y, y ) = κ + μ κμ where x, x ∈ R 2 are the position coordinates of P 1 and P 3 , respectively, and y, y ∈ R 2 are their respective linear momenta; the body P 2 is at the origin of the reference system. We perform a rescaling of coordinates as follows: (x, x ) → κ + μ μ 2 κ 2 (x, x ), (y, y ) → μ 2 κ 2 κ + μ (y, y ), t → μ 3 κ 3 κ + μ t (with t denoting the time). It does not alter the equations of motion provided that H 3b is changed to with α := κ + μ κμ(μ + 1) , β := κ + μ κ 2 (μ + 1) , γ := 1 μ + 1 , δ := κ(μ + 1) κ + μ , being the current mass parameters depending on μ and κ. The first two terms y 2 2 − 1 x are called Keplerian part and, assuming it takes negative values, it generate motions for the body of coordinates (x, y) on an ellipse. The Hamiltonian (2) has four degrees of freedom, and in the following step, we aim to reduce it to a three degrees of freedom system. 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 C = x × y · k + x × y · k where k = i × j is the unit vector orthogonal to 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).
We denote by E the ellipse generated by the Keplerian part of Eq.
(2) for a given initial datum (x, y). Finally, we define a couple of Delaunay variables ( , G, , g) ∈ R 2 × T 2 for the body P 3 with respect to P 2 as follows: where a is the semi-major axis of E, -G = x × y · k is the angular momentum of the third body P 3 , -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 couple coordinate (R, r ) ∈ R 2 for the body P 1 as: namely R is the radial velocity of P 1 and r is the Euclidean length of x . We underline that the coordinates (R, G, , r , g, ) refer to a frame moving with P 1 and they have been introduced to reduce the symmetry of rotations (we refer to Pinzari 2019 for an exhaustive definition of the variables). The Hamiltonian (2) written in the new coordinates appears as where the mean anomaly is related by the eccentric anomaly ξ by means of the Kepler Equation and e is the eccentricity of the ellipse E given by R, G, ξ, r , g) is the Newtonian potential and I = I ( , R, G, ξ, r , g) is the indirect term. We rewrite Hamiltonian (3) as where −1/(2 2 ) is the Keplerian term, and are, respectively, the kinetic term, the disturbing term, the Newtonian potential term and the indirect term. We remark that after the symmetry reduction, the Hamiltonian depends parametrically on the total angular momentum C.
Hamiltonian (3) is a 3 degrees of freedom system and, as it has been done in Di Ruzza and Pinzari (2022), we split the term inside parentheses as the sum of its -average (denoted by H C ) and the zero-average part (denoted by H C ): ( , R, G, , r , g) . (6) We claim two main assumptions. The first is that the Keplerian term is much more bigger than the zero-averaged term: where · is some norm of functions. Under condition (7), perturbation theory (see for example Arnold 1963) allows to conjugate Hamiltonian (6) to where O 2 represents a remainder term depending on all the coordinates. Let us now consider the system (8) when the remainder is neglected: the Keplerian term in eq. (8) becomes an additive term which does not change the equation of motions, so we can consider as a constant parameter and, without loss of generality, we can fix it equal to 1.
Moreover, we can reabsorb the parameter δ through a change of time, and the system can be reduced to a new Hamiltonian H C , which is given by where Hamiltonian H C (R, G, r , g) is a 2 degrees of freedom system. With the notations defined in (5), Hamiltonian (9) can be rewritten as Now, we state the second assumption, namely we ask for Remark 1 Note that term P = γ y · y in Eq. (5) is a -zero-average term, so it is part of the term H C together with the -zero-average part of the Newtonian potential U pot .
The main purpose of this article is to numerically study the dynamics of the systems governed by Hamiltonian (11) when Conditions (7) and (12) are satisfied. Left: orders of magnitude of terms K 0 in dark-green, K 1 in blue, U in red, P in light-green, multiplied by δ, referred to the Hamiltonian (4); right: orders of magnitude of terms K 0 , K 1 , U , respectively in dark-green, blue, red, referred to the Hamiltonian (9)

Energy terms comparison
Fixing G 0 , ξ 0 , g 0 , in the left panel of Fig. 2, we can see the orders of the Hamiltonian terms given in Eq.(5) propagated for a given period under the flow given by Eq. (4). Changing G 0 , the four terms modify a little without changing the orders of magnitude. Changing ξ 0 and g 0 we just get horizontal translations (in time phasing) of the figure. It is easy to verify that Assumptions (7) and (12) are fulfilled recalling that = 1 implies − 1 2 2 = 1 2 . In an analogous way, we define parameters and initial conditions for the averaged Hamiltonian (9), choosing as parameters μ = 2, κ = 0.1, which implies α = 3.5, β = 70; C = 10 2 ; = 1 and ad hoc initial conditions as: for some G 0 ∈ (−1, 1) and g 0 ∈ [0, 2π]. In the right panel of Fig. 2, we can compare the orders of magnitude of the terms K 0 , K 1 , U in such case. In this section, we analyze the motions taken place in this configuration, in the averaged system.

Analysis of the three-dimensional phase space
Let us start by describing the motion in the averaged system ruled by Eq. (9). In our computation, the Newtonian averaged potential U (Eq. (10)) is expanded in power series of the small parameter a/r = 1/r and the series is truncated at order 10, which, after a check of the error of truncation, is chosen as good compromise between precision and working time of the machine. We fix the total energy, namely we fix the value of the Hamiltonian (9) taken in suitable initial conditions so to reduce the phase space from dimension 4 to 3. At the beginning, we consider initial conditions in Eq. (14) with G 0 = 0.8 and g 0 = π/2 and analyze the three-dimensional phase space (g, G, r ). Such a space appears foliated by regular invariant two-dimensional manifolds where the motions take place. To better visualize the threedimensional structure, we propose, in both Figs. 3 and 4 , two different simulations: in the upper panels, we show the propagation of the trajectories with U truncated at order 2 and in the bottom panels the propagation of the trajectories with U truncated at order 10 (it has been done just for a better visualization). Let us describe the first case: as it can be seen in the upper left panel of Fig. 3, the trajectories belong to cylinders orthogonal to the (g, G)-plane; the variable g circulates in [0, 2π] when the variable G ∈ I c and librates when G ∈ I l for suitable real intervals I c , I l which will be analyzed in the following. In the upper right panel, a zoom on the variable G allows to see the cylinder which librating trajectories belong to. At the center of the librating island, an elliptic trajectory, represented as a magenta segment, lies in the three-dimensional space: this trajectory is associated with the elliptic fixed point of coordinates (g, G) = (π, 0.9985968919390512) of the Poincaré map defined in the next section; the variables g, G of such trajectory remain constant, while the variable r oscillates along the real interval I r which will be specified in the following. In the upper left panel of Fig. 3, it is also highlighted in orange a particular trajectory named last trajectory. This trajectory is peculiar because it belongs to a unidimensional curve parallel to the (g, G)plane, namely the variable r remains constant during the motion. We called it last trajectory because, when g = π, the variable G takes the minimum value G min of the interval I c and for values smaller than G min no motion takes place for the fixed value of the energy: we mean that, for the fixed value of the energy, initial conditions (g, G, r ) are not admissible for G < G min .
In the upper left panel of Fig. 4, the projection on the (G, r )-plane of the trajectories is represented. In such figure, we can well see the last trajectory in orange and the elliptic trajectory in magenta and how the range of the variable r increases from the former to the latter. In the right panel of Fig. 4, the projection on the (r , R)-plane of the trajectories is represented. In this figure, we highlight again with an orange point the last trajectory whose variables r = r last and R = 0 remain constant and in magenta the elliptic trajectory where the variation of the variables r and R is maximum. In the bottom panels of Figs. 3 and 4 (where the same simulations with U truncated at order 10 are performed), we can see the Configuration such that the variable g = π , meaning that the pericenter P of the orbit of the body P 3 is on the line of the bodies P 1 and P 2 , beyond P 2 manifolds just little deformed and displaced with respect to the previous case. The elliptic trajectory shown in magenta in the bottom right panel of Fig. 3 does not appear as a segment but as a closed curve with a very small variation in the variables g, G as well as the last trajectory shown in the bottom left panel of Fig. 3 has a very small but nonzero variation in the variable r . The same imperceptible but present differences can be seen in Fig. 4 between the upper and the bottom panels.
Let us show a physical interpretation of such results. The elliptic trajectory corresponds to the equilibrium configurations g = π, namely the pericenter P of the orbit of the body P 3 is on the line of the bodies P 1 and P 2 , as depicted in Fig. 5. This stable equilibrium corresponds to a 1 : 1 resonance between the slow orbital revolution of the body P 1 and the precession of the pericenter P of the orbit of the body P 3 . In such stable equilibrium configuration, the variable r , which is the distance between the bodies P 1 and P 2 , oscillates along the interval I r . Librations of g around the value π correspond to librational motions about the 1 : 1 resonance just described. Circulations of g correspond to complete revolutions of the body P 1 with respect to the pericenter P of the orbit of the body P 3 ; in such case, varying the value of the angular momentum G from G min to 1, the variation of the variable r varies from 0 (when G = G min in the last trajectory) to the length of the interval I r .

The two-dimensional phase space and the Poincaré section
All the computation from now on have been done with the averaged Newtonian potential truncated at order 10.
Fixed the energy, we want to reduce the phase space from dimension 3 to dimension 2; we fix a bi-dimensional plane in the three-dimensional space (g, G, r ) parallel to the (G, g)plane such that the variable r = r last , and we construct a Poincaré map as first return map on this fixed plane. Fixed g = π, we call I c = [G last , G c ) the real interval such that the trajectories circulate and I l = (G l , 1) the real interval such that the trajectories librate, where G last < G c < G l < 1.
The Poincaré map shown in Fig. 6has a very regular structure where we can recognize regular invariant tori circulating for G ∈ I c and closed curves librating for G ∈ I l . In the left panel, we can see all the admissible range of the variable G, say I G and we see that, when g = π, for G < G last no orbit appears, as discussed in the previous section; in orange we see the section of the last trajectory and with the magenta point we can see the section of the elliptic trajectory; this point corresponds to the elliptic fixed point of the Poincaré map computed with a Newton algorithm. In the right panel, we see a zoom on the variable G Configuration such that the variable g = 0, meaning that the pericenter P of the orbit of the body P 3 is between the bodies P 1 and P 2 where the region between circulating and librating trajectories is represented. Let us remark that, for a given value of G close to 1, the librations end at g = 0 (that is also g = 2π). As depicted in Fig. 7, this corresponds to a configuration where the pericenter P of the orbit of the body P 3 is on the line between the bodies P 1 and P 2 and it is an unstable configuration.
We are particularly interested in the analysis of the region close to the separatrix which splits the phase space from the circulating to the librating zones. We wonder if there exists one or more orbits with initial conditions (π, G, r ), with G ∈ (G c , G l ) and for some r ∈ I r close to the separatrix and in the affirmative case we ask about the behavior of such orbits, namely if they are regular or if they show some kind of chaotic motion.
The Newton algorithm does not allow to find any hyperbolic fixed points as well as the analysis of the three-dimensional structure of the orbits does not provide elements to prove the existence of any suitable periodic orbit which could be associated with a hyperbolic fixed point. In the following section, we analyze the region close to the separatrix by means of an important tool of chaos indicator, namely the Fast Lyapunov Indicator.

Fast Lyapunov indicator analysis
Let us briefly recall some definitions and apply them to our problem. Let us consider the differential equationẋ = F(x) and let v(t) be the solution of the variational equation We define as Fast Lyapunov Indicator (FLI) of an initial condition (x 0 , v 0 ) at time T (see, for example Lega et al. 2016), the quantity and we denote by N the logarithmic variation of the tangent vector as In our case, x = (R, G, r , g) and the differential equation is given by the Hamilton equations: where H C is defined in (9). Through the use of the chaos indicators, we intend to provide an accurate analysis of the transition zone between the circulating and the librating orbits, proving, numerically, the existence of weakly chaotic orbits in which circulating and librating phases coexist. Let us begin with this analysis. As in the previous section, we fix an energy value reducing the phase space from dimension 4 to 3; then, we fix a plane parallel to the (g, G)-plane, such that, r = r last . We produce a two-dimensional grid of initial conditions (g, G) and we complete each two-dimensional initial condition (g 0 , G 0 ) to a four-dimensional initial condition (R 0 , G 0 , r 0 , g 0 ), where r 0 belongs to the plane and R 0 are such that all the initial conditions are isoenergetic; we propagate them under the Hamiltonian flow given by the map (16) computing the FLI of each orbit for a given time T . Then, we provide plots showing the FLI values with a color code according to their values and projected on the section (g, G) providing a stability chart. Stable orbits correspond to dark region, while chaotic orbits appear in yellow. Figure 8, upper panel, shows a grid of 1000×1000 equally spaced initial conditions (g, G) where g ∈ [0, 2π); G ∈ [G last , 1) and T = 10 4 (the period, namely the time of first return on the Poincaré section of the elliptic three-dimensional phase space trajectory defined in Sect. 3.2, is about t = 1000, thus T is the time of about 10 periods). We can see that the admissible domain turns out to be very regular and stable; only in the highest part of the plot a thin curve appears with color coded less stable. The highest part is highlighted with a light-blue box. In the bottom panel, we enlarge this zone taking a grid of 1000 × 1000 initial conditions (g, G) where g ∈ [0, 2π) and G ∈ [0.994, 1). In this case a sort of separatrix splitting the phase space into two regions is well visible. Comparing Fig. 6 with Fig. 8, it is evident that the separatrix splits the regions of circulating and librating orbits.
Let us continue by analyzing some particular orbits. As it is well-known, in stable orbits, FLIs have a logarithmic growth in time, while chaotic motions present FLIs with linear growth in time (see, for example, Froeschle et al. 2000). In Fig. 9, left panel, we consider orbits in all the admissible domain and represent them through their Poincaré section as computed  figure   Fig. 9 On the left, we highlight with different colors some orbits on the Poincaré section. On the right, we represent the evolution of the FLI versus time of initial conditions associated with the orbits represented in the left panel with the same colors, respectively  Fig. 9, on the left, we highlight with different colors some orbits on a smaller region of the Poincaré section. On the right, we represent the evolution of the FLI versus time of initial conditions associated with the orbits represented in the left panel with the same colors, respectively before. We highlight some regular orbits, respectively, from bottom to top, in orange the last orbit, in red one rotational orbit, in blue a rotational orbit closer to the separatrix, in sky-blue a librational orbit close (but not very close) to the separatrix and, finally, with a magenta point the elliptic periodic orbit. In the same figure, in the right panel, we report the variation in time, for T = 10 4 of the FLIs of the five orbits using the same color to represent the same orbit. What we can see is that all the selected orbits have a small growth from zero to a small relative value compared with the values of the legend in Fig. 8 and then that values remain constant for all the considered time T . The elliptic orbit (magenta) turns out to be the more stable orbit and the two orbits approaching to the separatrix have the FLI values little bit higher. Next, we do the same analysis to the orbits closer to the separatrix. In Fig. 10, in the left panel, we plot the Poincaré section of some orbits: in blue a rotational orbit and in light-green a librational one, both of them with very small FLI values (about 0.13) shown with the same color in the right panel. Another librational curve plotted in dark-green closer to the separatrix has a little bit higher FLI value (about 0.4). Finally, we plot three orbits (one in dark-red, one in pink and one in dark-blue) appearing overlapped and very close to the separatrix in the left panel of Fig. 10. Let us show in the right panel the evolution in time of the FLIs of these three orbits; the dark-red orbit corresponding to a rotational motion and the pink corresponding to a librational one have a greater growth with respect to the previous curves and both of them reach the same value (about 2) in a time equal to T . The dark-blue curve in the left panel is completely hidden by the previous two curves when G is close to 1, and in the right panel we can see the growth of its FLI in time is too much more greater than all the others. This orbit and the variation of its FLI in time has been studied extensively in Sect. 4.1.

Analysis of orbits and respective FLI close to the separatrix
Let us concentrate the analysis on the three orbits plotted in dark-red, pink and dark-blue depicted in Fig. 10. We already told about the behavior of the dark-red and pink, being, respectively, rotational and librational orbits. The peculiarity of the dark-blue orbit is that it changes its dynamics from librational to rotational and vice versa. Let us show in Fig. 11, the three-dimensional projection on the space (g, G, r ) of the three isoenergetic orbits we are talking about. In the upper left panel, we plot the dark-red rotational orbit: the variable g circulates from 0 to 2π; the variable G oscillates very close to 1, in particular its value is almost 1 when g ∈ (0, π/2) ∨ (3/2 π, 2π) and different from 1 when g ∈ (π/2, 3/2 π, ); the variable r oscillates in the admissible range I r (even if it is not well visible in this plot). In the upper right panel, the pink librational orbit is depicted: the variable g oscillates around π, in particular librates from π/2 to 3/2 π; the variable G oscillates very close to 1, assuming the value almost 1 for a part of its period and less close values to 1 in the remaining part of the period; the variable r , as before, oscillates in the admissible range I r . In the bottom left panel, we plot in dark-blue an orbit where rotational and librational motions coexist. In particular the variable g changes its oscillation in time passing from a rotational motion to a librational one, and vice versa. In the bottom right panel, the previous three orbits are plotted together to see the relation to each other. To better see how the three orbits behave in the three-dimensional space (g, G, r ), we plot, in Fig. 12, a two-dimensional (upper) and a threedimensional (bottom) projections on the (g, G)-plane and (g, G, r )-space, respectively: the blue orbit lies in a two-dimensional manifold constrained between the red and the pink ones. Orbits close to the blue one have strong sensitivity about initial conditions and this is an indication of the chaotic behavior of such kind of orbits. But, being a two degrees of freedom system, the chaos is bounded by the manifolds on which the rotational and librational orbits lie.

Fig. 12
Upper and bottom, respectively, a zoom of a two-dimensional (on the (g, G)-plane) and a three-dimensional (on the (g, G, r )-space) projection of the three orbits represented in Fig. 11

Analysis of the orbits changing their dynamics
Let us continue the analysis of the orbit with double behavior changing from rotational to librational dynamics and vice versa. We propagate this orbit (the dark-blue in previous figures) for a much more longer time, say T = 3 · 10 6 . As we can see in the left panel of Fig. 13, the three-dimensional projection shows very well the double nature of this orbit. The variable g changes in time from circulating from 0 to 2π to librating around π and it is well represented in the left panel of Fig. 15, where g versus time is plotted. In the left panel of Fig. 13, we can see that, when the variable G is close to 1, the orbit is splitted into rotational and librational regimes; in particular, there exists a value of the variable r , say r m , such that when G is close to 1, the orbit is in a rotational regime for each r > r m and in the librational regime for each r < r m with r ∈ I r .
In the right panel of Fig. 13, we plot the variation of FLI of such orbit versus time. If the trend appeared to be with linear growth for a shorter time as shown in Fig. 10 (blue line in the right panel), for a longer time, the growth stops and the behavior of the FLI is completely constant. This seems to be in contrast with the chaotic behavior of this orbit, but, as told  before, the orbit, although chaotic, cannot move in space because its manifold is constrained by the other manifolds. To analyze this phenomenon, we have done a study of the logarithmic variation of the tangent vector defined in Eq. (15).
In Fig. 14, we plot the variation of N (t) of the considered orbit. In the upper left panel, where we plot the function N versus time, what is very evident is a "sort of" periodic behavior over time of this function. The function N has peaks which highlight the chaotic behavior of the orbit in given moments, and the remaining time, has low values typical of regular orbits. We are interested in understanding when these peaks come out. Thus, we plot, in Fig. 14, upper right, bottom left and bottom right the variation of the function N versus the variable g, G, r , respectively. With respect to the variable g, the higher values are taken for g = 0, π, 2π; with respect to the variable G, we have just one peak, in fact all the higher values are concentrated in the same value of the variable G very close to 1; with respect to the r variable, we also have one peak and all the higher values are concentrated in a value of r , say r = r m . Comparing these plots with the orbit plotted in Fig. 13, we can see how the higher values of the function N , are concentrated in the zone where the orbit changes its behavior from rotational to librational and vice versa. At the same time, the values of the function N are very low in the phases of rotational and librational motion. This shows the extremely regular behavior of the orbit when it is in the rotational and librational regimes and a chaotic behavior when the transitions from one regime to another occur.
In Fig. 15, we highlighted this aspect. In the left panel, the variation of the variable g versus time is plotted for a time t = 3 · 10 6 and we see the change from the circulation from 0 to 2π to the libration around π (and vice versa) of such variable. In the right panel, we overlap two different quantities: in pink, we plot the variation of variable g versus time, and in dark-green, we plot the variation on the function N versus time, both for a longer time t = 10 7 . The plot shows as the green peaks correspond to the transitions of the variable g from a regime to the other.
We want to show that this "sort of" periodicity of transitions strictly depends on the initial conditions. We say a "sort of" periodicity because the number of oscillations between one regime and the other is not constant. Let us recall the FLI map represented in the bottom panel of Fig. 8, and as it is shown in the left panel of Fig. 16, we fix a vertical segment in the (g, G)-plane depicted in sky-blue, such that g = π and G ∈ [0.99434, 0.9944] and we compute the FLI for 10 6 initial conditions on such segment. In the right panel of Fig. 16, we plot the FLI at time T = 10 4 versus the variable G, and we get a zone (although small) Fig. 16 Left: FLI map of a grid of initial conditions containing the separatrix. A small vertical segment intersecting the separatrix in g = π is represented in sky-blue. Right: Value of FLI of initial conditions belonging to the sky-blue segment on the left panel of initial conditions where the FLI grow more than the initial conditions in the rest of the domain. There are different peaks and, in that zone, the orbits present the double rotational and librational dynamics.
We compare two different orbits associated with different initial conditions; we show, for example, the two initial conditions corresponding to the first two peaks in the right panel of Fig. 16. We call the two orbits A and B with, respectively, the initial conditions: The distance in the four-dimensional space (R, G, r , g) of conditions A and B is about 1.9642 · 10 −6 . We propagate the two initial conditions to get orbits A and B and, in Fig. 17, we plot the evolution of the variable g versus time of the two orbits for a short time. In the upper panels, in blue, the orbit A is represented and in the bottom panels, in red, the orbit B is depicted. In the left and right panels we plot two different time spans, from 0 to 15,000 and from 400,000 to 415,000, respectively. In the highlighted boxes, we show that the number of circulations from 0 to 2π and the librations around π changes from orbit to orbit; namely, the period between transitions is not constant for different initial conditions.

Changing orbits analyzed according to the variation of the variable r
The changing orbits analyzed in the previous section have initial conditions in the chaotic zone such that the value of the variable r is at about the middle of the admissible interval I r , which corresponds to the value r = r last . We want to analyze how the changing orbits behave when the value of the variable r changes. We have therefore produced new FLI maps in a vertical plane, i.e. in the (G, r )-plane for a given constant value of g. At the beginning, we plot a grid of 1000 × 2000 initial conditions (G, r ) in the plane g = π for (G, r ) ∈ I G × I r .
We show the FLI map in Fig. 18. This plot is interesting for a global view of the admissible domain, where we could see a very regular structure unless a vertical line close to G 0.994 and the bottom and the upper borders. Then, we have made a study of the filtered FLI (see Guzzo and Lega 2014), i.e. we considered again a grid on the (g, G)-plane and we have computed the FLI of each initial datum only when its orbit passed through determined regions of the phase space. In this way, we identified neighborhoods of the phase space where the FLI are higher than in the other areas. We identified, for example, the region (g, G) ∈ (1.61, 1.63) × (0.99998, 1) and in Fig. 19 , we show the FLI map computed with a grid of 1500 × 1000 initial conditions (g 0 , G 0 ) ∈ (1.56, 1.67) × (0.99995, 1).
At this point, we compute the FLI on a grid of 1000 × 2000 initial conditions (G 0 , r 0 ) in the domain [0.9999, 1) × I r , with g = 1.62 and we plot in Fig. 20the FLI map. We can see a vertical strip of initial conditions close to 1 where the FLI assume the higher values and a chaotic structure appears. We want to show the existence of changing chaotic orbits in all this vertical strip and how they change according to the initial value of the variable r .

Inferior orbit
Let us choose an initial condition in the vertical chaotic strip of Fig. 20 such that the variable r has the minimum admissible value in I r . Let us call the associated orbit "inferior orbit". The initial condition is given as follows: In the left panel of Fig. 21, we plot the three-dimensional projection of the orbit; as we can see the double dynamics is still present, although it appears different from the orbit represented in Fig. 13. In this case, during the phase of librational dynamics, the variation in r when the variable G is close to 1 appears very very small, while in the rotational regime the variation in r appears to span in all the range I r . In the right panel of Fig. 21, we plot the variation in time of the FLI associated with this orbit and, as before, we see that, for a long In the upper left panel of Fig. 22, we plot the function N versus time. Also in this case, we see peaks with a sort of periodicity. The peaks correspond to the instants when the chaotic behavior of the orbit comes out, namely when it changes its dynamics from rotational to librational or vice versa. In the upper right, bottom left and bottom right, we plot the function N versus the variables g, G, r , respectively. In particular, we remark that the highest values of the function N are concentrated in only one peak in G 1 when the function is plotted versus G and in r = r inf when the function is plotted versus r . In the left panel of Fig. 23, we plot the evolution of the variable g of the inferior orbit versus time. The variable changes its regime from circulating from 0 to 2 π to librating around π. The phase of libration appears with a sort of periodicity for a very short time and it is highlighted in Fig. 24. It is completely in agreement with the three-dimensional visualization of the orbit represented in the left panel of Fig. 21. Moreover, in the right panel of Fig. 23, we overlap the variation in time of the variable g depicted in pink with the variation in time on the function N colored in dark-green. We can appreciate the superimposition of the changes of the variable g with the peaks of the function N . We want to better show this overlapping: in Fig. 24, we represent a zoom of the right panel of Fig. 23. We can see that each green peak of Fig. 23 corresponds, if enlarged, to two different peaks and they coincide with the transition from circulating to librating and vice versa of the variable g.

Superior orbit
The same analysis has been done for orbits with initial conditions in the chaotic zone such that r = max I r = r sup . We choose and show the orbit with initial condition: Let us compare the left panels of Figs. 21 and 25. In the superior orbit, we can observe that, in the phase of rotational dynamics, the variation of variable r when the variable G is close to 1 is very very small, while in the librational regime the variation of r spans in the whole  Fig. 25 shows the growth of FLI of the superior orbit. It does not look to be constant, but if we increase the propagation time it becomes constant. As for the previous orbits, we analyze the variation of the function N (t). Figure 26shows, again, peaks where the chaotic nature of this orbit comes out. In particular, from the bottom panels we see that high values of N (t) are concentrated, in this case, when G is close to 1 and r = r sup in agreement with the three-dimensional visualization of the orbit.
In the left panel of Fig. 27 , we plot the variation of the variable g versus time and the double regime of circulating and librating motions is visible. Moreover, in the right panel of Fig. 27, we overlap, again, the variation in time of the variable g depicted in pink with the variation in time on the function N colored in dark-green. Also in this case, we can appreciate the superimposition of the transitions of dynamics of the variable g with the peaks of the function N . To better represent this overlapping, in Fig. 28, we show a magnification of the right panel of Fig. 27. We can see that each green peak of Fig. 27 corresponds, if enlarged, to two different peaks and they coincide with the transition of the variable g from librating to circulating and vice versa.

Conclusions and perspectives
We start from a standard formulation of the Hamiltonian of the full (meaning non-restricted) three-body problem (see Eq. (1)) and, introducing suitable transformations and new mass parameters, we rewrite the Hamiltonian, first, in the form of Eq. (4) for the non-averaged system, then, as Eq. (9) for the averaged case. In this context, we provide the existence of two wide regions of initial condition phase space where the motions are regular and stable and they can be librational or rotational. Close to the separatrix of these two regions, we find a zone of initial conditions where rotational and librational motions coexist for the same initial condition in a chaotic way.
These results could be related to the notion of "stickiness", which are studied in several dynamical systems in different papers as, for example (Contopoulos et al. 1997;Dvorak et al. 1998). Sticky orbits can be found in the case of restricted three-body problem, for example in Tsiganis et al. (2000), where the authors studied the trojan asteroid Thersites which is in a stable L 4 librating position; they show, through many simulations, a high probability of the asteroid to "escape" from its equilibrium, to unstable orbits or to jump from L 4 to L 5 and it is called, for this reason, a "jumping" trojan. The same was proved in Sidorenko (2018) where the author shows also one more asteroid jumping from the libration about L 4 to libration about L 5 (thus, the author provides two examples, one for Jupiter and one for the Earth). Also in the restricted four-body problem (see for example Soulis et al. 2008), the authors prove the existence of a structure of ordered motion surrounded by sticky and chaotic orbits as well as orbits which rapidly escape to infinity. Two important issues differentiate the results of the current article with the notion of stickiness: first, the last cited articles deal with the restricted problem while the current paper concerns with the full problem; second, usually, stickiness time corresponds to the time to escape to infinity, while in the current paper the chaotic region is constrained by stable orbits so that the chaotic motion can not escape and can not diffuse. Nevertheless, it would be worth investigating the concept of stickiness in the full problem and looking for the existence of escaping to infinity orbits in the non-averaged case or in the spatial case.
However, the results obtained in this paper could be very interesting also in a wider context. The choice of our configuration is due in order to have the right weight for each term which composes the Hamiltonian and they depend on the parameters μ, κ, r , C. Thus, one of the future goals is to determine, in a constructive way, the range of these 4 parameters where: -the system is stable; -the chaotic behavior appears; -application to real celestial problems can be possible (for example, real asteroid systems, natural and artificial satellite motions, exoplanetary problems, hierarchical triple star/black hole systems, the Lidov-Kozai mechanism appearing in the spatial case); -comparison between averaged and non-averaged Hamiltonian is possible.
The last one is already a work in progress, in fact, at the beginning of Sect. 3, we propose an initial configuration (Eq. (13)) such that both Assumptions (7) and (12) are fulfilled. Thus, we wonder if in the non-averaged system we have the onset of chaos and, being a three degrees of freedom problem we ask if and how the chaos could diffuse in the phase space.