On the classification of orbits in the three-dimensional Copenhagen problem with oblate primaries

TL;DR

This study classifies various orbit types in a three-dimensional restricted three-body problem with oblate primaries, revealing how initial conditions and oblateness influence orbital behavior and collision probabilities.

Contribution

It provides a detailed numerical classification of orbit types in the 3D oblate primaries problem, including the impact of oblateness on collision orbits.

Findings

Orbits are classified into four types based on initial conditions and energy.

Collision orbit fraction depends algebraically on oblateness coefficient.

Color-coded diagrams effectively visualize orbit distributions.

Abstract

The character of motion for the three-dimensional circular restricted three-body problem with oblate primaries is investigated. The orbits of the test particle are classified into four types: non-escaping regular orbits around the primaries, trapped chaotic (or sticky) orbits, escaping orbits that pass over the Lagrange saddle points $L_2$ and $L_3$, and orbits that lead the test particle to collide with one of the primary bodies. We numerically explore the motion of the test particle by presenting color-coded diagrams, where the initial conditions are mapped to the orbit type and studied as a function of the total orbital energy, the initial value of the $z$-coordinate and the oblateness coefficient. The fraction of the collision orbits, measured on the color-coded diagrams, show an algebraic dependence on the oblateness coefficient, which can be derived by simple semi-theoretical arguments.

Figures

Equations

$$\Omega (x,y,z)$$

$$+\frac{n^2}{2}\left(x^2+y^2\right)+\frac{\mu \left(1-\mu \right)}{2},$$

$$r_1$$

$$r_2$$

$$n=\sqrt{1+3\left(A_1+A_2\right)/2}.$$

$$\ddot{x}$$

$$\ddot{y}$$

$$\ddot{z}$$

$$J(x,y,z,\dot{x},\dot{y},\dot{z})=2\Omega (x,y,z)-\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)=C,$$

$$\dot{y_0}(x_0,y_0,z_0;C)=-\sqrt{2\Omega (x_0,y_0=0,z=z_0)-C}.$$

$$\mathrm{SALI}(\mathrm{t})\equiv \min ({\mathrm{d}}_{-},{\mathrm{d}}_{+}),$$

$$d_{-}$$

$$d_{+}$$

$$\dot{(\delta x)}$$

$$\dot{(\delta y)}$$

$$\dot{(\delta z)}$$

$$(\dot{\delta \dot{x}})$$

$$(\dot{\delta \dot{y}})$$

$$(\dot{\delta \dot{z}})$$

$$f_1(x,C)=2\Omega (x,y=0,z=z_0)=C,$$

$$f_2(y,C)=2\Omega (x=0,y,z=z_0)=C.$$

$$E=-\frac{1}{r_a+r_p}\pm {\left(\frac{2r_p{r}_a}{r_p+r_a}\right)}^{1/2},$$

$$r_a(r_p)=-r_p-\frac{E+r_p^2}{E^2-2r_p}-\frac{{\left(r_p^4+2E{r}_p^2+2r_p\right)}^{1/2}}{E^2-2r_p}.$$

$$S_{\text{collision}}\sim R^{\beta },$$

$$A=\frac{R_e^2-R_p^2}{5R_d^2},$$

Full text

††journal: International Journal of Non-Linear Mechanics

1 Introduction

A simple model in celestial mechanics is the planar restricted three-body problem (RTBP), that was introduced early on by Euler. In the model, two primary bodies rotate, around their center of mass in a plane and a test particle is restricted to this plane and moves only under the gravitational attraction of the primaries [64]. The chaotic dynamics of the test particle results from its non-integrability and it is a very active field of research, until today (e.g., [11, 6, 12, 19, 20, 21, 22, 23, 24, 25, 30, 31, 32, 33, 37, 47, 64, 65, 66, 60]). Both theoretical and numerical results are presented in the famous papers of Szebehely (e.g., [64, 65, 66, 32, 60]) and Hénon (e.g., [21, 22, 23, 24, 25]), where the two pioneers devoted major parts to find, describe, and classify periodic orbits. The RTBP became a prototypical model for deterministic chaos (e.g., [47]). Extensions considered relativistic dynamics (e.g., [11, 33]), quantum mechanics (e.g., [6]), and chemical (e.g, [12, 30]), or even astrophysical (e.g., [19, 20, 37]). The RTBP has been used for investigating the stability of several types of solar systems (e.g., [31]), chaos-assisted asteroid capture (e.g., [5]), as well as for modelling the orbital dynamics of a pair of black holes, orbited by a sun (e.g., [46]) – to mention only a few applications.

