Numerical Investigation for Periodic Orbits in the Hill Three-Body Problem

: The current work performs a numerical study on periodic motions of the Hill three-body problem. In particular, by computing the stability of its basic planar families we determine vertical self-resonant (VSR) periodic orbits at which families of three-dimensional periodic orbits bifurcate. It is found that each VSR orbit generates two such families where the multiplicity and symmetry of their member orbits depend on certain property characteristics of the corresponding VSR orbit’s stability. We trace twenty four bifurcated families which are computed and continued up to their natural termination forming thus a manifold of three-dimensional solutions. These solutions are of special importance in the Sun-Earth-Satellite system since they may serve as reference orbits for observations or space mission design.


Introduction
It is well-known that the restricted three-body problem may serve as a basic dynamical model for the study of many real systems in the field of Celestial Mechanics. One of its benefits, in this field, is that it may be used for several applications in our Solar system and especially in the Earth-Moon one (see [1,2], among others) while another asset is that it may be utilized for revealing the orbital dynamics of a natural/artificial satellite or an asteroid (e.g., [3][4][5]). In addition, it may be adopted for the consideration of the dynamical stability of an inner terrestrial planet which is in mean-motion resonance with an outer giant planet [6]. For the study of the restricted three-body problem, several authors concentrated their work on periodic orbits because they play a key role on the exploration of its dynamics due to their immediate connection with the characterization of nearby orbits (see [7][8][9][10][11] and references therein). Additionally, stable periodic orbits are important in planetary dynamics since they can host real planetary systems [12]. These orbits are also strongly connected with the librational motion either in two or three-dimensions. For example, Voyatzis et al. [13] determined resonant families of three-dimensional periodic orbits related to the dynamics of the Sun-Neptune and a trans-Neptunian object system in order to study the librations and the long-term evolution of orbits near them. Furthermore, several modifications of the restricted problem have been proposed in the past so as to make it more realistic to certain systems of Celestial Mechanics. Such modifications include the radiation and oblateness effects of the primaries [14][15][16][17] or the incorporation of some relativistic terms [18,19], while various works involve also a larger number for the primary bodies [20,21].
A special variant of the classical restricted three-body problem is the Hill one which was firstly proposed by Hill [22] for the study of the moon's motion. In this limiting case of the restricted problem the massless body is attracted by two primary bodies one of which is extremely larger from the second one, e.g. the Sun and the Earth. In the rotating coordinate system, the smaller primary is located at the origin while the positive Ox−axis points to the larger one which is always at infinite distance from the secondary. For this problem, Hénon [23] explored numerically the network of families of simple planar periodic orbits together with their horizontal stability properties while the same author also studied the vertical stability of these families [24]. The three-dimensional periodic solutions emanating from the two collinear equilibrium points of the Hill's problem were determined by Zagouras and Markellos [25]. For the same problem, Hénon [26] searched for families of multiple planar periodic orbits while low-energy escaping trajectories were numerically investigated by means of the Poincaré maps by Villac and Scheeres [27].
In the framework of the Hill's problem where the radiation of the larger primary is also taken into account Kanavos et al. [28] studied its equilibrium points and produced a general map of symmetric periodic orbits in the space of initial conditions while Papadakis [29] presented families of symmetric periodic orbits in the case of regularized variables through the Levi-Civita coordinate transformation. In addition, for the same modified model, de Bustos et al. [30] worked on the bifurcation analysis of the main families of simple periodic orbits. Also, Markellos et al. [31] introduced the version of Hill's problem with oblateness for which they determined the Hill stability of direct orbits around the smaller primary while the network of families of simple and multiple planar periodic orbits was computed by Perdiou et al. [32] and Perdiou [33], respectively. For the same problem, Papadakis [34] determined the map of the basic families using the regularized equations of motion. In the case where both the oblateness and radiation are incorporated in the model Perdiou et al. [35] presented the chart of the basic planar families together with their stability giving special attention to the stability of retrograde satellites whilst by the use of several numerical techniques Zotos [36] revealed the fractal basins of attraction associated to the collinear equilibria in the complex plane .
Furthermore, in the framework of the classical spatial Hill three-body problem, Batkhin and Batkhina [37] investigated the families of spatial periodic orbits which bifurcate from the Vertical Critical (VC) orbits of the basic families of simple planar periodic solutions and formed a network that connect those planar orbits with the determined spatial ones. In the same vein, our aim here is to numerically explore all the families of three-dimensional periodic orbits (up to their natural termination), which emerge through their bifurcations from the Vertical Self-Resonant (VSR) orbits of the previous mentioned basic planar families. In all the considered cases, we find that each VSR orbit gives rise to two branches of families of three-dimensional periodic orbits whose type of symmetry depends on their own multiplicity. In particular, if the multiplicity of the detected spatial orbits is odd, the member orbits of the one branch is of axisymmetric type while the members of the second one possess the plane symmetric type. In the case where the spotted spatial orbits have been ascertained to have even multiplicity, both the generated branches are constituted by doubly symmetric periodic orbits. This pattern was firstly observed by Robin and Markellos [38] for the vertical branches of the basic family of retrograde satellites in the circular restricted three-body problem. Our results focus to the VSR orbits which give rise to three-dimensional periodic orbits with multiplicity three and four which means that the generated spatial families consist of orbits having the triple or quadruple the period comparably to that of the VSR orbit, respectively. Twenty four such families are found, six of which consist of axisymmetric orbits, six of plane symmetric ones and twelve families are composed by doubly symmetric periodic orbits.
Our work is structured as follows. In Section 2, we recall the equations of motion of the Hill's problem and compute the families of simple planar symmetric periodic orbits together with their stability properties. Specifically, we determine accurately the VSR periodic orbits at which families of three-dimensional periodic orbits bifurcate where their orbits have the triple or quadruple multiplicity (period) with respect to the multiplicity of the detected planar VSR orbits. In Section 3, the families of spatial periodic orbits that bifurcate from the planar VSR orbits are computed and presented. Finally, in Section 4, we summarize our work and conclude.

