{'fraction_non_alphanumeric': 0.030915601890110823, 'fraction_numerical': 0.01317992075095883, 'mean_word_length': 4.409054441260745, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orbit classification in the Hill problem: I. The classical case

The case of the classical Hill problem is numerically investigated by performing a thorough and systematic classification of the initial conditions of the orbits. More precisely, the initial conditions of the orbits are classified into four categories: (i) non-escaping regular orbits; (ii) trapped chaotic orbits; (iii) escaping orbits; and (iv) collision orbits. In order to obtain a more general and complete view of the orbital structure of the dynamical system our exploration takes place in both planar (2D) and the spatial (3D) version of the Hill problem. For the 2D system we numerically integrate large sets of initial conditions in several types of planes, while for the system with three degrees of freedom, three-dimensional distributions of initial conditions of orbits are examined. For distinguishing between ordered and chaotic bounded motion the Smaller ALingment Index (SALI) method is used. We managed to locate the several bounded basins, as well as the basins of escape and collision and also to relate them with the corresponding escape and collision time of the orbits. Our numerical calculations indicate that the overall orbital dynamics of the Hamiltonian system is a complicated but highly interested problem. We hope our contribution to be useful for a further understanding of the orbital properties of the classical Hill problem.

nonlinear dynamics (e.g., [8-10, 42, 43]). In the case where the energy of escape has a finite value escape is possible. In particular, when the particles have values of energy higher than the critical energy of escape they are able to find the openings, or exit channels, in the equipotential surface and therefore escape to infinity. Systems with exit channels are also known as open or leaking Hamiltonian systems. The literature is replete with studies on such open Hamiltonian systems (e.g., [4, 12, 17, 30, 31, 35, 50-52, 60, 61, 63, 70]).
The topic of leaking Hamiltonian systems however is by far less explored than the closely related issue of chaotic scattering. In many previous works the chaos theory has been successfully used in order to investigate and explain the phenomenon of chaotic scattering (e.g., [5-7, 14-16, 18, 21-28, 44-48]). Needles to say that all the above-mentioned citations on both issues of chaotic scattering and open Hamiltonian systems are exemplary rather than exhaustive.
The restricted three-body problem (RTBP for short) is an excellent example of an open Hamiltonian system with escapes (e.g., [57,58]). Over the last decade or so a large number of studies have been devoted on orbit classification in the RTBP. It all started with the pioneer works of Nagler [33,34] where initial conditions of orbits were classified as bounded, escaping or collision. Moreover, bounded orbits were further classified into orbital families by taking into account the type of the motion of the test particle around the primary bodies. In the same vein, orbit classification in the RTBP with perturbations has also been performed in [64,65] where we investigated the influence of the oblateness, while in [66] we explored how the radiation pressure affects the orbital content of the dynamical system.
Similarly, orbit classification has also been conducted in celestial mechanics and dynamical astronomy regarding planetary systems. In particular, the orbital dynamics of the Earth-Moon system with two degrees of freedom using a scattering region around the Moon has been studied in [13].
In [69] the numerical investigation has been expanded into three dimensions thus exploring the orbital structure of the three degrees of freedom Earth-Moon system. Furthermore, very recently the escape and collision dynamics in the Saturn-Titan and in the Pluto-Charon binary planetary systems have been revealed in [67] and [68], respectively, by classifying initial conditions of orbits in the configuration space.
The Hill limiting case is in fact a simplified modification of the RTBP which focus on the vicinity of the secondary (e.g., [20,38,39]). This allow us to study the motion of the test particles in the neighborhood of the equilibrium points L 1 and L 2 . At this point it should be emphasized that the Hill approximation is valid only when the mass of the secondary is much smaller compared with the mass of the primary body. One can directly obtain the Hill model from the classical RTBP by translating the origin to the center of the secondary and also by rescaling the coordinates by a factor µ 1/3 . In [49] the authors performed an extensive investigation of the planar Hill problem by computing several homoclinic and heteroclinic connections. In this work, which is the first part of a series of papers, we shall perform a thorough and systematic orbit classification in the classical Hill problem. To our knowledge, there are no previous detailed and systematic numerical studies regarding orbit classification in the classical Hill problem. On this basis, the new information presented in this paper is exactly the contribution as well as the novelty of this research work. In the following papers of the series we will explore how several perturbations (i.e., the oblateness and the radiation pressure) influence the orbital structure of the dynamical system.
The present paper is organized as follows: In Section 2 we present in detail the derivation and the basic dynamical properties of the mathematical model. All the computational techniques we used in order to classify the initial conditions of the orbits are described in Section 3. In the following Section, a thorough and systematic numerical investigation takes place which allow us to reveal the overall orbital structure of the Hill system with two and three degrees of freedom. Our paper ends with Section 6 where the discussion and the main conclusions of our research are given.

Derivation of the mathematical model
The classical Hill problem is derived from the circular restricted three-body problem if we make some scale changes and if we also take the limiting case where the mass ratio µ = m 2 /(m 1 + m 2 ) tends to zero (µ → 0).
The time-independent effective potential function of the circular restricted three-body problem is where (X, Y, Z) are the coordinates of the test particle, while are the distances of the test particle from the two primaries. According to [54], in a rotating system of reference the equations of motion are , We now place the origin of the coordinates at the center of the secondary and we change the scale of lengths by a factor of µ 1/3 Here we would like to note that the notations (X, Y, Z) and (x, y, z), regarding the coordinates for the RTBP and the Hill problem, respectively are according to [37]. Applying the above-mentioned transformations to (1) we obtain with r = x 2 + y 2 + z 2 .
Taking the limit of the right-hand side of (5) for µ → 0 we arrive at the potential function of the classical Hill problem More details about the derivation of the potential function of the Hill problem are given in the Appendix.
The corresponding equations of motion are The Hill problem admits a Jacobi integral of motion where J is the new Jacobi constant which is related to the Jacobi constant C of the restricted three-body problem through the relation In order to locate the positions of the equilibrium points of the system we have to solve the system of differential equations The Hill problem has only two equilibrium points which are located on the x axis at x L = ±1/3 1/3 (see Fig. 1a). The value of the Jacobi integral at the equilibrium points is J L = 3 4/3 and this is a critical value. This is because for J < J L the zero velocity surfaces are open and two symmetrical channels (exits) are present near the equilibrium points L 1 and L 2 through which the test particles can escape from the interior region (−x L ≤ x ≤ x L ) of the system. Therefore we may argue that J L plays, in a way, the role of the energy of escape and this is exactly why it is a critical value of the Jacobi integral. In Fig. 1b we present a plot of the three dimensional equipotential surface 2W(x, y, z) = J, when J = 4.25.
We observe the presence of the two channels (or throats) through which the test particles can leak out. One could say that these two exit channels act, in a way, as hoses thus allowing the test particles to travel from the interior region of the system to the "outside world". Channel 1 indicates escape towards the negative x direction, while channel 2 indicates escape towards the positive x direction.

Computational methods
In order to reveal the orbital structure of the Hill problem we need to define sets of initial conditions whose properties will be identified. For obtaining a more complete and spherical view regarding the nature of orbits we are going to explore both the two-dimensional (2D) and the three-dimensional (3D) Hill systems. In this work we decided to classify orbits with initial conditions in the configuration space only. In particular, for several values of the Jacobi integral J we define dense uniform grids of initial conditions inside the corresponding zero velocity surface. For the 2D system (z = 0), grids of 1024 × 1024 initial conditions (x 0 , y 0 ) are numerically integrated. For all 2D orbitsẋ 0 = 0, while the initial value ofẏ is always obtained from the Jacobi integral of motion (8), asẏ 0 =ẏ(x 0 , y 0 ,ẋ 0 , J) > 0. In the same vein, for the 3D system, grids of 100 × 100 × 100 initial conditions (x 0 , y 0 , z 0 ) are numerically integrated. For all 3D orbitṡ x 0 =ż 0 = 0, while the initial value ofẏ is always obtained from the energy integral (8), asẏ 0 =ẏ(x 0 , y 0 , z 0 ,ẋ 0 ,ż 0 , J) > 0. In both cases, the initial conditions of the orbits lie inside the interior region (which is the scattering region) with The classification of the initial conditions of the orbits in the Hill problem is a rather demanding task if we take into account that the configuration space extends to infinity. In this study we shall classify initial conditions of orbits into four main categories: 1. Orbits that move in bounded, or non-escaping regular orbits inside the interior region. 2. Orbits that move in trapped chaotic orbits inside the interior region. 3. Orbits that escape from the interior region passing either through L 1 or L 2 . 4. Orbits that collide with the secondary located at the origin of the coordinates.
Our next task is to define appropriate numerical criteria in order to distinguish between the above-mentioned four types of motion. An orbit is considered to escape from the system when x < −x L −δ, or when x > x L +δ, where δ = 0.1. At this point, it should be clarified that the tolerance δ was included so as to avoid the unstable Lyapunov orbits [32], located near the equilibrium points, to be incorrectly classified as escaping orbits. A collision with the secondary occurs when R < R col , where R col = 10 −4 . Our preliminary calculations suggest that a considerable portion of the initial conditions correspond to bounded orbits, which can be either regular or chaotic. Therefore we decided to distinguish between initial conditions of regular non-escaping and trapped chaotic motion. Over the years many dynamical indicators have been introduced for distinguishing between regular and chaotic orbits. We chose to use the Smaller ALingment Index (SALI) method [53], which has been proved a very fast and accurate tool. The mathematical definition of SALI is the following where d − ≡ w 1 (t) − w 2 (t) and d + ≡ w 1 (t) + w 2 (t) are the alignments indices, while w 1 (t) and w 2 (t), are two deviation vectors which initially point in two random directions. For distinguishing between ordered and chaotic motion, all we have to do is to compute the SALI along a time interval t max of numerical integration. In particular, we track simultaneously the time-evolution of the main orbit itself as well as the two deviation vectors w 1 (t) and w 2 (t) in order to compute the SALI. The determination of the nature of an orbit is obtained from the final value of the SALI at the end of the numerical integration. In particular, if SALI > 10 −4 the motion is regular, while if SALI < 10 −8 the motion is chaotic. If the final value of SALI lies in the interval [10 −4 , 10 −8 ] then we have the case of a "sticky" orbit 1 and further numerical integration is required for obtaining the true nature of the orbit. In our study the maximum time of the numerical integration 1 Sticky orbits need an extremely long time interval in order to move away from the invariant sticky tori [36]. Therefore, sticky orbits behave, for long time interval, as regular ones before revealing their true chaotic nature.
(t max ) was set to be equal to 10 4 time units. Orbits that remain inside the interior region after integrating them for t max are considered as non-escaping regular or trapped chaotic, according to their final value of SALI.
The set of the equations of motion (7) as well as the corresponding variational equations, needed for the computation of the SALI, have been forwarded integrated using a double precision Bulirsch-Stoer FORTRAN 77 algorithm (e.g., [40]) with a variable time step. Throughout our numerical calculations the error regarding the conservation of the Jacobi integral of motion (8) was generally smaller than 10 −12 , while in some cases it was observed to be smaller than 10 −14 . For those initial conditions of orbits which move inside a region of radius R < 10 −2 , thus leading to collision with the secondary, the Lemaitre's global regularization method (e.g., [54]) has been applied. All graphical illustrations presented in this paper have been created using the latest version 11 of the software Mathematica (e.g., [59]).

Orbit classification
In this section we shall classify orbits with initial conditions in the configuration space into the four categories analyzed in the previous section. Parallel to the classification we shall also record the time scale (or time period) of the collision and the time scale of the escape. Our aim is to numerically explore the orbital properties of the Hill problem for several values of the Jacobi integral of motion.
In the following we shall present color-coded diagrams (CCDs) in which each pixel is assigned a specific color according to the particular type of the nature of the orbit. These CCDs are modern types of color-coded maps where the phase space is a complex mixture of basins of escape, collision basins and regions of bounded motion. By the term "basin" we refer to a set of initial conditions which lead to a certain final state (collision, escape or bounded motion). Our preliminary numerical calculations indicate that bounded motion is almost always present. Generally speaking, the vast majority of the bounded basins correspond to initial conditions of regular orbits, where an adelphic integral of motion is present. This additional integral poses new restrictions to the available phase space and therefore it prevents them from escaping to infinity.

Results for the 2D system
We begin our investigation with the two-dimensional (2D) system when z =ż = 0. In Fig. 2 we present, using a rich collection of CCDs, the evolution of the orbital structure of the configuration (x, y) plane as the value of the Jacobi constant decreases. The outermost black solid line denotes the zero velocity curve which is defined as 2W(x, y) = J. We Fig. 2 Evolution of the orbital structure of the configuration (x, y) space with decreasing value of the Jacobi integral of motion. The color code is as follows: non-escaping regular orbits (blue), trapped chaotic orbits (yellow), escaping orbits through L 1 (red), escaping orbits through L 2 (green), and collision orbits (cyan). (For the interpretation of references to colour in this figure caption and the corresponding text, the reader is referred to the electronic version of the article.) observe that when J = 4.326, that is an energy level just above the critical value J L , the vast majority of the configuration plane is covered either by stability islands or collision basins. The remaining area is a highly fractal 2 mixture of escaping orbits. Inside this fractal area it is impossible to predict through which channel a test particle will escape. This is true because even a slight change on the initial conditions of an orbit leads the test particle to escape from the opposite channel, which is of course a classical indication of chaotic motion. As the value of the Jacobi integral of motion reduces even further (remember that at the same time the total orbital energy increases) the following important changes take place on the configuration (x, y) plane: -The two symmetrical escape channels near the equilibrium points become more and more wide, as the energetically forbidden regions are reduced. -Several basins of escape start to emerge inside the escape regions, while the fractality of the same regions is reduced. In particular, for J < 4.2 fractal regions are only present in the vicinity of the basin boundaries.
-The area of the basins of escape grows rapidly and for J < 3.7 it covers more than 50% of the (x, y) plane. It is observed that for relatively low values of the Jacobi integral of motion (J < 0) escape basins dominate occupying almost the entire configuration space. -The extent of the main stability island located at the left side of the (x, y) plane (x < 0) is gradually reduced and when J < 0 there is no indication of non-escaping reg-ular motion whatsoever. Additional smaller stability islands appear and disappear with decreasing value of J. -The collision basins are mainly located around the stability island. For J < 2.5 however, they appear only as thin filaments inside the vast escape domains. Our calculations suggest that collision motion remains possible even for extremely low values of the Jacobi value. -With increasing orbital energy it is seen that the vast majority of the escaping orbits leak out through exit channel 2. In fact at the highest energy level studied (J = −0.5) the area covered by initial conditions of orbits that escape through L 2 is almost four times larger than the area of initial conditions of orbits that escape through L 1 . Therefore we may say that in the 2D system exit channel 2 seems to be much more preferable. -Our computations indicate that in the 2D Hill system there is no numerical evidence of trapped chaotic motion.
In the following Fig. 3 we demonstrate how the escape and the collision time of the orbits are distributed on the configuration (x, y) space for the same values of the Jacobi integral shown in Fig. 2. Light reddish colors correspond to fast escape or collision orbits, dark blue/purple colors indicate large escape and collision time, while white color denote regions of non-escaping regular motion and trapped chaotic motion. We note that the scale on the color bar is logarithmic varying always from -2 to 4. Inspecting the distribution of the escape time of orbits it is rather easy to associate the stable manifold of the non-attracting invariant set with medium escape time. Similarly, the largest observed escape rates are directly linked with the presence of sticky motion in the vicinity of boundaries of the several stability islands (see also [11]). It is interesting to emphasize that when J = 4.326 most of the orbits escape only after 1000 time units. However as we proceed to higher energy levels the average escape time is considerably reduced. As it was seen in Fig. 2, for J < 2.5 only a thin filament of collision orbits is presents. Our calculations suggest that these orbits collide within the first ten steps of the numerical integration.
It would be very informative to monitor how the percentages of all types of orbits evolve as a function of the value of the Jacobi integral of motion. Such a diagram is presented in Fig. 4. It was constructed by taking into account the information of the 16 CCDs previously described in Fig. 2. We see that at relatively high values of J (J > 4.2) about half of the configuration (x, y) plane is covered by initial conditions of non-escaping regular orbits. However as the value of J decreases the percentage of non-escaping orbits is also reduced and for J < 0 it completely vanishes. The pattern of the evolution of the rate of collision orbits is very similar. The main difference is that for J < 0 collision motion is still possible even though the corresponding rate is extremely low (about 0.5%). The percentage of escaping orbits which leak out through L 1 initially increases, while for J < 3 it seems to saturate around 30%. On the other hand, the rate of escaping orbits through L 2 increases, almost linearly, and at the highest energy level studied (J = −0.5) they dominate covering about 80% of the (x, y) plane. Taking into consideration the above-mentioned analysis, we may reasonably conclude that in the 2D Hill problem non-escaping regular orbits is the most populated type of orbits at low energy levels, while at high enough values of the energy escaping orbits dominate. In particular we observed that the probability of an orbit to leak out through L 2 is much more higher than through L 1 . Therefore we may claim that exit channel 2 is more preferable in the 2D system.
The CCDs presented in Fig. 2 can provide useful information regarding the orbital structure of the 2D system however, only for some specific values of the Jacobi integral of motion. In order to overcome this limitation we shall adopt the Hénon's method [19] thus using the section y =ẋ = 0, y > 0. This mean that all the initial conditions of the 2D orbits will be launched from the x axis, with x = x 0 , parallel to the y axis. Therefore, only orbits with pericenters on the x axis will be included, while the value of the Jacobi constant J can now be used as an ordinate. Using this method we are able to monitor the evolution of the orbital structure of the 2D Hill system using a continuous spectrum of energy values, rather than a few discrete ones. The orbital structure of the (x, J) plane when J ∈ [−0.5, J L ) is presented in Fig. 5a,  while in Fig. 5b the distribution of the corresponding escape and collision time of the orbits is given.
With a closer look at Fig. 5a we can distinguish a white vertical line located at x = 0. We decided to exclude the initial condition x 0 = 0 for all energy levels because in this case the radius R = x 2 0 + y 2 0 is zero (remember that in the (x, J) plane all orbits have y 0 = 0) and therefore is smaller than the defined collision radius R col = 10 −4 . This mean that all orbits with x 0 = 0 start inside the assumed radius of the secondary and this has no physical meaning. It is seen that at high values of the Jacobi constant, specifically for x < 0, a main stability islands is present. As we proceed to lower values of J however its area is constantly reduced and when J < 0 it disappears. Escape basins on the other hand are always present, while their area grows rapidly as we move away form the critical energy level J L . Basins of escape composed of initial conditions of orbits that escape through L 2 are much more prominent with respect to those corresponding to escape through L 1 . It is interesting to note that Fig. 5a indicates how the fractality of the several basin boundaries varies not only as a function of the Jacobi constant but also of the spatial variable. In particular, we can observe that the fractality of the basin boundaries, which is of course related to the unpredictability, migrates from the upper right side of the (x, J) plane for relatively high values of J (or in other words low values of the total orbital energy) to the lower left side of the same plane for low values of the Jacobi constant. According to Fig. 5b near the vicinity of those fractal basin boundaries the highest values of the escape time of the orbits have been reported.
In Fig. 6 we illustrate how the percentages of all types of orbits with initial conditions on the (x, J) plane of Fig.  5 evolve as a function of the value of the Jacobi integral of motion. We see that the evolution of the rates of all types of orbits are more or less similar to that observed earlier in Fig. 4. The main differences between the two plots, which correspond to initial conditions of orbits in different types of planes, that occur for relatively high energy levels (J < 0) are the following: (i) the rate of escaping orbits through L 1 in the (x, J) plane seems to saturate at a lower value (at about 5%) than in the (x, y) plane which was equal to 20%, (ii) the percentage of escaping orbits through L 2 is more dominant in the (x, J) plane as it occupies more than 95% of the same plane, while in the configuration (x, y) plane it was found to be lower than 80%.

Results for the 3D system
Our numerical investigation continues with the full 3D system. Fig. 7 presents a collection of CCDs using three dimensional distributions of orbits with initial condition inside the (x, y, z) space for the same energy levels of Fig. 2. We observe that in this case the grids composed of the initial conditions of the orbits are in fact three dimensional solids which means that only their outer shell is visible. Three dimensional distributions of initial conditions of orbits have also been recently used in order to reveal the orbital structure of the 3D Earth-Moon system [69]. The corresponding distributions of the escape and collision time of the 3D orbits are given in Fig. 8, where for the initial conditions of the nonescaping regular orbits as well as the trapped chaotic ones we have used transparent white color in order to be able to inspect, in a way, the interior of the solids.
If we want to examine in more detail the interior region of the 3D grids we have to use a tomographic-style approach as in [60]. According to this approach we take the cuts on the (x, y) plane for several levels of the z coordinate. This method however is very costly as for every 3D grid we need several 2D cuts on the (x, y) plane. For saving space we decided to present for every energy level a cut on the (x, z) plane. Therefore all 3D orbits with initial conditions on the (x, z) plane have y 0 =ẋ 0 =ż 0 = 0, while the initial value ofẏ is always obtained from the Jacobi integral of motion (8). There is also one more reason justifying our choice regarding the type of the 3D cuts. The z coordinate is directly related with the 3D motion of orbits, therefore we decided to use a continuous spectrum of z values rather than a few cuts for specific values of z 0 as in [60]. Furthermore, the CCDs of the 2D system presented earlier in Fig. 2 are in fact the z = 0 cuts of the corresponding 3D solids. At this point, we would like to stress out that as far as we know this is the first time that 3D distribution of initial conditions of orbits are systematically classified in the Hill problem.
The evolution of the orbital structure of the (x, z) space, as the value of the Jacobi integral decreases, is revealed through the rich collection of the CCDs presented in Fig. 9, while the corresponding distribution of the escape and collision time of the orbits is given in Fig. 10. The outermost black solid line denotes the zero velocity curve which is defined as 2W(x, y = 0, z) = J. It is seen that the orbital structure is very different with respect to that observed in the 2D system. For high values of the Jacobi constant non-escaping regular motion dominates, while initial conditions of escaping orbits are mainly confined to the right side of the (x, z) planes. A very interesting phenomenon is the presence of trapped chaos. Trapped chaotic orbits have also been observed in other types of 3D open Hamiltonian systems. For instance, in [62] we detected a considerable amount of trapped chaotic orbits when we investigated the orbital properties of an open tidally limited star cluster, using a three dimensional dynamical model. On the other hand, in the 3D Earth-Moon system which was explored in [69] we did not find any numerical evidence of trapped chaos. The most important changes that occur on the (x, z) plane as the value of the Jacobi constant decreases are the following: -The area composed of initial conditions corresponding to non-escaping regular orbits is constantly reduced. It is seen that for J < 2.5 bounded basins are present only at the left side of the (x, z) plane, while for negative values of the Jacobi integral of motion they completely disappear. -Collision basins exist mainly inside the vast regular regions thus creating a complicated mixture of initial conditions of both collision and non-escaping regular orbits. We found that for J < 1.9 there is no indication of collision motion whatsoever. -Initial conditions corresponding to trapped chaotic orbits are visible only for relatively high values of the Jacobi constant, very close to the critical energy value J L . Our computations suggest that for J < 4.2 the amount of trapped chaotic orbits is extremely low and the corresponding initial conditions appear as lonely points randomly scattered in the (x, z) plane and therefore they are not visible. -As we move away from the critical value of the energy to higher levels of energy the basins of escape start to expand. Initially they are mainly located to the right side of the (x, z) plane however as the value of J decreases they migrate also to the left side of the same plane. Finally at extremely low levels of J and particularly for negative values basins of escape cover the entire (x, z) plane. -Our numerical calculations strongly indicate that the two exit channels, even though they are symmetrical, they are not equiprobable. Indeed it was found that the test particles have the tendency to leak out mainly from the exit channel near L 2 . Similar behavior was also observed in the 2D system. It is interesting to note that for extremely low values of the Jacobi integral of motion (J < 0) almost the entire (x, z) plane is dominated by initial conditions of orbits which escape through channel 2. However, a small basin of escape corresponding to channel 1 is still present near the origin.
Exploiting the numerical integration of the three dimensional distributions of initial conditions of orbits shown in Fig. 7 we managed to demonstrate in the diagram of Fig. 11 the evolution of the percentages of all types of orbits as a function of the value of the Jacobi integral of motion. It is seen that at high enough values of J, just above the energy of escape, about 80% of the configuration (x, y, z) space is covered by initial conditions of non-escaping regular orbits. As the value of J decreases however the rate of regular orbits also decreases until it completely disappears for J < 0. The evolution of the rate of collision orbits is quite similar as it starts at about 18% for J > 4.3 and then it drops tending asymptotically to zero for relatively high values of the energy. Here it should be pointed out that additional numerical calculations (not shown here) suggest that collision motion remains mathematically possible, even for extremely low values of the Jacobi constant, even though the corresponding percentage is lower than 0.01%. The percentage of escaping orbits through L 1 initially increases until J = 1.3 where it exhibits its maximum value at around 33.5%, while for lower values of J the tendency is reversed. On the other hand, the rate of orbits that leak out through L 2 is constantly increases. At the highest energy level studied, that is for J = −0.5, escaping orbits through channel 2 cover about 86% of the configuration space, while at the same time orbits that escape through exit channel 1 occupy only 13% of the same space. Thus is becomes more that evident that the escape process in the 3D Hill system is a very complicated phenomenon.
The CCDs presented in both Figs. 7 and 9 reveal the orbital properties of the 3D system however for a fixed value of the Jacobi integral of motion. In order to obtain a more Fig. 9 Evolution of the orbital structure of the (x, z) space with decreasing value of the Jacobi integral of motion. The color code is the same as in Fig. 2. (For the interpretation of references to colour in this figure caption and the corresponding text, the reader is referred to the electronic version of the article.) complete view of the orbital structure of the 3D Hill system, we shall try to gather information of a continuous spectrum of energy values. To achieve this we will adopt the approach used previously for the 2D system. In particular, for specific values of z 0 we define dense uniform grids of initial conditions on the (x, J) plane with y 0 =ẋ 0 =ż 0 = 0, while in all cases the initial value ofẏ is obtained from the energy integral (8). Fig. 12a shows the case when z 0 = 0.012. We observe that the orbital structure is quite simi-lar to that shown in Fig. 5a. This however is anticipated because the particular value of z 0 is relatively low and therefore the motion of the orbits takes place very close to the twodimensional (x, y) plane. When z 0 = 0.2 the pattern changes drastically. The collision basin located at the upper left part of the (x, J) plane seems to weaken, since initial conditions of non-escaping regular orbits emerge inside it. Furthermore it is seen that basins composed of bounded regular orbits expand also at the right side of the plane, while the extent of the basins of initial conditions of orbits that escape through L 1 is reduced. The outermost black solid line is the limiting curve which is defined as In panel (c) where z 0 = 0.4 it is interesting to note the rich fractal boundaries between the several basins of escaping and non-escaping regular orbits. On the contrary in panel (d) where z 0 = 0.6 the degree of fractality is considerable lower since the basin boundaries are much more smooth. In addition, the basins of collision orbits seem to be confined only to the upper left side of the (x, J) plane. It should be clarified that in this case (the 3D system) the vertical white line corresponding to initial conditions of orbits inside the collision radius of the secondary is not present. This is true because all four values of z 0 are such (z 0 > 0.01) so the distance R = x 2 0 + y 2 0 + z 2 0 is always larger than the collision radius, even when x 0 = y 0 = 0. Looking the four Fig. 11 Evolution of the percentages of all types of orbits with initial conditions in the configuration (x, y, z) space, as a function of the Jacobi constant J. (For the interpretation of references to colour in this figure caption and the corresponding text, the reader is referred to the electronic version of the article.) panels of Fig. 12(a-d) it becomes evident that the energetically permissible area on the (x, J) plane is reduced as the initial value of the z coordinate increases. The distribution of the corresponding escape and collision time of the orbits with initial conditions on the (x, J) plane is presented in Fig. 13(a-d). We see in panel (c) that the escape time of the orbits with initial conditions in the vicinity of the fractal basin boundaries are extremely high corresponding to more than 1000 time units. On the other hand, all orbits with initial conditions very close to the equilibrium point L 2 (outer right side of the plots) escape almost immediately within the first couple of time steps of the numerical integration.
We close our numerical investigation by presenting in Fig. 14(a-d) the evolution of the percentages of all types of orbits with initial conditions on the (x, J) plane for the four cases discussed in Fig. 12(a-d). One may observe that in all cases escaping orbits completely dominate for extremely low values of the Jacobi integral of motion, or in other words for high values of the total orbital energy. For high values of J on the other hand, non-escaping regular orbits is the most populated type of orbits in the first three cases, except the case where z 0 = 0.6. In this case escaping orbits dominate but this is actually a numerical artifact created by the strange geometry of the limiting curve (see Fig. 12d).
Here we would like to clarify that the evolution of the percentages of trapped chaotic orbits was not included in all four cases because the corresponding percentage was measured to be always extremely low (lower than 0.1%). Taking into account all the results given in Fig. 14(a-d) we may reasonably conclude that the behavior according to which for high enough values of the energy the probability an orbit to escape through exit channel 2 is considerably much higher than trough L 1 , that initially observed in the 2D system, is still valid also in the 3D system.

Trapped chaos
For 3-dof Hamiltonian systems the general situation for moderate perturbation is the following: There are many surviving KAM tori which fill a part of the energy shell whose measure is larger than zero and even can be large. These invariant surfaces are 3 dimensional and therefore in the 5 dimensional energy shell they do not divide anything. I.e. the complement of all these tori is connected. On the other hand, there is an infinity of resonance zones which contain unstable behaviour. In principle, all these resonance zones (fine chaos strips) are connected and form a single web of resonances (Arnold web). For perturbation different from zero, i.e. for systems which are not integrable, the measure of this web of resonance zones also is larger than zero. It is a fractal web of an infinity of small layers. Accordingly, at very large times an orbit starting in any initial point of this web can reach any other point of the web, even when it usually may take an extremely long time. The motion in the web looks similar to diffusive motion and for moderate perturbation the diffusion constant is very small. The surviving 3 dimensional tori form a kind of fractal of obstacles for the general motion in the space between the tori. Now if we imagine the situation where there is also a small escape channel for escape. This escape channel is in contact with the web of chaotic resonance zones. Therefore, any orbit starting in any of the resonance zones will finally escape after a possibly extremely large time. However, when we observe for a finite time only (this may be a large time, even large compared to the life time of the universe) then we only observe the diffusive like chaotic motion inside of the web and for most orbits we do not see the final escape. For all practical purposes we have the impression having observed trapped chaos.
In the CCDs presented in Fig. 9 we observed that for energy levels just above the critical energy of escape J L there is an amount of trapped chaotic orbits. The phenomenon of trapped chaos has been also observed in other types of Hamiltonian systems (e.g., [29,62]). According to classical chaos theory however, chaotic orbits do not admit a third integral of motion. Therefore for time t → ∞ they will fill all the available phase space inside the zero velocity surface, unless of course they escape. In other words, all trapped chaotic orbits will eventually leak out passing through one of the two escape channels, given enough time of numerical integration. In Fig. 15 we present a characteristic exam- To check whether all trapped chaotic orbits, observed in the corresponding CCD (when J = 4.326), eventually escape or not, we numerically integrated the initial conditions of the 3D trapped chaotic orbits for 10 6 time units. Our re- sults are illustrated in Fig. 16a where one can observe that all these initial conditions, that were initially classified as trapped chaotic ones for 10 4 time units, they have escaped creating a highly fractal escape mixture. The distribution of the corresponding escape time of the 3D orbits is given in panel (b) of Fig. 16, where this time the values on the color bar vary from 4 to 6 (remember that the scale of the accompanying color bar is in logarithmic scale). The same procedure was carried out for all the other sets of trapped chaotic orbits reported in the CCDs for lower values of the Jacobi constant. We found that all of these orbits eventually escape through one of the two channels however after extremely long time of numerical integration (t esc ∈ (10 4 , 10 6 )). Therefore, we may reasonably conclude that in the 3D Hill sys-tem long-lived (or temporarily trapped) chaotic orbits exist for values of energy above but very close to the energy of escape.
In 2-dof systems we have 3 dimensional energy shells and the 2 dimensional KAM tori separate these energy shells in disjoint inner and an outer parts. In general, there are chaotic layers also inside of the KAM tori. But these layers are disconnected from the outside. Therefore, chaos in these inside layers can never reach the exit and consequently it can never escape. Thus it is truly trapped chaos for infinite time in the strict mathematical sense. Here it should be emphasized that none of the examined initial conditions of the 3D trapped chaotic orbits in the CCDs has z 0 = 0. This fact supports our findings of the 2D system, where there was no numerical indication of 2D (temporarily) trapped chaotic motion whatsoever (e.g., [1]).

Discussion and conclusions
The aim of this paper was to reveal the overall orbital structure of the classical Hill problem. By performing a thor-ough and extensive numerical investigation using large sets of initial conditions, in both the two-dimensional and the three-dimensional configuration space, we managed to classify orbits into four categories: (i) non-escaping regular, (ii) trapped chaotic, (iii) escaping, and (iv) collision. We also related the several basins of escape and collision with the corresponding escape and collision time of the orbits. As far as we know, this is the first detailed and systematic nu- merical analysis regarding orbit classification in the classical Hill problem (especially in the 3D configuration space) and therefore this is exactly the contribution as well as the novelty of our study.
A Quad-Core i7 2.4 GHz PC was used for the numerical integration of the sets of the initial conditions of the orbits. For the completion of every set of initial conditions of orbits we needed between about 2 minutes and 4 days of CPU time, depending of course on the collision and escape time of the orbits in each case. Here it should be explained that, for saving time, when an orbit escaped or collided the numerical integration was stopped and proceeded to next available initial condition.
In this paper we provide quantitative information regarding the orbital structure of the classical planar (2D) and spatial (3D Hill problem. The main novel results of our orbit classification can be summarized as follows: 1. At high enough values of the Jacobi integral of motion, or in other words at low values of the total orbital energy non-escaping regular orbits as well as collision orbits dominate in both the 2D and the 3D system, while initial conditions of escaping orbits were found mainly near the equilibrium points. 2. Escaping orbits are the most populated type of orbits at relatively high values of the total orbital energy, where basins of escape were found to occupy almost the entire configuration space. 3. As the value of the Jacobi constant decreases the percentages of both non-escaping regular orbits and collision orbits are reduced. Our computations suggest that for negative values of J there is no indication of bounded motion whatsoever, while on the other hand collision motion is still possible. 4. In the 3D system, and especially for values of energy above yet very close to the critical energy of escape, we identified a portion of trapped chaotic orbits. Additional numerical calculations revealed that all these orbits eventually do escape, while having huge escape periods (more than 10 4 time units). It should be emphasized that the phenomenon of (temporarily) trapped chaos was not observed in the 2D system. 5. In both the 2D and the 3D systems we observed that the two symmetrical escape channels located near the equilibrium points are not equiprobable 3 . In fact, for low values of the Jacobi constant the vast majority of orbits choose to escape through exit channel 2, while the percentage of exit channel 1 remains very low.
For very low values of the mass ratio µ the restricted three-body problem is well approximated by the the Hill problem. It should be emphasized that there are numerous previous works on the field of capture and escape on the RTBP (e.g., [2,41]) as well as on the Hill problem (e.g., [55,56]). Furthermore, similar systematic classification of initial conditions of orbits, with respect to the present one, has been performed in the RTBP in general (e.g., [33,34,[64][65][66]) and also in specific planetary systems (e.g., [13,69]) (Earth-Moon system), (e.g., [67]) (Saturn-Titan system), and (e.g., [68]) (Pluto-Charon system). Therefore it would be very interesting if we could determine the similarities and the differences between the RTBP and the Hill problems. However there are two main reasons which make this task impossible: (i) the scattering region is not the same in all these works. In some cases the scattering region is extended around both primaries, while in some other cases it is confined very close 3 The two saddle points L 1 and L 2 lie along the x axis and the zero velocity surfaces are symmetric with respect to the three axes. It is easy to prove that if (x(t), y(t), z(t)) is a solution of the Hill problem, (−x(t), −y(t), −z(t)) is also a solution, of course with appropriately chosen momenta (because of the presence of the Coriolis force). Therefore, for every orbit escaping through L 1 there is a symmetric-related orbit escaping through L 2 , because the dynamics of the system must be symmetric. On this basis, the observed preference for escape through exit channel 2 is due to the particular choice of the initial conditions. Similar choices of initial conditions of orbits can be found in several earlier related papers (e.g., [13,33,34,[67][68][69]) to the secondary. (ii) the criteria for distinguishing between the several types of the orbits are different. For instance, in the present Hill problem there is only one type of escaping orbits, while on the other hand in the previous cases there were two types of escape (escape toward the primary realm and escape toward the exterior region). Therefore, we believe that if we try to compare all the different results we might end to erroneous and/or misleading conclusions due to the two above-mentioned reasons.
Taking into account the novel outcomes of our detailed numerical investigation we could argue that out task has been successfully completed. We hope that the present results to be useful in the field of the orbital dynamics of the Hill problem. In the next two papers of the series we will explore the orbital structure of the Hill problem where the perturbations of the oblateness (Part II) and the radiation pressure (Part III) are present.