Usually, collision orbits can be considered as some kind of “leakage” of the phase space (e.g., [4, 7, 13, 14, 15, 18, 29, 40, 41, 48, 68]). In the classical version of the restricted three-body problem the shape of the two primary bodies is assumed, in most of the cases, to be spherically symmetric. However, in nature, and particularly in space, there are many celestial bodies, such as Saturn and Jupiter, which have in fact oblate shape (e.g., [9]) or even a triaxial one. Moreover, completely irregular shapes also exist. Such celestial bodies with irregular shapes are mainly asteroids and meteoroids (e.g., [36, 42]). The main effect on a celestial body, which is directly related to the oblateness and triaxiality, is a perturbation which leads to deviation from the classical two-body motion. The influence of the oblateness coefficient, when the more massive primary body is an oblate spheroid, has been studied in a series of works (e.g., [10, 26, 27, 28, 34, 35, 44, 49, 50, 51, 52, 55, 56, 57, 58, 62, 63]).

By today it is a well known fact that almost all the planets of our solar system are in fact oblate spheroidals. Such deviations from spherical symmetry have been taken into consideration in many recent studies. For instance for the Saturn-Dione-satellite and Saturn-Tethys-satellite systems it was proved in [43] that the convergence between the theoretical data and the corresponding ones from observations is much more better, when the oblateness of Saturn is included in the equations. In the same vein, in [61] the oblateness of Neptune was taken into consideration for modeling the motion of the spacecraft in the Neptune-Triton system.

In this work, we numerically investigate two aspects of the RTBP. First, the three dimensional version of this model is studied, where the test particle is not restricted to the plane but it can move freely in three dimensions. This extension was also investigated in [72]. Second, being the main novelty, we determine the impact of oblateness of the primary bodies on escape orbits, bounded orbits, and collision orbits. Our study complements previous works, where bounded, escape and collision motion was studied in a series of papers for the 2-dof RTBP (e.g. [69, 70, 71]).

The present paper is structured as follows: the basic properties of the dynamical system are presented in Section 2. All the computational methodology which will be used is described in detail in Section 3. The following Sections contains all the numerical results, regarding the orbit classification in the 3-dof RTBP. Our paper ends with Section 6, where the concluding remarks are emphasized.

2 Description of the dynamical system

The dynamical system consists of two primary bodies, $P_{1}$ and $P_{2}$, which perform circular Keplerian orbits around their common mass center [64] (see Fig. 1). The third body moves under the combined gravitational attraction of the two primaries. Considering that the mass of the third body $m$ is considerable smaller, with respect to the masses of the two primary bodies, we may assume the motion of the primaries is not perturbed, in any way, by the test particle.

In our system the gravitational constant $G$, the distance $R$ between the primaries and the sum of their masses are equal to unity. The dimensionless masses of the primary bodies are $m_{1}=1-\mu$ and $m_{2}=\mu$, where of course $\mu=m_{2}/(m_{1}+m_{2})\leq 1/2$ is the mass parameter. Furthermore, the centers of both primary bodies lie on the $x$-axis, at $(x_{1},0,0)$ and $(x_{2},0,0)$, where $x_{1}=-\mu$ and $x_{2}=1-\mu$. In the present study we shall consider the Copenhagen case, that is the case where the primaries have equal masses $m_{1}=m_{2}$ and therefore $\mu=1/2$. We consider a dimensionless, barycentric, rotating system of coordinates $Oxyz$, in which the $Ox$ axis always contains the two primary bodies, while the center of mass coincides with the origin $(0,0,0)$.

According to [53, 43, 16, 1] the time-independent effective potential function of the RTBP with oblate primaries is

[TABLE]

where

[TABLE]

are the distances of the third body from the respective primaries. Moreover, $A_{i}$, $i=1,2$ are the oblateness coefficients, while $n$ is the angular velocity which is defined as

[TABLE]

In our study we shall consider the case of relatively low values of oblateness coefficient, where $A_{1}=A_{2}=10^{-4}$. This choice can be justified if we take into consideration that the values of the oblateness in our Solar System are very low (see e.g., Table 1 in [54]).

The equations describing the motion of the test particle, in the corotating frame of reference, read

[TABLE]

The Jacobi integral of motion is described by the Hamiltonian

[TABLE]

where $\dot{x}$, $\dot{y}$, and $\dot{z}$ are the velocities, while $C$ is the conserved value of the Hamiltonian.

In the classical RTBP (that is when $A_{1}=A_{2}=0$) there are five equilibrium points, which are also known as Lagrange points. All five equilibrium points are located on the configuration $(x,y)$ plane with $z=0$. The central point $L_{1}$ is located between the two primaries, $L_{2}$ is located at the right side of primary $P_{2}$ (with $x>0$), while $L_{3}$ is located at the left side of primary $P_{1}$ (with $x<0$). In addition $L_{4}$ has $y>0$, while the libration point $L_{5}$ has $y<0$. In our case, where $A_{1}=A_{2}=10^{-4}$ the positions of the five equilibrium points are the following:

Equilibrium point $L_{1}$: $(x,y)=(0,0)$.

• 2.

Equilibrium point $L_{2}$: $(x,y)=(1.19839762,0)$.

• 3.

Equilibrium point $L_{3}$: $(x,y)=(-1.19839762,0)$.

• 4.