Planar Motion
In this section we study the simple planar symmetric periodic orbits, together with their stability properties, of the planar Hill problem. In particular, we recompute the network of its basic families, as it was determined by [23], in order to detect and determine accurately the planar symmetric periodic orbits at which families of three-dimensional periodic orbits bifurcate having period up to the quadruple of that of the planar orbits.

Hill's Equations of Motion
The equations of motion of the Hill's problem are (for a detailed description of the problem we may refer here to the famous book by Szebehely [39]): where r = x 2 + y 2 + z 2 is the distance of the particle from the secondary while the relevant potential function has the form: The problem possesses the following integral: where Γ is the Jacobi-like constant related to the Jacobi constant C of the restricted three-body problem by the relation C = 3 + µ 2/3 Γ. We recall that both Ox, Oy are axes of symmetry and that two collinear equilibrium points exist; L 1 which is located on the negative axis and L 2 on the positive one.

Vertical Stability
To study the stability of motion of the massless body we introduce the new variables x 1 = x, x 2 = y, x 3 = z, x 4 =ẋ, x 5 =ẏ, x 6 =ż and write equations (1) in the following form: where The coordinates of the third body, along any solution, depend uniquely on the initial conditions and time, namely x i = (x 01 , x 02 , x 03 , x 04 , x 05 , x 06 ; t), i = 1, 2, . . . , 6, and their partial derivatives with respect to the initial conditions fulfill the variational equations: or equivalently by using matrix notation: while the involved partial derivatives in P can be easily computed. The variational equations take the following special form: where we have abbreviated v ij = ∂x i /∂x 0j and f ij = ∂ f i /∂x 0j . Of importance here are the variations which describe the stability of planar periodic orbits when small perturbations perpendicular to the plane of motion occur. These variations correspond to the vertical stability indices as they were defined by [40] and provide information for families of three-dimensional periodic orbits bifurcating from the plane. Explicitly, the corresponding indices are: and stability occurs when |s v | < 1, where s v = (a v + d v )/2, while for planar symmetric periodic orbits it holds a v = d v and the latter condition reduces to: According to Hénon [40] when a v = 1 the basic family of planar periodic orbits bifurcates with a family of three-dimensional periodic orbits of the same period while the case a v = −1 matches to a corresponding period doubling bifurcation, i.e. the three-dimensional periodic orbits of the bifurcating family have the double period w.r.t. that of the planar periodic orbits of the basic one. The respective planar periodic orbits are known as VC. Condition (11) also indicates the existence of higher order resonances, i.e. families of three-dimensional periodic orbits are generated from the plane whose orbits have multiple period with respect to the period of the planar periodic orbits of the original family. These are the VSR orbits and correspond to a value of the vertical stability index: where p, q ∈ Z with p < q and p/q irreducible fraction while q is the multiplicity of the bifurcating spatial family with q = 1, 2 since these values correspond to the critical cases described above. So, if T is the period of the VSR orbit, the bifurcating spatial periodic orbit will have the period qT.

