Bifurcation of frozen orbits in a gravity field with zonal harmonics

We propose a methodology to study the bifurcation sequences of frozen orbits when the 2nd-order fundamental model of the satellite problem is augmented with the contribution of octupolar terms and relativistic corrections. The method is based on the analysis of twice-reduced closed normal forms expressed in terms of suitable combinations of the invariants of the Kepler problem, able to provide a clear geometric view of the problem.


Introduction
Among the manifold versions of the perturbed Kepler problem, the investigation of the gravity field expanded in multipole terms has traditionally received great attention for its relevance in applications. Therefore, several analytical tools have been developed to highlight the most important phenomena. Perturbation theory with the construction of normal forms is the standard method since the first pioneering studies (Brouwer, 1959;Kozai, 1962). The case in which only zonal terms are included in one of the settings in which we can obtain explicit approximations of the regular dynamics since the normal form is integrable. However, the presence of several parameters, both dynamical (or "distinguished" in the language of the theory of integrable systems) and physical like the multipole coefficients, hinders a global description of the dynamics. More efficient geometric and group-theoretic tools have been exploited to study the bifurcation of invariant objects when these parameters are varied (Coffey et al., 1994(Coffey et al., , 1986Cushman, 1983;Palacián, 2007). Here we study the bifurcation sequences of frozen orbits when the 2nd-order fundamental model of the satellite problem is augmented with further features of a typical planetary gravity field. We consider the contribution of the octupolar term (Coffey et al., 1994;Vinti, 1963) and the relativistic correction due to the quadrupolar term (Heimberger et al., 1990). We implement a twice-reduced normal form (Cushman, 1988;Pucacco, 2019;Pucacco and Marchesiello, 2014) which allows us to obtain in an efficient way the conditions for relative equilibria corresponding to the family of periodic orbits with fixed eccentricity and inclination. The method is tested in the 2nd-order J 2 -problem in which known results are reproduced (Palacián, 2007) and then applied to the above-mentioned perturbations. For the J 4 -problem, interesting features around the parameter values of the "Vinti problem" are highlighted with an additional family of stable frozen orbits. For the relativistic J 2 -correction, the treatment extends and completes several results obtained by Jupp and Brumberg (1991). The plan of the paper is as follows: in Section 2 we recall the model problem based on the normal from obtained after averaging with respect to the mean anomaly; in Section 3 we review the reduction methods adapted to the symmetries of the present model, discuss the version adopted here to cope with the structure of the Brouwer class of hamiltonians and show how it works in locating relative equilibria; in Section 4 we illustrate the results in concrete cases; in Section 5 we conclude with some hint for possible developments and future works. 1 arXiv:2210.12182v1 [math-ph]

Oct 2022
We are discussing some aspects of the general problem described by a Hamiltonian of the form H(L, H, G, , g, h) = ∞ j=0 j H j (L, H, G, , g, h), where H 0 is the Kepler Hamiltonian and the canonical Delaunay variables have the following expression in terms of the standard Keplerian elements (a, e, i, , ω, Ω) In the above equation, is a formal parameter, called book-keeping parameter, suitably chosen to order the hierarchy of perturbing terms (see Efthymiopoulos, 2012). Therefore, we have a perturbed Kepler problem.
Specifically, in the even zonal artificial satellite problem, we assume to start with the "original Hamiltonian" in standard Cartesian form, where q = {x, y, z}, p = {ẋ,ẏ,ż}, p = |p|, V CGF is the classical gravity field and c is the speed of light. We include the classical gravity field V CGF expanded in terms of the zonal harmonics of even degree 1 where µ = GM P is the product of Newton constant and the mass of the "planet", R P is its radius and the P k are the Legendre polynomials with sin θ = z r , r = x 2 + y 2 + z 2 .
We also add the first-order relativistic corrections following e.g. Weinberg (1972).
To simplify the structure of the Hamiltonian, we then perform a closed-form normalisation like in (Coffey et al., 1994) and (Heimberger et al., 1990). This method, inspired by works of Deprit (1981Deprit ( , 1982, has the advantage of avoiding expansions in the eccentricity and inclination (Cavallari and Efthymiopoulos, 2022;Palacián, 2002). The model in (4) is rich enough to convey several interesting dynamical features keeping the closed form structure at the lowest level of complexity. In fact, after the Delaunay reduction and the elimination of the ascending node, we deal with a secular Hamiltonian in closed form which depends on only one degree of freedom, corresponding to the pair G and g (the argument of the perigee): K(L, H, G, g) = j j K j (L, H, G, g), with L and H formal integrals of the motion. The zero-order term is clearly The first-order term is 1 In this work, we focus on the even zonal problem. Thus, only the even zonal harmonics are considered in the expansion of the gravitational potential. The complete expansion, including also tesseral terms, can be found in (Kaula, 1966).

2
The second-order term K 2 consists of two contributions: The first is related to the propagation at second order of the J 2 term in the normalising transformation (Deprit, 1969;Efthymiopoulos, 2012), The second is associated directly with the average of the H 2 term: In this work, we do not consider terms of order higher than j = 2. Hamiltonians of this type are generally denoted as "Brouwer's" ones (Brouwer, 1959;Cushman, 1983). They are characterised by the independence on the mean anomaly and the longitude of the node h (with corresponding conservation of the actions L and H), whereas the argument of perigee appears only with the harmonic cos 2g. These symmetries will all be exploited in the geometric approach described in the following. The two relativistic terms proportional to J 2 /c 2 appearing in (9) and (10) have the same structure. However, in the literature (Heimberger et al., 1990;Schanner and Soffel, 2018), they are usually kept separate and are respectively referred to as the indirect and direct term related to the non-trivial relativistic contribution of the quadrupole of the gravity field of the central body. The ordering of the perturbing terms is performed by assuming (with a certain degree of arbitrariness) the J 2 and c −2 terms to be of order and the J 4 term of order 2 , like the J 2 2 and J 2 × c −2 terms. We remark that, with a slight abuse of notation, we have denoted with the same symbols the Delaunay variables appearing in (1) and (6). We have to recall that actually they are respectively the original and the new variables related by the normalising transformation. In the present work, we are not interested in the explicit construction of particular solutions. Therefore, we will not detail the backtransformation from the new to the original coordinates. Moreover, we are not going to investigate any issue connected with the convergence of the expansions. We rely on the asymptotic properties of these series and their ability to provide reliable approximations, especially in the cases of Earth-like gravity fields. For sake of completeness, the different parts of the normalised Hamiltonian K = K 0 +K 1 +(T 2 + H 2 ) 2 , expressed in terms of the orbital elements (a, e, i, ω), are given by T 2 = 3µJ 2 2 R 4 P 128 a 5 η 7 − 5 η 2 + 36 η + 35 sin 4 i + 8 −η 2 + 6 η + 10 sin 2 i

Geometric reduction
The secular Hamiltonian in closed form in (6), while computed with an ingenious combination of tools based on the Lie transform method (Deprit, 1969;Efthymiopoulos, 2012) and the elimination of the parallax (Deprit, 1981), is nonetheless standard in being essentially an average with respect to the mean anomaly (Deprit, 1982;Palacián, 2002). However, it is liable to be treated with a group theoretically approach. It can be interpreted as a suitable combination of the invariants generating the SO(3) symmetry of the Kepler problem. In fact, the dynamics ensues from the reduction of the Hamiltonian defined on the space of the trajectories having, for the unperturbed Kepler problem with negative energy, the structure of the direct product of two spheres. The additional symmetries of the closed form of the perturbed problem are exploited to identify a regular reduced phase space with the topology of the 2-sphere. In practice, we will use a further transformation leading to a singular reduction on a surface with equivalent topology, which produces a clearer geometric view of the bifurcation sequence of frozen orbits. Here, we provide a quick reminder of the invariant theory of the Kepler problem and then apply the reduction process to perturbed Kepler problems described by Brouwer's Hamiltonians.

Invariants of the Kepler problem
Let us call G the angular momentum and A the Laplace-Runge-Lenz vector, given by with i = arccos (H/G) the orbital inclination. By defining we get the Poisson structure of the generators of SO(3) and phase-space defined by the direct product of the two 2-spheres It can therefore be imagined as the invariant space of the states characterised by given eccentricity, inclination, and arguments of perigee and node, but nonetheless equivalent for what pertains to the mean anomaly. In the unperturbed problem, the state is a given still point of the invariant space. The state point is kept moving on it by the action of the perturbation.

Reduction of the axial symmetry
Perturbed Kepler problems described by Hamiltonians of the form (6) are characterised by axial symmetry with H as formal third integral. In Cushman (1983) and Coffey et al. (1986) it is shown that, if 0 < |H| < L, the two-dimensional phase space of such problems is still diffeomorphic to a sphere. Two 4 different sets of variables, both functions of the Keplerian invariants x k , y k (k = 1, 2, 3) and suitable to analyse the dynamics, are proposed. The variables (π 1 , π 2 , π 3 ) are defined as where k = (0, 0, 1) T (see Cushman, 1983). The phase-space is then Instead, in Coffey et al. (1986), the variables (ξ 1 , ξ 2 , ξ 3 ) are introduced, defined as or, in terms of Delaunay variables, In this case, the phase-space is the sphere of radius (L 2 − H 2 )/2: The relation between the π k and the ξ k is The advantage of both these sets of variables with respect to the Delaunay variables is well explained in Coffey et al. (1986) with an imaginative metaphor. In simpler words, we can say that the Kepler reduction allows us to translate the closed form dynamics in terms of the invariants of the unperturbed problem (formal conservation of L) and the further reduction generated by the invariants ξ k is readily apt to account for the axial symmetry associated with the formal conservation of H. Recalling the description of the states of the space defined in (12), we now have that the states of (15), given a value of H, are characterised by the eccentricity and the perigee but are nonetheless equivalent for what concerns h. The dynamical evolution of the system is then determined by the intersections of the reduced phase-space S with the Hamiltonian expressed in terms of the invariants, e.g. K(ξ 1 , ξ 2 , ξ 3 ). Whenever one uses the (G, g) chart to analyse the dynamics of the closed form for given values of L and H, one excludes circular and equatorial orbits. Indeed, when either the orbital eccentricity or the orbital inclination is zero, the argument of the perigee g is not defined, thus the Delaunay variables result unsuitable to evaluate the stability of such orbits, if they are periodic as typically happens in the artificial satellite problem. Following Cushman (1983), in Iñarrea et al. (2004) it is shown that when K possesses independent symmetries of the type R 1 :(π 1 , π 2 , π 3 ) → (−π 1 , π 2 , π 3 ), R 2 :(π 1 , π 2 , π 3 ) → (π 1 , −π 2 , π 3 ), the phase-space can be further reduced, and the variables σ 1 , σ 2 , defined as are introduced, where σ 2 = G. We propose here to exploit a further set of variables, which is particularly suitable when the normalised Hamiltonian possesses symmetries of the type We introduce the variables (X, Y, Z) defined as which turn the spherical phase space S into a lemon space: This kind of reduction was proposed for the first time by Hanßmann and Sommer (2001). It is an example of singular reduction (Cushman and Bates, 1997) as opposed to the regular setting generated by the invariants ξ k . This occurs here due to the appearance of cusps in the reduced phase-space L contrary to the smoothness of the 2-sphere S. However, as it will appear clear in the following, this fact does not pose any practical issue in the induction process implemented hereafter. Even though the phase-space is still three-dimensional, we see that, in the case in which symmetries (16) are fulfilled (such as in the problem of the geo-potential when only even zonal harmonics are retained), the transformed closed form does not depend on the variable Y : K = K(X, Z). In particular, in the case of the Brouwer's Hamiltonian (6), K depends linearly on X, i.e. it is of the form where a is the set of parameters characterising the problem, including the "distinguished parameter" E. For such a problem, the analysis of the intersection of the reduced phase-space L with the function (18) is simplified by the extra symmetry of the Brouwer's Hamiltonian since, rather than working in the full 3D-space, all significant information can be obtained by projection on the (Z, X) plane. As a matter of fact, when expressed in Delaunay variables, (X, Y, Z) are equal to and, considering the structure of the normalised Hamiltonian presented in the previous section, the possibility of using the general form (18) appears immediately justified.

Equilibrium points
Relative equilibria of the reduced systems correspond to periodic orbits of the original closed form in (6), which in turn are approximations of the periodic orbits of the model problem in (1). Our main concern refers to frozen orbits which play a major role in shaping the phase-space structure of the system. They can be identified by locating "contacts" between the surfaces defined by the Hamiltonian function (18) and the lemon space L (Pucacco and Marchesiello, 2014) or in some peculiar case we will encounter in what follows if the Hamiltonian possesses a 1-dimensional level set whose intersection with 6 the phase-space produces additional (unstable) critical points. In the present subsection we describe the general procedure to locate equilibria, postponing to the next section the details of each case. Considering G as a function of Z, G = Z + (L 2 + H 2 )/2, the Poisson structure of the (X, Y, Z) variables is Henceforth, given a Hamiltonian of the form K in (18), the equations of motion are Since we are typically interested in elliptic trajectories, which implies G = 0, there exist equilibrium points whenever or The variables X, Y, Z are particularly useful in the first case when conditions (20) are fulfilled. On the (X, Z) plane, the contour of the lemon space L is C = C + C − , with For any values of the parameters a, condition (20) is fulfilled if X = 0. Thus, the normalised Hamiltonian K always possesses the equilibrium points From (19), Z = −E implies G = H; thus, the equilibrium point E 1 represents the family of equatorial orbits. Instead, Z = E implies G = L: the equilibrium point E 2 represents the family of circular orbits. Condition (20) is also fulfilled whenever a level curveX(Z; a, k) is tangent to the contour C, with k a given level set of the Hamiltonian K(X, Z; a). We can therefore have an equilibrium point of whereX is defined by recalling (18). From (23) and (24), we obtain that Z + is a zero of the function s + (Z; a) equal to Function s + (Z; p) can have multiple zeros corresponding to acceptable equilibrium solutions. In the following, we will refer to them as equilibrium points of type E + . On the other hand, we can have an equilibrium point of coordinates In this case, Z − results to be a zero of the function s − (Z; a) given by Similarly as before, equation s − (Z; a) = 0 can have multiple acceptable solutions. In this case, we are going to talk about equilibrium points of type E − . From the first of (19), we have that equilibrium points of type E + correspond to the families of periodic orbits with g = 0, π, while those of type E − correspond to the families of periodic orbits with g = ± π 2 . In the second case of (21), if there existX ∈ R andZ ∈ R fulfilling these conditions, one must verify whether the two resulting equilibrium pointsĒ 1 = (X,Ȳ 1 ,Z) andĒ 2 = (X,Ȳ 2 ,Z), with belong to L, i.e. whetherȲ 1 ,Ȳ 2 ∈ R. It is interesting to notice that for every Y the level curves of the Hamiltonian, {K = k}, given by (24), have a singularity at Z =Z as ∂K ∂X (Z) = f (Z; a) = 0. The value Z =Z gives a vertical asymptote that is a vertical plane in the 3D space X, Y, Z. The condition gives an oblique asymptote, a tilted surface in the 3D space X, Y, Z. The two surfaces cross in a straight line, orthogonal to the (Z, X) plane, which "pierces" the lemon in the symmetric fixed pointsĒ 1 ,Ē 2 .
Remark 1. Each equilibrium point E 1,2 in X, Y, Z, corresponds to one equilibrium in the variables (ξ 1 , ξ 2 , ξ 3 ), respectively equal to Instead, each equilibrium point of type E ± and of typeĒ 1,2 correspond to two equilibria. We have the following list of correspondences:

Stability of the equilibria
To study the stability of the equilibrium points, it is more convenient to come back to the variables ξ 1 , ξ 2 , ξ 3 (Coffey et al., 1994). The transformed closed form is Let us set ξ = (ξ 1 , ξ 2 , ξ 3 ) T . We havė We recall that G = G(ξ 3 ) = ξ 3 + L 2 +H 2 2 . Let us call ξ E an equilibrium point and δξ = ξ − ξ E a small displacement from it. The linearised system around the equilibrium is Since the ξ ∈ S, the solution of the previous differential system must identically satisfy the constraint Thus, we obtain the reduced system To evaluate the stability of the equilibrium point we have to compute the eigenvalues α 1,2 of DF R (ξ E ), by solving the characteristic equation with Tr DF R (ξ E ) and detDF R (ξ E ) the trace and the determinant of DF R (ξ E ). By using transformation (17), we obtain that the characteristic equation for the equilibrium point E 1 is while the one for E 2 is with s + (Z; a) and s − (Z; a) given in (25) and (27). Note that whenever the parameters a are such that an equilibrium point of either type E + or E − coincides with E 1 (i.e. Z = −E is a zero of either s + (Z; a) or s − (Z; a)) E 1 becomes degenerate. The same holds true for E 2 . For an equilibrium point of type E + of coordinates (X(Z + ; E), 0, Z + ), it can be proved that the characteristic equation is with k + the value of the Hamiltonian such thatX(Z + ; a, k + ) =X(Z + ; E). Similarly, for an equilibrium with k − such thatX(Z − ; a, k − ) = −X(Z − ; E). Since for Z = ±E,X(Z; E) > 0, the stability of the equilibrium points of type E + and E − can be determined by comparing the concavities of the level curveX(Z; a, k) and of the contour C of L at their point of tangency. Finally, the characteristic equations forĒ 1 andĒ 2 are AsȲ 2 1 =Ȳ 2 2 , the two characteristic equations coincide:Ē 1 andĒ 2 have the same stability. WhenĒ 1 andĒ 2 coincide sinceȲ 1 =Ȳ 2 = 0, the resulting equilibrium point is degenerate.
Remark 2. To evaluate the stability we can also exploit the Poincaré-Hopf index theorem: Let M be a compact manifold and w a smooth vector field on M with isolated zeros. The sum ι of the indices of the zeros of w is equal to the Euler characteristic of M (Milnor, 1965). As the phase space is a sphere in the coordinates (ξ 1 , ξ 2 , ξ 3 ), its Euler characteristic is equal to 2. Whenever E 1 and E 2 are Lyapunov stable, in the linearised reduced system they are centres; thus, their indexes ι are both equal to +1. Instead, when one of them is Lyapunov unstable, it corresponds to a saddle with ι = −1. Each equilibrium point of type E ± corresponds to two equilibrium points in (ξ 1 , ξ 2 , ξ 3 ) (see Remark 1), both either Lyapunov stable or unstable. The bifurcation of a first stable pair, implies a stability/instability transition of one of the cusps so that the indexes are (+1 − 1 + 1 + 1). The bifurcation of a second unstable pair implies that the cusp regains stability and the indexes are (+1 + 1 + 1 + 1 − 1 − 1). Due generalisation applies in the case of the points of typeĒ 1,2 .

The J 2 -problem
We are going to apply the variables X, Y, Z to analyse a classical and well-known problem in the framework of the artificial satellite theory: the study of the secular Hamiltonian in which only the second zonal harmonic of the gravitational potential is retained, i.e. the J 2 terms. From (7), (8), (9) and (10), the resulting closed form is To simplify its analysis, we make the system dimensionless by performing the following choice of units: we take the orbital semi-major axis a as unit of length and the unit of time such that µ = 1. Let us call ρ = H/L. In the adimensional system the Delaunay actions L and H become where ρ = G cos i, with i the orbital inclination. Moreover, the action G coincides with η = √ 1 − e 2 , being e the orbital eccentricity. Since in the adimensional system the planet's radius R P < 1, let us set λ plays here the role of small parameter, of the same order as the book-keeping parameter . We drop the constant Keplerian term and we perform a transformation of the time variable t → τ , defined as Thus, the secular Hamiltonian in closed form becomes In these units, the X, Y, Z variables are and we also have E = (1 − ρ 2 )/2. The introduction of the variables X, Y, Z leads to a closed form with the same structure of K in (18), with and a = (ρ; λ). Note that if the terms proportional to λ are neglected, the problem has one equilibrium solution for Since ρ = G cos i, the orbit has then a stationary pericentre at the so-called critical inclination: In the following, we study the J 2 -problem for |ρ| ∈ (0, 1) and λ ∈ (0, 1): we discuss the existence and the stability of frozen orbits by analysing the corresponding properties of the equilibrium points of the reduced system. First of all, we show that the equilibrium point E 1 , representative of the family of equatorial orbits, is always stable. Then, we analyse the stability of the equilibrium point E 2 , representative of the family of circular orbits. In particular, we determine the values ρ + and ρ − of |ρ| at which pitchfork bifurcations occur: for |ρ| between ρ − and ρ + E 2 is unstable, otherwise it is stable; moreover, there exist a stable equilibrium point of type E + for |ρ| < ρ + and an unstable equilibrium point of type E − for |ρ| < ρ − . At last, we show that the equilibrium pointsĒ 1 andĒ 2 do not exist for any ρ and λ.
For the J 2 problem and the other problems analysed in the following, all the equilibrium points of type E + are indicated with an odd integer number larger than 1 as a subscript; similarly, the subscript of the equilibrium points of type E − is an even integer number larger than 2. We recall that ρ = G cos i and G = √ 1 − e 2 . In the procedure we follow, we select a planet and we fix the value of the semi-major axis, on which λ depends through the dimensionless R P . Suppose to select a value of ρ such that an equilibrium point of type E + exists and to compute the value of the action G of such equilibrium point; from the selected ρ and the value of G we can obtain the orbital eccentricity and inclination of the family of orbits represented by the equilibrium point itself. The same holds for all the equilibrium points. In the same way, since the eccentricity of the orbits represented by E 2 is equal to zero, from ρ + and ρ − we can compute the values of the orbital inclination at which the stability of the circular orbits changes.

Stability of E 1
From (28) it results that the stability of the equilibrium point E 1 depends on the sign of s + (−E; a) and s − (−E; a). We have It is straightforward that s + (−E; a) < 0 and s − (−E, a) < 0 ∀λ ∈ (0, 1) and ∀|ρ| ∈ (0, 1). Thus, the equilibrium point E 1 is stable and it does never coincide with equilibrium points of either type E + or E − . In analogy to E 1 , from (29) we obtain that the stability of E 2 depends on the sign of s + (E; a) and s − (E; a). We have The function s + (E; a) has two real zeros ρ = ±ρ + , where Similarly, s − (E; a) possesses two real zeros ρ = ±ρ − , with It holds that ρ + > ρ − , ∀λ ∈ (0, 1). Thus, • for ρ + < |ρ| < 1 and for 0 < |ρ| < ρ − , E 2 is stable; Moreover, at |ρ| = ρ + the degenerate E 2 coincides with an equilibrium point of type E + , while at |ρ| = ρ − it coincides with an equilibrium point of type E − . If we approximate ρ + and ρ − as series in λ, we obtain Thus, the bifurcations occur nearby the zero-order solution (35) when Z = E.

Existence and stability of the equilibrium points of type E + and E −
The equilibrium points of type E + correspond to the zeros of s + (Z; a). Performing the change of variable Z → G using (34), we obtain The zeros of s + (G 2 − 1+ρ 2 2 ; a) are the zeros of S + (G; a). This is a polynomial function of degree 8 in G. Thus, finding its zeros is not straightforward. However, some pieces of information can be inferred by inverting the roles of G and ρ: we consider G as a parameter and ρ becomes the independent variable of the problem. We obtain S + (G; a) = 0 for ρ 2 = ρ 2 The solutions are admissible if 0 < ρ 2 E + 1,2 < G 2 . Since G ∈ (0, 1] and λ ∈ (0, 1), it is easy to verify that C + < 0 and is admissible ∀λ ∈ (0, 1) and ∀G ∈ (0, 1]. For G = 1 we obtain ρ 2 see the Appendix. It follows that for each |ρ| ≤ ρ + there exists only one value of G solving S + (G; a) = 0. Thus, for each |ρ| ≤ ρ + there exists only one equilibrium point of type E + , which we call E 3 and which coincides with E 2 for |ρ| = ρ + . The value of the only meaningful solution of S + (G; a) = 0 can be approximated with a perturbation method. We observe that the solution of the "unperturbed" problem with λ = 0 is just the critical value (36). Then, we can look for a solution of the form At third order in λ, we find We apply the same technique to verify the existence of equilibrium points of type E − . They correspond to the zeros of a function S − (G; a) equal to ∀G ∈ (0, 1] and ∀λ ∈ (0, 1) we have C − > 0 and 13 Thus, ρ 2 E − 1 < 0 and is not an admissible solution. Instead, ρ 2 is admissible ∀G ∈ (0, 1] and ∀λ ∈ (0, 1). For G = 1, ρ 2 see the Appendix. As a consequence, also in this case we obtain that there exists one equilibrium point of type E − for any |ρ| ≤ ρ − . We call it E 4 . For |ρ| = ρ − , it coincides with E 2 . The solution of S − (G; a) = 0 can be approximated in analogy with what seen above. At third order in λ, we find For ρ − < |ρ| < ρ + , the equilibrium E 3 is stable as a consequence of the Poincaré-Hopf theorem. By applying this last theorem, we also obtain that for 0 < |ρ| < ρ − , one equilibrium points between E 3 and E 4 is stable, while the other is unstable. Since E 3 does not undergo any bifurcation at |ρ| = ρ − , it is stable, while E 4 is unstable.

About the existence ofĒ 1 andĒ 2
The coordinatesX andZ of the equilibrium pointsĒ 1 andĒ 2 arē In order to haveZ ∈ [−E, E] it is necessary that |ρ| < 1/ √ 15. Let us call Y the square of Y coordinates of the equilibrium points,Ȳ 1,2 . It holds We have to verify whether there exist values of |ρ| It follows that Y < 0, ∀ρ ∈ (0, 1/ √ 15) and ∀λ ∈ (0, 1). Thus, the equilibrium pointsĒ 1 andĒ 2 never exist for the J 2 -problem. On the left, they are shown in the (Z, X) plane. The black line represents the contour C of the lemon space. The blue line represents the level curve tangent to C at the stable equilibrium point E 3 , while the red one represents the level curve tangent to C at the unstable E 4 . The dashed black line corresponds to Z =Z, for which the level curves have a singularity. On the right, the level curves are shown in an enlargement of the (g, G) plane surrounding the equilibrium points. The blue dots are the stable equilibrium points corresponding to E 3 ; the red curve is the separatrix of the equilibrium points at g = ± π 2 corresponding to E 4 .

Summary and comparison with previous works
Here we summarise the results for the J 2 -problem and we compare them with those previously obtained by Coffey et al. (1986) and Palacián (2007). At |ρ| = ρ + and |ρ| = ρ − , with ρ + and ρ − defined in (37) and (38), there are two pitchfork bifurcations.
In particular, we have that • for ρ + < |ρ| < 1 there exist only the equilibrium points E 1 and E 2 and they are stable; • at |ρ| = ρ + there is a bifurcation: E 2 is degenerate and coincides with E 3 , while E 1 is still stable; E 3 is an equilibrium point of type E + ; • for ρ − ≤ |ρ| < ρ + there exist the equilibrium points E 1 , E 3 , which are stable, and E 2 which is unstable; • at |ρ| = ρ − there is a bifurcation: E 2 is degenerate and coincides with E 4 of type E − , while E 1 and E 3 are still stable; • for |ρ| < ρ − there exist the equilibrium points E 1 , E 2 and E 3 , which are stable, and E 4 , which is unstable.
In Fig.1 (left panel), we show the level curves of the closed form in the (Z, X) plane when |ρ| < ρ − . The blue and red lines are tangent to the contour C of the lemon space respectively at the equilibrium points E 3 and E 4 . Note that the blue line has a concavity larger than that of C at their tangency point as E 3 is stable (see 30). Also the concavity of the red line is larger than that of C at their tangency point, which in this case implies that E 4 is unstable as follows from (31). In the right panel of Fig.1 we show the level curves in an enlargement of the (g, G) plane containing the equilibrium points. Here, E 3 corresponds to the stable equilibrium points at g = 0, π (blue dots). Instead, E 4 corresponds to the two unstable equilibrium points at g = ±π/2: the separatrix is in red. By using (44) and (48) (39) and (40), we obtain With the same transformation, by using (44) and (48), the values G ± (L, H) for the two bifurcated families can be expressed as series in J 2 . In conclusion, by exploiting the (X, Y, Z) variable and the geometrical approach we have recovered the results found in (Coffey et al., 1986) and in (Palacián, 2007).

The J 4 -problem
We study now the zonal problem containing both the J 2 and the J 4 terms. From (7), (8), (9) and (10), the closed form is Let us set After introducing it in the Hamiltonian, we adopt the same adimensional system and perform the same transformations described in Sect.4.1. Also for this problem, we obtain a secular Hamiltonian in closed form with the structure of K in (18), with g(Z, a) = −5ρ 2 + 2Z + 1 √ 2 (ρ 2 + 2Z + 1) , and a = (ρ; λ, j 4 ).
In the following, we discuss the dynamical behaviour of the problem for |ρ| ∈ (0, 1) and j 4 ∈ [−6, 6]. This range of j 4 is coherent with the book-keeping scheme used for the computation of the normalised Hamiltonian in Sect.2 and its extent allows us to include Earth and Mars. Our results are both the outcomes of analytical considerations and numerical studies. In this case, to simplify the analysis, we fix the value of λ, taking λ = 0.001. We expect the main features of the dynamics to be qualitatively similar also for other values of λ sufficiently small. First, we analyse the stability of the equilibrium points E 1 and E 2 . Then, we discuss the existence of the equilibrium points of type E + and E − and the existence ofĒ 1 andĒ 2 . Finally, we discuss their stability and we trace a bifurcation diagram. We find out that the stability of E 1 depends on both j 4 and ρ. In particular, there are ranges of j 4 in which pitchfork bifurcations occur: they affect the stability of E 1 and can give rise to either a stable equilibrium point of type E − or an unstable equilibrium point of type E + . For each j 4 ∈ [−6, 6], we also have pitchfork bifurcations affecting the stability of E 2 . We determine the values ρ + and ρ − of |ρ| at which they occur. As in the J 2 problem, E 2 is unstable if the value of |ρ| lies between ρ − and ρ + , otherwise it is stable. Following the pitchfork bifurcation occurring at |ρ| = ρ + , an equilibrium point of type E + is generated; similarly, for |ρ| < ρ − there exists an equilibrium point of type E − . The stability of these points depends on both j 4 and ρ. An interesting result is that there are ranges of j 4 where their stability changes as a consequence of further pitchfork bifurcations, which affect the existence of the equilibrium pointsĒ 1 andĒ 2 . We show that, when existing, these last ones are always unstable. Finally, we find out that for some j 4 saddle-node bifurcations also occur. They can give rise to either a pair of equilibrium points of type E + or a pair of equilibrium points of type E − . Independently of the type, one of the point of the pair is stable, while the other is unstable.
We recall once again that, for the selected planet and the fixed value of the semi-major axis (i.e. for the given j 4 and λ), the values of ρ and G of one considered equilibrium point allow us to determine the eccentricity and the inclination of the family of orbits represented by the equilibrium point itself.
In the following, we perform a general analysis not taking into account some physical limitations, for example, the fact that the orbits corresponding to a given equilibrium point may be collisional.

Stability of E 1
The stability of E 1 depends on the sign of the product s + (−E, a)s − (−E, a). For the J 4 -problem, we have . • for ρ < |ρ| < 1, E 1 is stable; • for |ρ| < ρ , E 1 is unstable.
More manageable expressions are given by the series expansions At first order they coincide with those found by Coffey et al. (1994). For j 4 = j 4 bif1 , with it holds ρ + = ρ − ; if j 4 > j 4 bif1 , ρ − > ρ + , while for j 4 < j 4 bif1 , ρ − < ρ + . As in the J 2 -problem, when the value of |ρ| is between ρ − and ρ + E 2 is unstable; at either |ρ| = ρ + or |ρ| = ρ − , it is degenerate and coincides respectively with an equilibrium point of type E + and E − . For all the other values of |ρ| E 2 is stable.

Existence of the equilibrium points of type E + and E −
We use here the same strategy applied for the J 2 -problem. After the change of variables Z → G, we obtain We have thatŜ + = 0 for ρ 2 =ρ 2 E + 1,2 , wherê For λ sufficiently small, it turns out that 0 <ρ 2 E + 2 < G 2 , ∀G ∈ (0, 1] and ∀j 4 ∈ [−6, 6]. For G = 1 it holdsρ 2 E + 2 = ρ + . Moreover, we numerically verified thatρ 2 E + 2 is increasing with respect to G. Consequently, for each |ρ| ∈ (0, ρ + ) there exists one equilibrium point, E 3 , which coincides with E 2 for |ρ| = ρ + . By applying the same perturbation method used in Section 4.1.3, we can determine the value of G corresponding to E 3 . At third order in λ we obtain The other solutionρ 2 E + 1 is only admissible for some values of j 4 . The analysis of theρ 2 E + 1 is complex and we are forced to fix the value of λ at 0.001. Anyway, we expect similar outcomes for all values of λ sufficiently small. For j 4 > j 4 bif2 , where j 4 bif2 ∼ 0.5695, there exists a range of values of G such that 0 <ρ 2 is not monotone with respect to G. Let us set For |ρ| = ρ there exists one equilibrium point E 5 of type E + . Instead, for |ρ| < ρ there are multiple equilibrium points of type E + ; they are typically two and we call them E 7 and E 9 . Also for j 4 < j 4 bif9 , it holdsρ 2 E + 1 > 0; through a numerical study, we observed that the function is increasing with G and thatρ 2 E + 1 ≤ G 2 up to a certain value of G lower than 1, for which it holdsρ 2 E + 1 = ρ 2 . Thus, for j 4 < j 4 bif9 and |ρ| ∈ (0, ρ ) there exists an equilibrium point E 11 , which coincides with E 1 for |ρ| = ρ . Concerning the equilibrium points of type E − , we have For λ sufficiently small, solutionρ 2 E − 2 is admissible ∀G ∈ (0, 1] and ∀j 4 ∈ [−6, 6]. For G = 1 we havê Moreover, we numerically verified that ∂ρ 2 E − 2 /∂G > 0. Thus, for each |ρ| ∈ (0, ρ − ) there exists the equilibrium point E 4 which coincides with E 2 for |ρ| = ρ − and whose value is Concerning the other solutionρ 2 E − 1 , its admissibility depends on j 4 . Here too, we set λ = 0.001. Through an analysis similar to the one done forρ 2 E + 1 , we reach the following conclusions: • for j 4 > j 4 bif5 , with j 4 bif5 ∼ 0.2755, at |ρ| = ρ there exists one equilibrium solution E 6 of type E − , while for |ρ| < ρ there exist typically two equilibrium solutions of type E − which we call E 8 and E 10 ; here, • for j 4 < j 4 bif6 and for |ρ| < ρ there exists an equilibrium point E 12 , which coincides with E 1 for |ρ| = ρ .

Stability analysis and bifurcation diagram
In the following we discuss the evolution of the dynamics. We set λ = 0.001, but we expect similar outcomes for all values of λ sufficiently small. We show in Fig.2 the bifurcation diagram, where the colour lines represent the values of |ρ| for which a bifurcation occurs. Some enlargements of interesting regions of the diagram are given in Fig.3. In Table 1 we summarise the bifurcations sequence and list the existing points for different ranges of j 4 ∈ [−6, 6]. Through a stability analysis based on the Poincaré-Hopf theorem, we obtain that • |ρ| = ρ − and |ρ| = ρ + are pitchfork bifurcations which cause a variation of stability of E 2 and affect the existence and the stability of the equilibrium points E 3 and E 4 ; • |ρ| = ρ , |ρ| = ρ ,bis ,|ρ| = ρ and |ρ| = ρ ,bis are pitchfork bifurcations, influencing the stability of E 3 and E 4 and the existence ofĒ 1 andĒ 2 : whenĒ 1 andĒ 2 exist, they are unstable, while E 3 and E 4 are stable; • |ρ| = ρ and |ρ| = ρ are pitchfork bifurcations affecting the stability of E 1 and the existence and stability of E 11 and E 12 ; • |ρ| = ρ and |ρ| = ρ are saddle-node bifurcations; they have no consequence on the stability of existing equilibrium solutions, but give rise to an even number of equilibrium points of type E + and E − , half of which are stable, while the other half is unstable.
To explain how to read the bifurcation diagram, let us fix a value of j 4 in the range (j 4 bif1 , 6), which is of interest for the Earth (j 4 ∼ 1.3) and Mars (j 4 ∼ 4). It holds ρ − > ρ + > ρ > ρ . For each |ρ| E 1 is stable. Moreover, • for |ρ| > ρ − , E 2 is stable; • at |ρ| = ρ − , there is a bifurcation: E 2 is degenerate and E 4 coincides with E 2 ; • for ρ + < |ρ| < ρ − , E 2 is unstable and E 4 is stable; • at |ρ| = ρ + , there is a bifurcation: E 2 is degenerate and coincides with E 3 ; E 4 is stable; • for ρ < |ρ| < ρ + , E 2 and E 4 are stable, while E 3 is unstable; • for |ρ| = ρ , E 2 and E 4 are stable and E 3 is unstable; there also exists the equilibrium point E 6 which is degenerate; • for ρ < |ρ| < ρ , E 2 and E 4 are stable and E 3 is unstable; there exist the equilibrium points E 8 and E 10 : one of them is stable, the other is unstable; • for |ρ| = ρ , E 2 and E 4 are stable; E 3 is stable; one between E 8 and E 10 is stable, while the other is unstable; there also exists the equilibrium point E 5 which is degenerate;   • |ρ| < ρ , E 2 and E 4 are stable; E 3 is stable; one between E 8 and E 10 is stable, while the other is unstable; there exist the equilibrium points E 7 and E 9 : one of them is stable, the other is unstable.
In Fig.4, we show the level curves in a neighbourhood of the bifurcations |ρ| = ρ − and |ρ| = ρ + . It is interesting to compare the phase portrait in Fig. 4d with the one shown in Fig.1 for the J 2 -problem: the concavities of the colour curves tangent to the contour of the lemon space are opposite. Indeed, in this case, E 3 is unstable and E 4 is stable. In Fig.5, we show the level curves in a neighbourhood of the two bifurcations |ρ| = ρ and |ρ| = ρ . Another significant range of values of j 4 is (j 4 bif2 , j 4 bif1 ). Here, it holds ρ + > ρ − > ρ > ρ > ρ > ρ .
In Fig.6, we show the levels curves in a neighbourhood of the bifurcations |ρ| = ρ , and |ρ| = ρ . At the bifurcation |ρ| = ρ , in the (Z, X) plane there is a level curve intersecting the contour of the lemon space at Z =Z: the intersection point is E 4 , coinciding withĒ 1 andĒ 2 . In all the range of values of |ρ| such thatĒ 1 andĒ 2 exist, there is a level curve for which Z =Z is not a singularity. When |ρ| = ρ the intersection point is E 3 . Note that for values of j 4 lower and higher than j 4 bif1 , the stability of E 3 and E 4 is different when they appear after the occurrence of the bifurcations |ρ| = ρ + and |ρ| = ρ − . A similar result was also found by Coffey et al. (1994). Here, the authors argued that this change of stability occurs at j 4 = 1, i.e. when we deal with the so called Vinti problem. Instead, we observe that the variation of the stability occurs at j 4 bif1 given by (51), which depends on λ.
To conclude, let us remark once again that the above analysis is general and it does not care about particular physical limitations. For example, one can notice that for low |ρ|, the value of G characterising the equilibrium points is typically small. This implies a large eccentricity. There is then the risk that the resulting distance of the pericentre is smaller than the central body's radius. In such a case, the resulting equilibrium cannot physically exist. For example, if we consider the case of Mars, the equilibrium points resulting from the bifurcations |ρ| = ρ and |ρ| = ρ do not exist for λ = 0.001. We study now the zonal problem containing both the J 2 and the relativistic terms. From (7), (8), (9) and (10), the closed form is We neglect here the J 4 terms to make evident the effects of the relativistic contribution. We adopt the same non-dimensional system described in Section 4.1. Let us set with λ defined in (32). We recall that λ was considered of the same order as the book-keeping parameter . Since the normalisation of the initial Hamiltonian was performed by assuming c −2 of order as well (see Section 2), j C should have a value in the neighbourhood of 1 or lower. If this was not the case, the book-keeping scheme used to compute the closed form would not be suitable anymore. Let us also remark that in the adimensional system the value of c and, thus, that of j C depend on the units of length and time, i.e. on the semi-major axis of the orbit of interest. We introduce λ and j C in the Hamiltonian. Then, we neglect the constant terms and we perform the time transformation (33). Also for this problem the resulting normalised Hamiltonian has the same structure of (18), with g(Z, a) = −5ρ 2 + 2Z + 1 √ 2 (ρ 2 + 2Z + 1) 40Z 3 + −84ρ 2 + 20 Z 2 + −74ρ 4 − 44ρ 2 − 10 Z − 11ρ 6 + 273ρ 4 − ρ 2 − 5 + 4 2 ρ 2 + 4Z + 2 −5ρ 2 + 2Z + 1 2 + j C λ(−5ρ 2 + 2Z + 1) 2 √ 2 (ρ 2 + 2Z + 1) 7 2 (−29ρ 2 + 18 2ρ 2 + 4Z + 2 − 58Z + 43), , and a = (ρ; λ, j C ). If we neglect the terms of first order in λ, we find two potential equilibrium solutions at In the following, we make some considerations about the problem considering λ ∈ (0, 1). For this problem, we perform a qualitative analysis. We find out that for j C 1, the sequence of bifurcations is the same as in the J 2 problem. On the contrary for higher values of j C , the dynamical evolution is more complex and depends on the values of λ and j C . The existence of a pair of equilibrium points of type E + , one stable and the other unstable, is triggered by a saddle-node bifurcation. The unstable point can become stable following a pitchfork bifurcation, which affects the existence of the equilibrium pointsĒ 1 andĒ 2 . The stable one can disappear following a pitchfork bifurcation, which changes the stability of the equilibrium point E 2 . A similar sequence of bifurcations occurs also concerning the equilibrium points of type E − . If none of the bifurcations affecting the stability of E 2 occur, this point is always stable. The equilibrium point E 1 is always stable.
Thus, if ρ 2 > 7/12j C , both s + (−E; a) and s − (−E; a) are positive; instead, if ρ 2 < 7/12j C they are both negative. From (28), we can conclude that E 1 is always stable. Moreover, E 1 never coincides with an equilibrium point of either type E + or E − .

About the stability of E 2
We have It holds s + (E; p) = 0 for |ρ| =ρ + , with which is positive, thus admissible, if either λ ≥ 8 49 or λ < 8 49 and j C < 16−9λ 64−392λ . We also have and s − (E; a) = 0 for |ρ| =ρ − , with ρ 2 − > 0 if either λ ≥ 8 61 or λ < 8 61 and j C < 16−5λ 64−488λ . Let us remark that for λ < 8 61 it holds 16−5λ 64−488λ > 16−9λ 64−392λ . Thus, ifρ + > 0 is an admissible solutions, alsoρ − is admissible. For each λ and j C such that bothρ − andρ + are admissible zeros of s + (E, a) and s − (E, a), it holds Let us now consider equation (29). If λ and j C are such that neitherρ − andρ + are admissible zeros, then E 2 is always stable. Also for j C =j C , E 2 is always stable, except when |ρ| =ρ + =ρ − : in this case, it is degenerate. If λ and j C are such thatρ − is an admissible solution, whileρ 2 + ≤ 0, E 2 is stable for |ρ| >ρ − , it is degenerate at |ρ| =ρ − and is unstable for |ρ| <ρ − . Finally, if bothρ − andρ + are admissible solutions, E 2 is unstable when the value of |ρ| lies betweenρ − andρ + , it is degenerate if either |ρ| =ρ − or |ρ| =ρ + and it is stable for all the other values of |ρ|. When |ρ| =ρ + , E 2 coincides with an equilibrium point of type E + . When |ρ| =ρ − it coincides with an equilibrium point of type E − . 4.3.3 About the existence of the equilibrium points of type E + and E − To discuss the existence of equilibrium points of type E + and E − , we use here the same strategy adopted for the problems previously analysed. We have We obtainS + (G; a) = 0 for ρ 2 =ρ 2 Note that for G 2 < 1/4j C ,D + > 0; instead, for G 2 > 1/4j C ,Ã + < 0. Thus, ∀λ, ∀j C , ∀G,ρ 2 and it is not admissible as solution. Whileρ 2 E + 1 > 0 if G, λ and j C are such thatD + > 0. Since is generally not monotone with respect to G. However, for j C = 0 it is equal to the same solution found for the J 2 -problem, i.e.ρ 2 (see Section 4.1.3). As a consequence, we expect that for sufficiently small values of j C ,ρ 2 E + 1 is an increasing function of G in the range of interest, i.e. G ∈ (0, 1]. In this case, for |ρ| <ρ + , there exists only one equilibrium point of type E + . Instead, for higher values of j C , such thatρ 2 E + 1 is not monotone, the outcome is different. Let us call ρ the value of |ρ| such that We have that • for |ρ| > ρ , there is no equilibrium point of type E + ; • for |ρ| = ρ , we have one equilibrium solution, which we call E 13 ; • for |ρ| < ρ there exist multiple equilibrium solutions, typically two which we call E 15 and E 17 .
Let us suppose that the Z coordinate of E 17 is larger than that of E 15 . When λ ≥ 8 49 or when λ < 8 49 and j C < 16−9λ 64−392λ , E 2 coincides with E 17 for |ρ| =ρ + . Thus, for |ρ| <ρ + , the number of equilibrium solutions reduces to one: there will exist only E 15 . In conclusion, we can infer that reducing the value of j C , the value G = G , corresponding to the maximum point of |ρ 2 E + 1 |, increases. For a fixed λ, it exists a value of j C such that G = 1, i.e. for which ρ =ρ + . Thus, for lower values of j C , the bifurcation |ρ| = ρ disappears and the only existing equilibrium point of type E + is E 15 for |ρ| <ρ + . As far as the equilibrium points of type E − , we have One can observe that for G 2 < 1/4j C ,D − > 0 and that for G 2 ≥ 1/4j C ,Ã − < 0. Thus, ∀λ, ∀j C and ∀G,ρ 2 Thus, there exist values of G, j C and λ such that is not monotone with respect to G. We find an outcome similar to the one obtained for the equilibrium points of type E + . Let us consider sufficiently high values of j C such thatρ 2 E − 1 is not monotone and let us set We have that • for |ρ| > ρ , there is no equilibrium point of the type of E − ; • for |ρ| = ρ , we have one equilibrium solution, which we call E 14 ; • for |ρ| < ρ there exist multiple equilibrium solutions, typically two which we call E 16 and E 18 .
Suppose that E 18 has a larger Z coordinate than E 16 . When λ ≥ 8 61 or when λ < 8 61 and j C < 16−5λ 64−488λ , at |ρ| =ρ − E 18 coincides with E 2 and for |ρ| <ρ − it disappears. For a fixed λ, by considering decreasing values of j C the value of G, G = G , corresponding to the maximum point ofρ 2 E − 1 , increases. Below the value of j C for which ρ =ρ − , G does not belong to the admissible range of values for G. In these cases, there only exists the equilibrium point E 16 for |ρ| <ρ − .

About the existence ofĒ 1 andĒ 2
The coordinatesX andZ of the two equilibrium points of type E + arē To haveZ ∈ [E, E], ρ 2 < min 24j C +1 15 , 7 12j C . Let us set Y =Ȳ 2 1,2 . In general, for given j C and λ, it can exists a subset of values of ρ such that Y > 0, i.e. such thatĒ 1 andĒ 2 exist. The endpoints of this range are values of ρ for whichĒ 1 andĒ 2 coincide with either an equilibrium point of type E + or E − . Let us call ρ the value of |ρ| such thatĒ 1 andĒ 2 coincide with an equilibrium point of type E + and ρ the the value of |ρ| such that they coincide with an equilibrium point of type E − . We can conclude that necessarily ρ < ρ and ρ < ρ . For j C → 0 we obtain instead the same outcome found for the J 2 -problem: for sufficiently small values of j C , there does not exist any value of ρ for whichĒ 1 andĒ 2 exist. 4.3.5 About the stability of the equilibrium points of type E + and E − and ofĒ 1 andĒ 2 Let us consider value of j C sufficiently high, such that E 15 , E 16 , E 17 and E 18 exist. We can assume that these equilibrium points are close to the equilibrium solutions (56) of the problem at order zero in λ. With this hypothesis, we can estimate their stability. To this aim, we need to assume j C ρ 2 < 1/80. At order zero in λ we obtain the same equations for the equilibrium points E 15 and E 16 , i.e.
The same holds for E 17 and E 18 : From the equations we obtain that for 7/810 < j C ρ 2 ≤ 1/80, E 15 and E 18 are unstable, while E 16 and E 17 are stable; instead for j C ρ 2 < 7/810, E 15 and E 17 are both stable, while E 16 and E 18 are both unstable. From this zero-order analysis and by applying the Poincaré-Hopf theorem we can infer the actual dynamical evolution: • for |ρ| = ρ , there exists one equilibrium solution E 13 which is degenerate; • for ρ < |ρ| < ρ , there exist E 15 , which is unstable and E 17 , which is stable; • for |ρ| = ρ , E 15 coincides withĒ 1 andĒ 2 and it is degenerate; E 17 is stable; • for |ρ| < ρ , both E 15 and E 17 are stable.
Something similar occurs concerning the equilibrium points E 16 and E 18 : • for |ρ| = ρ , there is one equilibrium solution E 14 which is degenerate; Enlargements of the regions containing the equilibrium points are performed. On the right, the level curves are shown on corresponding enlargements in the (g, G) plane. The same colour code used in Fig.4 is employed. Enlargements of the regions containing the equilibrium points are performed. On the right, the level curves are shown on corresponding enlargements in the (g, G) plane. The same colour code used in Fig.4 is employed. : Level curves for the J 2 -problem with relativistic term with j C = 0.2 and λ = 0.001 and for a value of |ρ| <ρ + .Enlargements of the regions containing the equilibrium points are performed. On the right, the level curves are shown on corresponding enlargements in the (g, G) plane. The same colour code used in Fig.4 is employed.

Compliance with ethical standards
Conflict of interest: The authors declare that they have no conflict of interest.