Equilibrium point $L_{4}$: $(x,y)=(0,0.86596768)$.

• 5.

Equilibrium point $L_{5}$: $(x,y)=(0,-0.86596768)$.

The values of the Jacobi constant at the Lagrange points $L_{i},i=1,...,5$ are denoted by $C_{i}$. For the Copenhagen problem with oblateness we have: $C_{1}=4.2508$, $C_{2}=C_{3}=3.70738405$, and $C_{4}=C_{5}$ = 3.00032499.

The relation $2\Omega(x,y,z)=C$ defines a three-dimensional surface, which is known as the zero velocity surface (ZVS) (see e.g., [17]). The ZVS allows us to determine the energetically possible orbits of the test particle, while its projection on the configuration $(x,y)$ plane is called Hill’s regions. In the same vein, the boundaries of the Hill’s regions are called zero velocity curves (ZVCs), since they are the locus where the kinetic energy is zero. The geometric shape of the ZVSs, and of course of the corresponding ZVCs strongly depends on the Jacobi constant. More precisely, four distinct energy regions exist regarding the configurations of the Hill’s regions:

Energy region I: $C>C_{1}$: All communication channels are closed, so the test particle moves either very close or far away from the primaries.

• 2.

Energy region II: $C_{1}>C>C_{2}$: Only the channel around $L_{1}$ is open and therefore orbits are allowed to communicate between the primary bodies.

• 3.

Energy region III: $C_{2}>C>C_{4}$: Both channels around $L_{2}$ and $L_{3}$ are open, which implies that the test particle can freely escape, through two directions.

• 4.

Energy region IV: $C<C_{4}$: The energetically forbidden regions are confined, while motion over the entire $(x,y)$ plane is possible.

The evolution of the geometry of the ZVSs, for several values of the Jacobi constant, is presented in Fig. 2.

3 Computational methodology

For revealing the orbital properties of the Hamiltonian system we have to determine the nature of several sets of initial conditions of orbits. In 3-dof systems the dynamics along the direction of the $z$ axis is very important. On this basis, we choose to classify initial conditions on the $(x,z)$ plane. More precisely, for several values of the Jacobi constant $C$, we create uniform grids of $1024\times 1024$ initial conditions $(x_{0},z_{0})$, inside a square area with $-4\leq x,z\leq 4$. For all orbits $y_{0}=\dot{x_{0}}=\dot{z_{0}}=0$, while the initial value of the velocity $\dot{y}$ is calculated as follows

[TABLE]

This particular choice of the initial conditions has a significant advantage, which will be presented and explained in the following Section.

The nature of the initial conditions of the orbits can be categorized as follows:

Bounded orbits. 2. 2.

Escaping orbits. 3. 3.

Orbits that collide with one of the primaries.

Furthermore, bounded orbits will be further classified into three sub-categories:

Non-escaping regular orbits. 2. 2.

Trapped sticky orbits. 3. 3.

Trapped chaotic orbits.

It should be noted, that there is also the type of orbits which pass over the inner equilibrium point $L_{1}$. However, we will not include this type of orbits in our study. This choice is justified because our classification reveals the final state of the orbits. On the other hand, orbits that pass over $L_{1}$ can later collide with one of the primaries, or escape from the system through $L_{2}$ or $L_{3}$, or even display trapped chaotic motion around the primaries. Therefore, passing through $L_{1}$ is only a first stage, while in this work we are interesting about the final state of the orbits. Moreover, we follow the same orbit classification which has been adopted in previous related papers (e.g., [38, 39, 69]).

For determining the regular or chaotic character of the orbits we use the Smaller Alignment Index (SALI) [59], which is defined as

[TABLE]

where

[TABLE]

are the alignments indices, while ${\vec{w_{1}}}(t)$ and ${\vec{w_{2}}}(t)$, are two deviation vectors which initially are orthonormal and point in two random directions. The variational equations needed for the computation of the two deviation vectors are given by

[TABLE]

For distinguishing between bounded and unbounded (escaping) motion we consider a limiting sphere with center at the origin $(0,0,0)$ and radius $R_{d}$. Then the case of bounded motion is present if the test particle stays confined, for integration time $t_{\rm max}$, inside the limiting sphere. On the other hand, escaping motion occurs if the test particle intersects the limiting sphere with velocity pointing outwards. For our calculations we choose $t_{\rm max}=10^{4}$ dimensionless time units and $R_{d}=10$. Moreover, a collision occurs if the test particle crosses the sphere with radius $R_{\rm col}$ around the primaries, where $R_{\rm col}=10^{-4}$. For both the parameters $R_{d}$ and $R_{\rm col}$ we adopted the values also used in [38, 39, 69, 71].

For the numerical integration of the equations of motion (4) and the corresponding variational equations (9) we used a double precision Bulirsch-Stoer FORTRAN 77 algorithm [45]. In all our calculations the numerical error, related to the conservation of the Jacobi integral (5) was generally smaller than $10^{-12}$, while in most of the cases it was smaller than $10^{-14}$. All the graphical illustration has been created using the latest version 11.2 of the software Mathematica*®* [67].