Families of Planar Periodic Orbits
In order to detect appropriate initial conditions for the accurate computation of the basic families of the Hill's problem we use the grid search method as it was described by Markellos et al. [41]. This method is appropriate for the detection of planar symmetric periodic orbits and has been used by several authors in order to sketch the skeleton of the basic families of periodic orbits in several dynamical models of two degrees of freedom (see, e.g., [42][43][44], among others).
For its application, we briefly recall that, we scan the plane of initial conditions (Γ 0 , x 0 ) so as to obtain a thin grid. Since we look for symmetric periodic orbits each node possesses initial conditions of the form (x, y, z,ẋ,ẏ,ż) = (x 0 , 0, 0, 0,ẏ 0 , 0), whereẏ 0 (> 0) is determined from the Jacobi equation (3) using the specific values (Γ 0 , x 0 ) of the node. The equations of motion (1) are integrated numerically up to the k−th intersection of the orbit with the Ox−axis. Then, we move to the next node of the grid which corresponds, e.g., to the same value of the Jacobi constant Γ 0 and a slightly different value x 0 + δx 0 of x, and using these new initial conditions we integrate the equations of motion again up to the k−th intersection of this orbit with the Ox−axis and look for a change of sign of the value ofẋ at this intersection. If such a change is found a symmetric periodic orbit, having 2k− intersections in total with the Ox−axis, has been detected with an initial value of x between the two values of x 0 and x 0 + δx 0 for which the change of sign was spotted. The second perpendicular crossing of an orbit with the Ox−axis ensures that the orbit will be a closed curve in the phase space according to the mirror theorem [45]. The same procedure is followed for all nodes of the grid.
The resulting, by the grid search technique, points on the (Γ 0 , x 0 ) plane correspond to roughly determined initial conditions which are then used for the accurate computation of planar symmetric periodic orbits. According to the above mentioned discussion, a symmetric orbit will be periodic if at the half number k− of its total crossings 2k− with the Ox−axis (or equivalently at the half period) it fulfills the following periodicity condition: If this condition is not satisfied we seek appropriate corrections δx 0 and δẏ 0 of the initial conditions so as to fulfill it. Linearizing (13) in the corrections we obtain: from which by fixing one of them, e.g. δx 0 , we obtain the remaining one in the form: When a symmetric periodic orbit has been determined with initial conditions (x 0 , 0, 0, 0,ẏ 0 , 0), this means that the RHS of (14) is equal to zero, we modify ∆x 0 to an arbitrarily chosen small value ε, e.g. (14) and get in this way the prediction for a next periodic orbit in the form: where now in expansion (14) we have set δx 0 ≡ ∆x 0 and δẏ 0 ≡ ∆ẏ 0 to distinguish between prediction and correction steps. The iterative application of the above procedure represents the predictor-corrector algorithm for the accurate computation of families of planar symmetric periodic orbits.  By applying the grid search and the aforementioned predictor-corrector techniques we have computed the basic families of planar symmetric periodic orbits of the classical Hill problem. The resulting family characteristics in the plane of initial conditions (Γ 0 , x 0 ) are shown in Figure 1. We note here that, for the Hill problem the respective families have been firstly computed and presented by Hénon [23] but we recompute them here in order to detect appropriate initial conditions for the VSR orbits (which had not been determined in that paper) which will be used as seed for their accurate computation. Also, for each computed family we keep the name which had been given by Hénon.
We now deal with the VC as well as the VSR orbits at which families of three-dimensional periodic orbits bifurcate. Our aim is to compute them with a predetermined accuracy, e.g. 0.5 × 10 −8 , and use them as appropriate starting points for the determination of the corresponding bifurcating families of spatial periodic orbits. So, at the half number of the total crossings of a planar symmetric orbit with the Ox−axis, i.e. at its half period t = T/2, the orbit will be considered to be periodic and VC or VSR if it simultaneously fulfills the following conditions: where d may take the values −1, −1/2, 0 and 1, which correspond to the specific bifurcations of the family of planar periodic orbits with a family of spatial periodic orbits which initially have the double, triple, quadruple or the same period with that of the planar one, respectively. These particular values of d are obtained directly from the relevant relations (11) and (12). If the above conditions do not hold we look for appropriate corrections δx 0 and δẏ 0 of the initial conditions such that to fulfill (17). Expanding the arising equations in Taylor series and keeping only the linear terms we get: from which we obtain the corrector system. In Figure 2 we present the vertical stability diagrams of all the computed families. In particular, in Figure 2(a) we show the vertical stability of the Lyapunov family a emanating from the collinear equilibrium point L 2 . Due to the symmetry of the classical Hill problem with respect to both axes the corresponding stability diagram of the second Lyapunov family c emanating from the other collinear equilibrium point L 1 is the same; therefore it is not presented here separately. As we see, there are three VC orbits, denoted by a1 v , a2 v and a5 v , from which two families of 3D periodic orbits of the same period and one family of such orbits of double period bifurcate. Also, one VSR orbit exists (denoted by a4 v ) from which a family of spatial periodic orbits of triple period bifurcates and one VSR orbit (a3 v ) from which a family of 3D periodic orbits of quadruple period bifurcates.
In Figure 2(b) the vertical stability diagram of family g is shown. For this family we observe that there are six VC and VSR periodic orbits. Specifically, at the VC orbit g3 v a family of 3D periodic orbits of double period bifurcate while at g6 v the 3D members of the bifurcated family are of the same period. Note here that, the stability curve of this family is tangential to the line a v = −1 at its VC orbit g3 v providing that this VC orbit is a double root. In addition, in this frame we see the existence of two VSR orbits (g2 v and g4 v ) at which two families of spatial periodic orbits of triple period bifurcate as well as two VSR orbits (g1 v and g5 v ) at which families of spatial periodic orbits of the quadruple period bifurcate.
Finally, in Figure 2(c), the vertical stability of family g is depicted. In this figure, it can be seen that there are five VC and six VSR periodic orbits. Two VC orbits (g 2 7 v and g 2 8 v ) give rise to two families of spatial periodic orbits of the same period while three VC orbits (g 3 v , g 4 v and g 2 11 v ) give rise to three such families whose members have the double period. Additionally, from the VSR orbits g 2 v , g 5 v and g 2 10 v three families of spatial periodic orbits of triple period bifurcate while from the VSR orbits g 1 v , g 6 v and g 2 9 v three families of such orbits bifurcate whose member orbits have the quadruple period. We note also here that, all the VC orbits had been determined in [24] but none of the VSR orbits and we choose to present them also here for the completeness of our study as well as for reader's convenience. Additionally, we note here, that family f does not possess any VSR orbits with q = 3 or q = 4 in relation (12), therefore we do not show its vertical stability diagram.
Accurate initial conditions of the determined VC and VSR orbits of families a, g and g are given in Table 1. In particular, in this table, we give the half of the orbit's period T/2, the position and velocity components x 0 andẏ 0 at the initial vertical intersection of the VC or VSR orbit with the Ox−axis (this means that the remaining components of the position and velocity vectors are y 0 = 0 andẋ 0 = 0), respectively, as well as the value of the Jacobi constant Γ 0 as it is computed through equation (3) of the integral of motion for (x 0 , 0, 0,ẏ 0 ). Finally, we present the value of the position x cut at the second orbit's vertical intersection with the Ox−axis while in the last three columns of the table we incorporate the vertical stability indices a v , b v , c v and d v , as they are determined by variations (10) and their values define the kind of the vertical bifurcation that occurs at a specific VSR planar periodic orbit. Table 1. VC and VSR orbits of families a, g and g . Due to the symmetry w.r.t. both axes the corresponding orbits of the Lyapunov family c emanating from the other equilibrium point L 1 are obtained directly from the Lyapunov family a; so, they are not given here.