4 Orbital dynamics

In this Section we shall present color-coded diagrams (CCDs) thus following the methods introduced and used in [38] and [39]. In these diagrams, each pixel corresponds to a specific initial condition and it is colored according to the classification of the particular orbit. Throughout the section we will use several technical terms for describing the orbital properties of the system. The most important terminology is explained in A, where we also describe how one can interpret the different aspects of the color-coded diagrams.

4.1 Energy case I $(C\geq C_{1})$

Our first case concerns the scenario where all the transport channels are closed. In Fig. 3(a-b) the CCDs reveal the orbital structure of the configuration $(x,z)$ space for two values of the Jacobi constant $C$. The black solid lines denote the zero velocity curves, while the inaccessible energetically forbidden regions are marked in gray. When $C=5$, it is seen in panel (a) of Fig. 3, that inside the disks around the primaries the test particle can perform only two types of motion. More specifically, it can either move in bounded orbits around each primary, or collide with them. The SALI method indicates that the vast majority of the bounded orbits corresponds to non-escaping regular motion, while thin layers of trapped chaotic motion are also present, mainly in the vicinity of the basin111A local set of initial conditions of orbits which with the same final state is called a “basin”. boundaries of the stability islands. Outside the energetically forbidden regions vast basins of escaping orbits exist, while there are also stability islands. Additional numerical calculations suggest that these regular orbits circulate around both primary bodies. In panel (b) of Fig. 3 we observe that when $C=C_{1}$ there are several differences regarding the orbital structure of the $(x,z)$ plane, with respect to panel (a). In particular: (i) the amount of collision orbits inside the central regions has been increased; (ii) the presence of trapped chaotic orbits in the central regions is weaker; (iii) the stability islands in the exterior region (which correspond to regular motion around both primaries) have come closer to the energetically forbidden regions.

Because the two primary bodies have equal masses $(\mu=1/2)$, as well as the same degree of oblateness $(A_{1}=A_{2})$, the dynamical system has several symmetries. Some of them are in fact displayed in the configuration $(x,z)$ plane, where as we can see in Fig. 3, the structures of the CCDs are symmetrical with respect to both the horizontal and the vertical axes. These symmetries are revealed due to the specific choice of the initial conditions of the orbits, explained in the previous section. Indeed, using the specific choice of the initial conditions the several types of basins display the expected symmetries. For example, the collision basins to two primaries are exactly the same (meaning that the amount of orbits which lead to collision with primary $P_{1}$ is exactly the same with the amount of orbits which collide with primary $P_{2}$). The same property applies also for the bounded basins, corresponding to the two primaries. As far as we know, this is the only choice of initial conditions for which the CCDs, on the $(x,z)$ plane, display such symmetries, regarding the several types of basins.

In the following Fig. 4(a-b) we illustrate how the escape and collision time of orbits are distributed on the configuration $(x,z)$ space for the two values of the Jacobi constant discussed in Fig. 3(a-b). Light reddish colors correspond to fast escaping/collision orbits, dark blue/purple colors indicate large escape/collision time, while white color denote initial conditions of both regular non-escaping and trapped chaotic orbits. Note that the scale on the accompanying color bar is logarithmic. We see that some orbits collide with the primaries almost immediately. Looking more carefully at panel (b) of Fig. 4 we clearly observe that for a small portion of escaping orbits the corresponding escape times are at least two orders of magnitude larger than those of the vast majority of the escaping basins. These initial conditions of escaping orbits with high escape periods are mainly located in the vicinity of the boundaries between escaping and bounded basins.

4.2 Energy case II $(C_{2}\leq C<C_{1})$

The next case considers the second Hill’s region configurations in which the test particle is allowed to move between the two primary bodies, inside the interior region, through the open neck around $L_{1}$. The orbital structure of the configuration $(x,z)$ space is unveiled in Fig. 5(a-b) through the CCDs. In panel (a), where $C=4.22$, one can see that the orbital structure is very similar to that observed earlier in Fig. 3b for $C=C_{1}$. Around the saddle point $L_{1}$ there is a fractal222We would like to clarify that when we use the term “fractal” for describing an area it means that this particular local region displays a fractal-like geometry. For more information the fractal dimension [2, 3] should be computed. mixture of initial conditions which correspond either to collision or bounded motion. As the value of the Jacobi constant decreases, which means that the corresponding total orbital energy increases, there are two significant changes regarding the orbital content of the configuration $(x,z)$ plane. More precisely, in panel (b) of Fig. 5, which corresponds to $C=C_{2}$, it is seen that initial conditions of orbits that collide to primary 1 are present in the realm of primary 2 and vice versa. The other noticeable alteration concerns the geometry of the stability islands located in the exterior region.

The corresponding distribution of the escape and collision times of orbits on the configuration $(x,z)$ space is shown in Fig. 6(a-b). We observe that the escaping orbits, with initial conditions in the exterior region, can be divided into two categories: (i) those, just outside the energetically forbidden regions, with relatively low escape periods ($t_{\rm esc}<50$ time units), corresponding to dark blue colors and (ii) those with considerable higher values of escape rates ($t_{\rm esc}>100$ time units) which correspond to milder intermediate colors.

4.3 Energy case III $(C_{4}\leq C<C_{2})$

Energy case III concerns the third Hill’s region configurations where the test particle, with initial conditions in the interior region, can escape from the system through the escape channels in the ZVSs, located near the Lagrange points $L_{2}$ and $L_{3}$. The following Fig. 7(a-d) shows the orbital structure of the configuration $(x,z)$ space, for values of the Jacobi constant. When $C=3.7$ it is seen in panel (a) of Fig. 7 that the pattern of the configuration $(x,z)$ plane is almost identical to that of panel (b) of Fig. 5. Even though the necks around the saddle points $L_{2}$ and $L_{3}$ are open there is no numerical indication of communication between the interior and the exterior realms. In particular, initial conditions of collision and escaping orbits are present only in the interior and in the exterior region, respectively. It seems that for such energy levels the two escape channels, around $L_{2}$ and $L_{3}$, are not wide enough so as to allow orbits to escape, even if this option is theoretically allowed. As the value of the Jacobi constant decreases the structure of the $(x,z)$ plane changes drastically (see panels (b-d) of Fig. 7).

The most important changes that occur are the following: (i) the channels around the saddle points $L_{2}$ and $L_{3}$ become more and more wide; (ii) the stability islands in the exterior region split off into four pieces; (iii) initial conditions of escaping orbits appear in the interior region, while at the same time initial conditions of collision motion populate the exterior realm; (iv) the initial conditions in the exterior region, corresponding to collision motion, seem to form waves which pass through the openings of the exterior stability islands. The formation of the wave-shaped structures is due to the specific choice of the initial conditions. Evidently, using other type of initial conditions (see e.g., [8]) the shape of the collision basins will be different.

In Fig. 8(a-d) we illustrate the corresponding distribution of the escape and collision time of orbits on the configuration $(x,z)$ space. This type of plot allows us to monitor in more detail the evolution of the structure of the central stability islands. We observe that when $C=C_{4}$ there exist only two stability islands in the central region, corresponding to non-escaping orbits that circulate around each primary body.

4.4 Energy case IV $(C<C_{4})$

The last energy case involves the scenario when $C<C_{4}$. Fig. 9(a-d) reveals the orbital structure of the configuration space through the CCDs. In panel (a) of Fig. 9 it is seen that when $C=2$ the vast majority of the configuration $(x,z)$ plane is covered by initial conditions of escaping orbits. Inside the vast escaping basin we can distinguish a plethora of smaller basins of initial conditions corresponding to collision motion. Furthermore, bounded motion is still possible. Indeed, stability islands of non-escaping regular orbits circulating around each and both primaries are present. According to panel (b) of Fig. 9, when $C=0.5$ there are three main differences with respect to the previous case: (i) the energetically forbidden regions have disappeared, (ii) the amount of collision orbits has been reduced from about 24% to 6% and (iii) non-escaping regular motion around each primary body is no longer possible. With decreasing value of the Jacobi constant the portion of initial conditions which correspond to collision motion continues to reduce. Finally when $C=-2$ (see panel (d) of Fig. 9) collision motion is almost negligible, while the $(x,z)$ plane is dominated only by escaping and non-escaping regular orbits around both primaries. The corresponding distribution of the escape as well as the collision times of orbits on the configuration space is depicted in Fig. 10(a-d). One can see similar outcomes with that presented earlier in the previous subsections.

4.5 An overview analysis

The CCDs in the $(x,z)$ space provide sufficient information for only a fixed value of the Jacobi constant. It would be very illuminating if we could scan a continuous spectrum of Jacobi constants $C$ rather than a few discrete values. In order to do this, we use new types of two-dimensional planes in which the value of the Jacobi constant is the ordinate, while the abscissa will be the $x$ or the $y$ coordinate of the orbits [23]. In particular, in the $(x,C)$ plane all orbits are launched from the $x$-axis, with $z=z_{0}$ and $y_{0}=\dot{x_{0}}=\dot{z_{0}}=0$, while the initial value of $\dot{y}$ is always derived through the Hamiltonian (5), in the form of Eq. (6). Similarly, in the $(y,C)$ plane the initial conditions are started form the $y$-axis, with $z=z_{0}$ and $x_{0}=\dot{y_{0}}=\dot{z_{0}}=0$, while the initial value of $\dot{x}$ is obtained as $\dot{x_{0}}(x_{0}=0,y_{0},z_{0};C)=\sqrt{2\Omega(x_{0}=0,y_{0},z=z_{0})-C}$.