Orbit
T

Spatial Periodic Orbits
The VSR orbits give rise to families of three-dimensional periodic orbits which may possess all possible types of symmetry. In particular, the multiplicity q, as it is defined in (12), determines the symmetry properties of the generated spatial periodic orbits. These orbits possess, at least initially, the period qT, where T is the whole period of a VSR orbit. Additionally, at each VSR orbit two branches of spatial periodic orbits bifurcate, i.e. from VSR orbits families always occur in pairs. This has been firstly identified by Robin and Markellos [38] who also described the mechanism where the spatial periodic orbits branch out from the plane. They pointed out that at a VSR orbit, i.e. q ≥ 3, exactly two branches of three-dimensional periodic orbits bifurcate while their symmetry properties depend on their own multiplicity. Specifically, in case that the generated spatial periodic orbits are of: 1. Odd multiplicity q, i.e. their period is (2n + 1)T, q = 2n + 1, n = 1, 2, . . . , with T being the VSR orbit's period, two such families branch out from the corresponding VSR orbit; one family consists of axis-axis symmetric periodic orbits while the members of the other family are plane-plane symmetric. More precisely, each branch has one of the following symmetries: (a) Ox−Ox symmetry, i.e. they are axisymmetric spatial periodic orbits. They start on the Ox−axis with initial conditions of the form (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ) while at their half period T/2 return on this axis with final conditions of the form (x cut , 0, 0, 0,ẏ cut ,ż cut ). (b) Oxz−Oxz symmetry, i.e. they are plane symmetric spatial periodic orbits. They start on the Oxz−plane with initial conditions of the form (x 0 , 0, z 0 , 0,ẏ 0 , 0) while at their half period T/2 return on this plane perpendicularly with final conditions (x cut , 0, z cut , 0,ẏ cut , 0).
2. Even multiplicity q, i.e. their period is (2n)T, q = 2n, n = 2, 3, . . . , with T being the VSR orbit's period, two families branch out from the corresponding VSR orbit; both of them consist of orbits which are doubly symmetric according to the following characteristics: (a) Ox−Oxz symmetry, i.e. they are doubly symmetric spatial periodic orbits. They start on the Ox−axis with initial conditions of the form (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ) while at the quarter of their period T/4 return perpendicularly on the Oxz−plane with the final conditions (x cut , 0, z cut , 0,ẏ cut , 0).
(b) Oxz−Ox symmetry, i.e. they are also doubly symmetric spatial periodic orbits. However, now they start on the Oxz−plane with initial conditions of the form (x 0 , 0, z 0 , 0,ẏ 0 , 0) while at the quarter of their period T/4 return on the Ox−axis with the final conditions (x cut , 0, 0, 0,ẏ cut ,ż cut ), respectively.
For the computation of these branches we can construct corrector-predictor algorithms based on the periodicity conditions. So, to compute a three-dimensional periodic orbit of type, e.g., Ox−Oxz double symmetry (case 2(a) above) we choose an initial state vector of the form x 0 = (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ). In this case, at the quarter T/4 of the period, namely the orbit meets the Oxz−plane, the following conditions must be hold:ẋ cut (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ) = 0, z cut (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ) = 0, (19) i.e., the orbit has to cross perpendicularly the Oxz−plane in order the periodicity to be established. In general, these are not satisfied, so we seek appropriate corrections δx 0 , δẏ 0 and δż 0 of the initial conditions. Then, the resulting equations are expanded in Taylor series up to first order terms obtaining: Since we have two simultaneous equations with three unknowns, we choose to fix one of them, e.g.
δż 0 = 0, and solve for obtaining the remaining corrections. In this case, these are: where the partial derivatives v 41 = ∂ẋ cut /∂x 0 , v 45 = ∂ẋ cut /∂ẏ 0 , v 61 = ∂ż cut /∂x 0 and v 65 = ∂ż cut /∂ẏ 0 are elements of the variational matrix given in (8). By iterating this process the three-dimensional periodic solution will be computed with the desired accuracy. If a spatial periodic orbit has been sought, i.e. conditions (19) are fulfilled with a predetermined accuracy, a next such orbit existing in its neighbourhood can be predicted. To do so, we use the linearized system (20) and slightly change ∆ż 0 to an arbitrarily chosen small constant ε in order to get the remaining predictions in the form: where in system (20) we have now set δx 0 ≡ ∆x 0 , δẏ 0 ≡ ∆ẏ 0 and δż 0 ≡ ∆ż 0 to distinguish the prediction and correction steps. Also, the new terms appeared in (22) are v 46 = ∂ẋ cut /∂ż 0 and v 66 = ∂ż cut /∂ż 0 . The stability of a three-dimensional periodic orbit can be determined through the variational equations (7) and using the following formulas [38]: with α = 2 − Tr(V) and β = (α 2 + 2 − Tr(V) 2 )/2, while V is the variational matrix given in (8), computed at the orbit's period T. In this case, the periodic orbit in question is stable if the defined parameters P, Q of (23) are reals for which it simultaneously holds |P| < 2 and |Q| < 2. Note that, we can also exploit the symmetry properties of the periodic orbit so as to determine the variational matrix at the half or the quarter of the period T, i.e. T 1 = T/2 or T 2 = T/4, depending on whether it possesses a simple or double symmetry, respectively. In particular, for the simple types of symmetry we get the variational matrix by using the following two formulas [38]: where the first one corresponds to the Ox−Ox type of symmetry and the second holds for the Oxz−Oxz plane one. Finally, for the orbits of double symmetry we have that: The first formula fits to the Ox−Oxz double symmetry while the second to the Oxz−Ox one. In all the above mentioned cases, L, M are the diagonal matrices L = diag{1, −1, 1, −1, 1, −1} and M = diag{1, −1, −1, −1, 1, 1}, respectively. Obviously, the determination of the variational matrix through transformations (24) and (25) offers economy to the computations of the numerical integration. So, by applying the above mentioned techniques we have determined twenty four families of three-dimensional periodic orbits, together with their stability properties, which bifurcate from the VSR orbits given in Table 1. Note that, the corresponding families where the VC orbits, of the same table, generate have been discussed by Batkhin and Batkhina [37], so we do not consider their study here. Specifically, Table 1 contains twelve VSR periodic orbits; each one of these planar orbits possesses two vertical intersections with the Ox−axis (given in this table as x 0 and x cut , respectively) so as a mirror configuration to be satisfied [45]. Twelve families of three-dimensional periodic orbits branch from the first vertical intersection of a VSR orbit while the rest twelve spatial families are generated from the second intersection of the same VSR orbit. All these branches of spatial periodic orbits have been computed and followed to their natural end.
The planar Lyapunov family a has two VSR orbits under consideration at which four families of spatial periodic orbits bifurcate. These are shown in Figure 3 where we have plotted their characteristic curves in the space of their initial conditions. Families f , correspondingly. In the superscript, the first number follows the running number of the corresponding VSR orbit while the second one indicates that the period of the spatial periodic orbits is, at least initially, qT, q = 3 or 4, where T is the period of the VSR orbit, i.e. it describes the period commensurability between the planar and spatial orbits.  More precisely, in the first frame of Figure 3 we present the characteristic curve of the planar family a by giving its initial conditions (Γ 0 , x 0 ), where the value of the Jacobi constant Γ 0 is determined by (3) for x 0 , y 0 = 0, z 0 = 0,ẋ 0 = 0,ẏ 0 andż 0 = 0 (first perpendicular cuts of family's orbits with the Ox−axis), while the second frame depicts the respective characteristic curve in the (Γ 0 , x cut ) plane, i.e. these conditions correspond to the second perpendicular cuts x cut , y cut = 0, z cut = 0,ẋ cut = 0,ẏ cut andż cut = 0 of the planar orbits with the Ox−axis (obviously it holds Γ 0 ≡ Γ c ). Also, according to the multiplicity q, family f have spatial members which are doubly symmetric. All the computed families go to three-dimensional collision orbits, so we consider that this is their natural termination and stop calculating them.
In Figure 4 the corresponding spatial branches which the planar family g generates are presented. This family possesses four VSR orbits therefore, eight such branches bifurcate. For their names we follow the same terminology described in the previous paragraph. So, g1 v gives rise to families All families of three-dimensional periodic orbits which bifurcate from the VSR orbits of the planar family g and their members have initially multiplicities 3 and 4 are presented in Figure 5. The six VSR orbits g i v , i = 1, 2, 5, 6, 9 and 10, of g which have been spotted indicate that, in total, twelve vertical branches intersect them. In particular, the branched families f g 2 , which are generated from the first vertical intersections of the VSR orbits g 2 v , g 5 v and g 2 10 v , respectively, consist of three-dimensional members which admit the Ox−Ox axis symmetry. Moreover, the first two families terminate on the Oxy−plane with coplanar orbits while the third one goes to a collision orbit in three-dimensions. The second vertical intersections of these VSR orbits give rise to the families f and f (9,4) g 2 , emerge from the VSR orbits g 1 v , g 6 v and g 2 9 v , respectively, and they are constituted by orbits of the double symmetry Ox−Oxz. The first family terminate on the plane while the rest two go to collision orbits. Finally, the same VSR orbits give also rise to the branches f , respectively, where the first two end on the plane while the third one goes, in its evolution, to collision orbits. The member orbits of the latter families are Oxz−Ox doubly symmetric.
The terminations of the spatial families which were found to end on the physical plane Oxy are shown in Figure 6 with the cyan, green and magenta coloured dots. In this figure, we have plotted the end up also on planar families with members of multiplicity three (the magenta dots which are shown to be located on the characteristic curves of g and g 2 ). The latter four cases of planar families have not been identified in our study since we have only considered the basic families of the Hill's problem. Also, the two cyan dots on the characteristic curve of family f withẏ < 0, i.e. family f cut , correspond to the termination orbits of families f . However, the first family bifurcates from a VSR orbit of the family f cut of retrograde satellites with q = 5 in relation (12) while the second one emanates from a VSR orbit with q = 6, two cases which have not also been considered in our study. The remaining five cyan dots located on the characteristic curves of g , g cut and g cut represent the terminations of families f 4) g and f (1cut,4) g which return on VSR orbits of these planar families (existing in Table 1) or on VSR orbits which are images of the latter ones under reflection in the origin and have the same value of the Jacobi constant Γ. In particular, families f originates, so these two families are essentially the same. Table 2 incorporates initial conditions for sample members of the six computed spatial families whose orbits are axisymmetric and have been generated from the VSR orbits a4v, g2v, g4v, g 2v, g 5v and g 2 10v. In particular, each entry involves the orbit's half period T/2, the position and velocity components (x 0 , 0, 0, 0,ẏ 0 ,ż 0 ) which the orbit has initially on the Ox−axis as well as the value of the Jacobi constant. The last column indicates whether the family contains some stable parts (S) or not (U). In Table 3 we present sample members of the six families which bifurcate from the other vertical intersection of the same planar VSR orbits. Since their orbits are planar symmetric, we give now the position and velocity components (x 0 , 0, z 0 , 0,ẏ 0 , 0) when the orbit starts perpendicularly from the Oxz−plane. Finally, in the Tables 4 and 5 we provide data for sample orbits of the families which the VSR orbits a3v, g1v, g5v, g 1v, g 6v and g 2 9v generate from their two vertical intersections with the Ox−axis. The presented data correspond to the quarter of the orbits' period T/4 since these families consist of spatial orbits which are doubly symmetric. Besides, the multiplicity of the generated spatial orbit's period was noticed to specify the symmetry properties of the orbit. Our results showed that the spatial orbits of the computed branches may have the following types of symmetry: axis-axis, plane-plane as well as double symmetry which combines both the other two types. The numerical continuation of these branches were accomplished through predictor-corrector algorithms which were based on mirror configurations where the periodicity conditions fulfill according to the spatial orbit's type of symmetry. Thirteen branches were found to terminate with coplanar periodic orbits while the remaining eleven ones were collided in their evolution with the smaller primary body. Additionally, among the twenty four computed families of three-dimensional periodic orbits, we found that only five of them include stable parts; these are Since the vertical stability diagrams of the basic planar families indicated the existence of higher order resonances (than q = 4 which was the maximum considered here) for the generation of vertical branches, a natural continuation of the present work would be to consider a similar study for their determination. Also, another interesting extension could be to incorporate additional forces in the Hill's problem, such as the radiation of the larger primary body, i.e. the Sun emits radiation, and examine the influence of these forces to the considered here spatial solutions.
Funding: This research received no external funding.