In Fig. 11(a-f) we present the orbital structure of the $(x,C)$ and $(y,C)$ planes, for three values of $z_{0}$, when $C\in[-6,6]$. The black solid lines are the limiting curves which distinguish between regions of energetically allowed and forbidden motion and they are defined as follows

[TABLE]

In panels (a) and (b), where $z_{0}=0.01$ we observe that between the energetically forbidden regions, regular bounded basins corresponding to quasi-periodic orbits around one of the primaries are present. In the exterior region we have a strong presence of stability islands which correspond to ordered orbits that circulate around both primaries. Undoubtedly, the most interesting orbital structure exists between these stability islands, whose shape resembles a claw. Indeed, in these regions several basins of collision coexist with a large portion of basins of escape. Outside the bounded basins of the exterior region a escaping motion always dominates. The distribution of the corresponding escape and collision time of the orbits is depicted in Fig. 12(a-f).

As we proceed to higher values of the initial coordinate $z_{0}$ the following changes, regarding the orbital structure of the $(x,C)$ and $(y,C)$ planes, take place: (i) the three energetically forbidden regions merge together, thus forming a unified forbidden region the area of which increases with increasing $z_{0}$ (ii) non-escaping regular motion around each primary is no longer possible, due to the unification of the energetically forbidden regions, (iii) the number of collision orbits, with initial conditions between the stability islands of the exterior region, is reduced, and (iv) the edges of the stability islands in the exterior region come closer and tend to merge near the boundaries of the energetically forbidden region.

In panels (a-f) of Fig. 13 we present the evolution of the rates of all types of orbits, as a function of $C$. Looking the diagrams of this figure we can draw the following conclusions:

Escaping motion dominate (with percentages larger than 90%) at very low $(C>6)$, very high $(C<-5.5)$ and also at intermediate energy levels $(-1.5<C<3.5)$. 2. 2.

Collision motion exist mainly when $-1.5<C<3.5$, while for all other values of the Jacobi constant the corresponding percentages tend to zero. 3. 3.

Bounded regular motion exhibits two peaks at the exact energy levels where the lowest rates of the escaping orbits are observed.

5 Influence of the oblateness coefficient

So far we have established the influence of the Jacobi constant (energy level) as well as of the initial value of the $z$-coordinate on the nature of the orbits. In [69] and [71] we demonstrated how the oblateness coefficient affects the character of the orbits in the planar version of the Copenhagen problem. Therefore it would be of great interest to determine the influence of the same dynamical quantity on the orbital structure of the spatial version of the Copenhagen problem.

In order to obtain this information we classified sets of initial conditions or orbits on the $(x,C)$ plane, for six values of the oblateness coefficient, for fixed $z_{0}=0.01$. Our results are depicted in Fig. 14(a-f), while the corresponding distributions of the escape and collision time of the orbits are given in Fig. 15(a-f). It is visible that as the primary bodies become more oblate several important changes occur regarding the orbital structure of the system. In particular:

The portion of the initial conditions of orbits, which lead to collision to one of the primaries, significantly increases. 2. 2.

The stability basins, corresponding to non-escaping regular motion around each primary, are reduced and for extremely high values of the oblateness ($A>0.01$) they completely disappear. On the other hand, the stability islands, corresponding to regular orbits around both primaries, seem to be completely unaffected by the shift of the oblateness coefficient. 3. 3.

The fractality of the phase space is considerably been reduced. Indeed, all the regions, containing fractal mixtures of different types of orbits, vanish due to the emergence of extended well-defined basins of collision.

Strikingly enough similar behavior, regarding the influence of the oblateness coefficient on the character of the orbits, has also been reported in [69] and [71]. Therefore we may conclude that the oblateness of the primaries affects, in a very similar way, the nature of the orbits in both 2-dof and 3-dof versions of the circular restricted three-body problem.

Our numerical calculations indicate that the oblateness affects also the escape as well as the collision time of the orbits. Specifically, we observed that both the escape and the collision time of the orbits are reduced (at least for certain energy ranges) with increasing value of the oblateness coefficient.

As can be seen in Fig. 13, the blue (collision with primary $P_{1}$) and the red (collision with primary $P_{2}$) curves lie directly on top of each other such that the blue is not visible. This is of course due to the symmetries of the problem ($\mu=1/2$, and $A_{1}=A_{2}$). As can be seen in Fig. 14, the phase space diagrams (i.e. the orbital structure) is indeed anti-symmetric regarding red and blue collision basins, as it should be for $A_{1}=A_{2}$ and $m_{1}=m_{2}=1/2$, given that the $(x,C)$ plane is displayed but the primaries are at positions $(x=-\mu,y=0)$ and $(x=1-\mu,y=0)$ in the corotating frame. Or in short, the primaries are located at $x=-1/2$ and $x=1/2$, hence the observed anti-symmetry, regarding the $x$-position (for $A_{1}=A_{2}$).

Very similar outcomes can be extracted also from the $(y,C)$ planes. However, mainly for saving space, we do not present the corresponding diagrams. In contrast, it is illuminating to discuss the relation between the area of the different types of basins and the value of the oblateness coefficient $A$. Thus, in Fig. 16 we depict the evolution of the area $S$ of the different types of basins on the $(x,C)$ plane (see Fig. 14), as a function of $A$. In this diagram we do not include the parametric evolution of the trapped chaotic orbits because the corresponding percentages are extremely low. We observe that the area of the non-escaping regular orbits remains almost unperturbed and only for high values of $A$ displays a small reduction. Moreover, the area of the basins of escape is constantly reduced with the increase in the value of the oblateness coefficient. Finally, the most interesting behavior is that of the collision orbits, where it is seen that the area of the corresponding basins exhibits a linear increase in a log-log scale, i.e., a power law. At this point, it should be noted that in [38] (see Fig. 6) and [39] (see Fig. 11) a similar behavior was observed, regarding the area of the basins of collision and the corresponding radius of collision. When $\log_{10}(A)=-1$ collision orbits (to both primaries) occupy half of the $(x,C)$ plane, while non-escaping regular and escaping orbits share the rest half of the same plane.

In Fig. 16 we find a power law $S_{\text{collision}}=A^{\gamma}$, with $\gamma=0.181$. In what follows we try to explain the power law behavior.

For the classical restricted three-body problem an approximation for the power law behavior of the area of the collision orbits in two-dimensional (2D) Poincaré surfaces of sections, $S_{\text{collision}}\sim R^{\beta}$, was derived, predicting the exponent $\beta=1/2$ (see [38]). The RTBP can be roughly approximated by the Kepler problem when the test particle is close to one primary body, just before a collision occurs, or when the rotation of the primaries is neglected. Following the approach used in [38] the total orbital energy, in the inertial system, is given by $E_{\text{in}}=E+L$. In our notation we have $C=-2E$.

From Kepler’s ellipse expression we obtain

[TABLE]

where $r_{p}$ denotes the perihelion and $r_{a}$ the aphelion distance. Solving Eq. (11) for $r_{a}$ yields

[TABLE]

A collision occurs when the test particle intersects the circle of radius $R$ around the Kepler singularity, $r_{p}\leq R$. Hence, for $R\ll 1$, $S_{\text{collision}}\approx 2\pi r_{a}(0)\left(r_{a}(0)-r_{a}(R)\right)$. Assuming $r_{p}\ll 1$, Eq. (12) shows approximately a square root behavior as a function of $r_{p}$, $r_{a}(r_{p})\approx-\frac{1}{E}+\frac{(2r_{p})^{1/2}}{E^{2}}$. For the collision orbit area $S_{\text{collision}}$, this translates to

[TABLE]

with the exponent $\beta=1/2$.

To link this result to the collision area as a function of the oblateness coefficient $A$, as presented in Fig. 16, it is useful to recall the definition of the oblateness coefficient,

[TABLE]

where $R_{e}$ is the equatorial radius and $R_{p}$ denotes the polar radius, while $R_{d}$ is the distance between the centers of the two primaries [54]. For a large equatorial oblateness, we find approximately that $A\sim R_{e}^{2}:=R^{2}$. Thus, using Eq. (13) this provides a rough approximation for the collision area as a function of the oblateness coefficient, giving the power law $S_{\text{collision}}\sim R^{\beta}\sim A^{\beta/2}=A^{\gamma}$, with $\gamma=1/4$.

The power law fit of Fig. 16 gives $\gamma=0.181$, which is a lower exponent that obtained from the rough approximation. This, however, may be due to the fact that the approximation does not account for the nonlinear dependencies between $A$ and the radii in Eq. (14), which may crucially determine the dynamics, collision orbits included. For the same reason, as it is clear from Fig. 16, the areas of the three classes in Fig. 16 depend on the oblateness coefficient and do not add up to a constant ($S_{\text{collision}}+S_{\text{bounded}}+S_{\text{escape}}\neq const.$).

6 Conclusions

In this paper we have studied the character of motion of a test particle in the three dimensional Copenhagen problem, where the two primary bodies are oblate spheroids. Our thorough and systematic numerical investigation allowed us to distinguish between bounded (with further classification between non-escaping regular and trapped chaotic orbits), unbounded and collision motion. We have studied the basins of escape and collision and related them to the corresponding escape and collision time of the orbits. Our analysis indicates that the three dimensional (3D) motion of the test particle in the system of two oblate spheroid primaries strongly depends on three main parameters: (i) the total orbital energy or equivalently the value of the Jacobi constant and (ii) on the initial value of the vertical $z$-coordinate of the orbits and (iii) the value of the oblateness coefficient.

To our knowledge this is the first time that the orbital dynamics, in the three dimensional Copenhagen problem with oblateness, is numerically investigated in a systematic manner, which represents our main contribution to existing literature.

The following list contains the most important conclusions of our numerical analysis.

Escaping motion completely dominates at both very low as well as very high levels of the total orbital energy. Furthermore, escaping motion is the only type of motion which is possible for all values of the energy. 2. 2.

For intermediate values of the energy we detected a very complicated mixture of all types of possible motion. In this regime all the types of basins exist and they display highly fractal basin boundaries. 3. 3.

The vast majority of the initial conditions corresponding to non-escaping bounded motion corresponds to regular orbits. In fact there are two main types of regular motion: (i) motion only around one of the primary bodies and (ii) motion around both primaries. Our numerical computations showed that regular motion mainly corresponds to simple loop orbits (not shown here), while stability islands of secondary resonances333In this study, with the term secondary resonances we refer to all regular orbits with a shape more complicated than a simple loop. are limited. 4. 4.

Collision motion was found to be the only type of motion which is strongly related with the initial value of the $z$-coordinate of the orbits. In particular, the higher the $z_{0}$ of the orbits the lower the corresponding amount of detected collision orbits. 5. 5.

Our calculations suggest that the fractality of the phase space depends on the initial value of the $z$-coordinate of the orbits. More precisely, we found fair numerical evidence that the phase space becomes less fractal as we proceed to higher levels of the $z$-coordinate. 6. 6.

The oblateness coefficient $A$ was found to influence mostly the collision orbits. In particular, the amount of the collision orbits increases rapidly, with increasing value of the oblateness coefficient. 7. 7.

The average escape and collision time of the orbits seems to be also connected with the oblateness coefficient. In fact, for specific energy regions, the escape and collision times were found to reduce with increasing value of $A$.

Using an Intel*®* Quad-CoreTM i7 2.4 GHz PC the required CPU time, for the classification of each grid of $1024\times 1024$ initial conditions, was varying between 14 hours and 5 days, depending of course on the percentage of escaping and collision orbits. We would like to point out that the total computational cost may be reduced to 1/2 (or even to 1/4) if we take into account all the symmetries of the system. Nevertheless, in our computations we conducted the full run just for being sure that everything works fine, regarding the choice of the initial conditions and the symmetries of the system.

All the results as well as the conclusions of our numerical exploration are considered as a step towards the continuous task of revealing and understanding the great complexity of the orbital dynamics of the spatial circular RTBP with perturbations (e.g., oblateness). It is in out future plans to continue the numerical exploration of the system with three degrees of freedom (3-dof) by considering other sets of initial conditions of orbits (e.g., $\left(x_{0},0,0,0,\dot{y_{0}},\dot{z_{0}}\right)$), so as to cover and unveil additional solutions of the system.

Acknowledgments

Our warmest thanks go to two anonymous referees for the careful reading of the manuscript as well as for all the apt suggestions and comments which allowed us to improve both the quality and the clarity of the paper.

Appendix A Terminology

In this Appendix we will provide the explanation as well as the definitions of some of the specialized terminology used in this paper. Fig. 17 will be used as a test figure in order to correlate the terms with the corresponding parts of the color-coded diagrams.

Energetically allowed regions: For every value of the Jacobi constant $C_{0}$ the areas of the phase space where $2\Omega(x,y,z)>C_{0}$ are called energetically allowed regions. In these regions the test particle is permitted to move, according to the equations of motion. On the other hand, the regions of the phase space where $2\Omega(x,y,z)<C_{0}$ are called energetically forbidden regions and the test particle is not allowed to enter them. The energetically forbidden regions are always shown in gray in the color-coded diagrams.

• 2.

Interior realm: The subregion of the energetically allowed areas between the two saddle points (that is for $x(L_{3})<x<x(L_{2})$) defines the interior realm or the central region of the phase space. In Fig. 17 the interior realm is defined as the area inside the black, dashed circle. All the energetically allowed areas outside this circle define the so-called exterior region or exterior realms.

• 3.

Transport channels: In Fig. 17 it is seen that there are two disjoint parts of energetically forbidden regions. These parts enclose the interior region of the phase space. Moreover, one can see, that near the equilibrium points $L_{2}$ and $L_{3}$ the energetically forbidden regions form two channels through which orbits from the interior region can enter the exterior region or vice versa. These channels, located in the vicinity of the saddle points, are the transport channels of the dynamical system.

• 4.

Stability islands: The green regions on the color-coded diagrams are formed by a collection of initial conditions corresponding to non-escaping regular orbits. These regions are called stability islands or bounded basins. In this dynamical system there are two main types of stability islands or bounded basins: (i) those which are located in the interior region and correspond to regular orbits that circulate around only one primary and (ii) those which are located in the exterior region and are composed of initial conditions of orbits that circulate around both primary bodies.

• 5.

Escaping and collision basins: We observe in Fig. 17 that a large portion of the $(x,z)$ plane is covered by initial conditions of escaping orbits which form the so-called basins of escape (cyan regions). In the same vein, the initial conditions of orbits which lead to collision with one of the primaries form the so-called basins of collision (blue and red regions). In some areas of the phase space the escape and collision basins are clearly distinguishable, while there are also regions where the initial conditions create a highly complicated fractal mixture of escaping and collision motion.