Orbital dynamics in the post-Newtonian planar circular Sun-Jupiter system
MMarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036
International Journal of Modern Physics Dc (cid:13)
World Scientific Publishing Company
ORBITAL DYNAMICS IN THE POST-NEWTONIAN PLANARCIRCULAR RESTRICTED SUN-JUPITER SYSTEM
EUAGGELOS E. ZOTOS
Department of Physics, School of Science,Aristotle University of Thessaloniki,GR-541 24, Thessaloniki,[email protected]
F. L. DUBEIBE
Facultad de Ciencias Humanas y de la Educaci´on,Universidad de los Llanos, Villavicencio, ColombiaGrupo de Investigaci´on en Relatividad y Gravitaci´on,Escuela de F´ısica, Universidad Industrial de Santander,A.A. 678, Bucaramanga 680002,Colombiafl[email protected]
Received 29 September 2017Revised 3 November 2017Accepted 14 November 2017Published 5 December 2017The theory of the post-Newtonian (PN) planar circular restricted three-body problem isused for numerically investigating the orbital dynamics of a test particle (e.g., a comet,asteroid, meteor or spacecraft) in the planar Sun-Jupiter system with a scattering regionaround Jupiter. For determining the orbital properties of the test particle, we classifylarge sets of initial conditions of orbits for several values of the Jacobi constant in allpossible Hill region configurations. The initial conditions are classified into three maincategories: (i) bounded, (ii) escaping and (iii) collisional. Using the smaller alignmentindex chaos indicator (SALI), we further classify bounded orbits into regular, sticky orchaotic. In order to get a spherical view of the dynamics of the system, the grids of theinitial conditions of the orbits are defined on different types of two-dimensional planes.We locate the different types of basins and we also relate them with the correspondingspatial distributions of the escape and collision time. Our thorough analysis exposes thehigh complexity of the orbital dynamics and exhibits an appreciable difference betweenthe final states of the orbits in the classical and PN approaches. Furthermore, our nu-merical results reveal a strong dependence of the properties of the considered basins withthe Jacobi constant, along with a remarkable presence of fractal basin boundaries. Ouroutcomes are compared with earlier ones, regarding other planetary systems.
Keywords : methods: numerical; celestial mechanics; chaos.PACS Number(s): 04.25.Nx, 05.45.-a, 45.50.Pk, 45.50.Jf1 a r X i v : . [ a s t r o - ph . E P ] M a r arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe
1. Introduction
In the strong field regime, the space-time describing binary systems is strongly non-linear and time-dependent due to the emission of gravitational waves. The complex-ity of this scenario is one of the reasons to explain why no exact solution to Einstein’sequations modeling a binary system in orbital motion has ever been found. Severalhighly sophisticated approaches to tackle this problem have been developed, includ-ing numerical solutions and approximation methods. In some systems with smallorbital velocities (compared to the speed of light) and weak gravitational fields,the flagship tool is the post-Newtonian (PN) approximation, developed since theyouth of general relativity. A natural scenario to observe non-negligible PN effectsover long time-scales is the Solar System. Indeed, it is a well-known fact that thePN corrections must be accounted for in cases such as in calculating the perihelionprecession of Mercury’s orbit.On the other hand, since Poincar´e’s demonstration about the impossibility tofind an exact solution for the motion of two planets around the Sun, the so-called three-body problem, and the existence of heteroclinic intersections in celestialbodies dynamics, it is generally accepted that the Solar System may be chaotic and therefore, there exist practical limits for the long-time predictions. Takinginto account that the Newtonian three-body problem is just a first approximationof a much more complex setting, it is desirable to refine the current understandingof the Solar System dynamics by including general relativistic corrections. To doso, some authors have derived the first-order PN equations of motion (1-PN) forthe circular restricted three-body problem, by using the Einstein-Infeld-Hoffmann (EIH) theory. Some studies have been carried out on the existence and stability of the equi-librium points of the relativistic problem of three bodies. The triangular librationpoints were studied in Ref. 4, finding that in contrast to the classical restrictedthree-body problem, in the relativistic system the equilibrium points L and L areunstable. Ref. 11 controverted this result, due to the existence of a region of lin-ear stability in the parameter space ( µ, /c ), with µ the mass parameter and c thespeed of light. A detailed analysis of the collinear points for several Sun-planet pairswas performed in Ref. 33, finding that for all cases examined the collinear points L , L , and L are unstable. Concerning the dynamics of the PN restricted three-bodyproblem, recent studies have shown that the PN terms can be understood as non-negligible perturbations to the classical system, and hence, for small distancesbetween the primaries the PN dynamics differ qualitatively from the Newtonianone. Undoubtedly, one of the most important aspects of the restricted three-bodyproblem (RTBP) is the classification of the initial conditions of the orbits. Know-ing the nature of the orbits has numerous modern applications, such as in spaceflight missions, in launching and positioning artificial satellites, while it also servesthe basis of several planetary and exoplanetary theories.
It all started witharch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036
Orbital dynamics of the Sun-Jupiter system the pioneer works of Refs. 27, 28, where the first thorough and systematic orbitclassification, in the classical RTBP, took place. During the last few years, theorbital dynamics of many planetary systems has been numerically investigated,through the classification of initial conditions in some particular scattering regions;the most notable works are the following: regarding the Earth-Moon system, forthe Saturn-Titan system and for the Pluto-Charon system. In this paper, we shall explore the orbital dynamics of the Sun-Jupiter system,when the PN correction terms are included in the equations of motion as well asin the variational equations. The scattering region will be located in the vicinityaround Jupiter and particularly between the Lagrange points L and L . Our maingoal is to classify initial conditions of the orbits of the test particle and determinethe corresponding types of basins (bounded, escape and collision).The paper is organized as follows: In section 2 we first identify the most im-portant properties of the dynamical system. Next, we outline the computationalmethods used for the classification of the initial conditions of the orbits in section3. In the following section, a thorough and systematic numerical investigation takesplace, thus revealing the orbital dynamics of the planar PN Sun-Jupiter system.Finally, our paper ends with section 5 where the discussion is given.
2. Properties of the dynamical system2.1.
Description of the mathematical model
With some necessary simplifying assumptions, the orbital motion of Jupiter aroundthe Sun can be considered circular (even though the eccentricity of its orbit isabout 0.04839266). Under this condition, the motion of a test particle (e.g. a comet,asteroid or spacecraft) can be modeled as a planar circular restricted three-bodyproblem (henceforth PCRTBP), in which the two bodies, the primary (Sun) m andthe secondary (Jupiter) m , move in circular orbits, with the same angular velocity,around their common centre of mass, at a fixed distance a . The test particle m ,whose mass is negligible compared to m and m , does not perturb the motion ofthe primary and secondary bodies, and moves on the same plane ( x, y ) under thecombined gravitational influence of the two main bodies.In order to facilitate the analysis of the system we use the Szebehely convention for the normalization of the constants m = 1 − µ, m = µ, x = − µ, x = 1 − µ = 1 + x , where µ = m / ( m + m ) ∈ [0 , / m > m . Withthis choice of units, a = 1 and the origin O is located at the centre of mass ofthe two main bodies. Moreover, the centers P and P of the two main bodies arelocated at ( x ,
0) and ( x , In aarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe synodic frame of reference, the equations of motion for the test particle, in accor-dance with Ref. 25, are given by the following set of differential equations¨ x = x + 2 ˙ y − m ( x − x ) d − m ( x − x ) d + (cid:15)c R x , (1)¨ y = y − x − y (cid:18) m d + m d (cid:19) + (cid:15)c R y , (2)where x i denotes, in the synodic frame of reference, the fixed position of the body m i ,while d i = (cid:113) ( x − x i ) + y . Furthermore, R x and R y are the relativistic correctionterms.The system of equations (1-2) admits one integral of motion, which is known asthe Jacobi integral (or Jacobi constant) and can be written as J ( x, y, ˙ x, ˙ y ) = m d + m d + 12 (cid:0) x − ˙ x + y − ˙ y (cid:1) + (cid:15)c J R = C, (3)where C is the Jacobi constant which is conserved.The parameter (cid:15) is a transition parameter where (cid:15) = 0 indicates classical New-tonian dynamics, while (cid:15) = 1 corresponds to PN dynamics. In what follows, we shallonly consider the PN case ( (cid:15) = 1), unless the contrary is explicitly stated.According to Ref. 12 the relativistic correction terms read R x = m m a (cid:18) x − x d + x − x d (cid:19) + y (cid:18) m d + m d (cid:19) × ( ˙ x − y ) ( x + 4 ˙ y ) + (cid:18) m ( x − x ) d + m ( x − x ) d (cid:19) × (cid:34) (cid:18) m d + m d (cid:19) − x ( y − ˙ x ) − ( x + ˙ y ) (cid:35) − (cid:18) m ( x − x ) x d + m ( x − x ) x d (cid:19) + 32 y (cid:18) m ( x − x ) x d + m ( x − x ) x d (cid:19) + 2 ω ( x + ˙ y ) + (cid:18) x y (cid:19) × (cid:18) m ( x − x ) x d + m ( x − x ) x d (cid:19) − (cid:18) m x d + m x d (cid:19) , (4)arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system R y = m m a y (cid:18) d + 1 d (cid:19) + 32 y (cid:18) m x d + m x d (cid:19) + 2 ω ( y − ˙ x ) + (cid:18) m d + m d (cid:19) × (cid:40) y (cid:20) (cid:18) m d + m d (cid:19) + 3 ˙ y ( x + ˙ y ) − y (cid:21) + ˙ x (cid:0) x + 3 x ˙ y + 2 y (cid:1) − ˙ x y (cid:41) − (cid:104) x (cid:16) ˙ x − y (cid:17) + 3 ˙ x ˙ y (cid:105) × (cid:18) m x d + m x d (cid:19) + (cid:18) x − y (cid:19) (cid:18) m x d + m x d (cid:19) − ( y − ˙ x ) ( x + ˙ y ) (cid:18) m ( x − x ) d + m ( x − x ) d (cid:19) , (5)and J R = 18 (cid:104)(cid:0) x + y (cid:1) − (cid:0) ˙ x + ˙ y (cid:1) (cid:105) − (cid:0) x + y (cid:1) (cid:0) ˙ x + ˙ y (cid:1) − ( x ˙ y − y ˙ x ) (cid:0) ˙ x + ˙ y (cid:1) + ω (cid:0) x + y (cid:1) −
12 ( x ˙ y − y ˙ x ) + 32 (cid:18) m d + m d (cid:19) (cid:0) x − ˙ x + y − ˙ y (cid:1) + 32 (cid:18) m x d + m x d (cid:19) − m m a (cid:18) d + 1 d (cid:19) − y (cid:18) m x d + m x d (cid:19) − (cid:18) m d + m d (cid:19) − x (cid:18) m x d + m x d (cid:19) , (6)The angular velocity ω was calculated by and is given by ω = m m − m + m ) a ( m + m ) , (7)which reduces to ω = [ µ (1 − µ ) − /
2, with the help of the mass parameter.It is important to note that in contrast to the classical Newtonian PCRTBP,the PN system depends on two parameters, the mass parameter µ and the speed oflight c . According to Ref. 29 for the Sun-Jupiter system the mass parameter is µ =0.000953817733371, and the speed of light in canonical units is c = 22945 . Equilibrium points
The equilibrium points (or Lagrange points) correspond to equilibrium points ofthe PCRTBP, and can be found by replacing ¨ x = ¨ y = ˙ x = ˙ y = 0 in the respectivearch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe equations of motion (1-2). In analogy with the classical Newtonian case, the PNPCRTBP has five equilibrium points (also known as Lagrange points): three ofthem, L , L , and L , are collinear points located on the x -axis, while the othertwo L and L are called triangular points and they are located on the vertices of anequilateral triangle. The central stationary point L is located between the primaryand the secondary, L is at the right side of the secondary (Jupiter), while L is atthe left side of the primary (Sun).We have solved numerically the resulting algebraic system ¨ x | ˙ x = ˙ y =0 = 0 and¨ y | ˙ x = ˙ y =0 = 0 finding approximate expressions for the Lagrange points in terms ofthe mass parameter µ . For the collinear points, we find that the classical expansionin powers of n/ µ in powers of n ,gives a very accurate solution with a maximum error of 10 − . Hence, the solutioncan be written as x ∗ L i = (cid:88) n =0 C in µ n + (cid:88) n =1 ˜ C in (cid:18) µ − µ (cid:19) n/ , (8)where i = 1 , ,
3, denotes the i -th collinear point. In a similar fashion, for thetriangular points we get x ∗ L j = (cid:88) n =0 T in µ n , y ∗ L j = (cid:88) n =0 ˜ T in µ n , (9)where j = 4 ,
5, denotes the j -th triangular point. The numerical coefficients enteringequations (8) and (9) are given in Appendix B.Additionally, in Table 1 we present the numerical values of the Lagrange pointsfor the Sun-Jupiter system in the framework of the PN-approximations (column 3)compared to the classical Newtonian values (column 2). It deserves mentioning thatthe numerical results presented in Table 1 satisfy the respective system of algebraicequation, with an accuracy of 10 − m, which indicates that the differences presentedin column 4 are not numerical artifacts.Our results suggest that the largest deviation is observed for the x -coordinateof L and L , while the smallest deviation takes place for the x -coordinate of L . Ingeneral, the displacements of the Lagrange points are such that the new positionsare closer to the primaries, i.e. L , L , L and L are moved toward Jupiter, while L is shifted toward the Sun.Coordinate Newtonian Gravity ( m ) PN-Approximation ( m ) Difference ( m ) x ( L ) 7.257656518990008 × × x ( L ) 8.319894593317031 × × -38.025 x ( L ) -7.787213864029970 × -7.787213864019399 × x ( L ) = x ( L ) 3.884635501925640 × × y ( L ) = − y ( L ) 6.741245897986063 × × -532.902arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system L and L , while the position of the centre of Jupiter ( P ) is indicated by a bluedot. (a): C = 1 . C = 1 . C = 1 . The values of the Jacobi integral at the Lagrange points L k , k = 1 , ..., C k and they are critical levels with values: C = 1.519379668835193, C =1.518743663753772, C = 1.500476898588919, and C = C = 1.499523545304652. Hill region configurations
The projection of the four-dimensional phase space onto the configuration ( x, y )plane is called the Hill region configurations, which constitute the energetically ac-cessible regions to the orbits for a given value C . The boundaries of these regionsare called Zero Velocity Curves (ZVCs) because they are the locus in the config-uration space where the kinetic energy vanishes. The value of the Jacobi constantstrongly dictates the structure of the corresponding Hill region configurations. Moreprecisely, there are five distinct cases: • Case I:
C > C : All channels are closed, so there are only bounded andcollisional motion (see panel (a) of Fig. 1). • Case II: C < C < C : Only the channel around L is open thus allowingthe test particle to enter the Sun realm (see panel (b) of Fig. 1). • Case III: C < C < C : The channel around L opens, so the test particlecan enter the exterior region and escape from the system (see panel (c) ofFig. 1). • Case IV: C < C < C : Both channels, around the Lagrange points L and L are open, therefore the test particle is free to escape through twodifferent directions. • Case V:
C < C : The energetically forbidden regions disappear, so motionover the entire configuration ( x, y ) space is possible.In Fig. 1(a-c) we present the structure of the first three Hill region configurations.arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe
3. Computational methods and criteria
For revealing the orbital dynamics in the PN Sun-Jupiter system, we need to numer-ically integrate the equations of motion (1-2), for several sets of initial conditions.For this purpose we consider dense uniform grids of 1024 × x , y ) regularly distributed on the configuration ( x, y ) plane inside the energet-ically allowed area defined by the corresponding Jacobi constant C . Following atypical approach, all orbits are launched with initial conditions inside a certainscattering region which in our case is x ( L ) ≤ x ≤ x ( L ) and − . ≤ y ≤ . x = 0, while the initial value of ˙ y is always derived by the Jacobiintegral (3) as ˙ y = ˙ y ( x , y , ˙ x , C ) > • Bounded orbits around the secondary (Jupiter). • Orbits that escape from the scattering region. • Orbits that collide with the secondary (Jupiter).Moreover, all bounded orbits will be further classified into three sub-categories: • Non-escaping regular orbits. • Trapped sticky orbits. • Trapped chaotic orbits.For distinguishing between regular and chaotic dynamics we use the SmallerAlignment Index (SALI) method which has been proved to be a very fast yetreliable tool. The mathematical definition of SALI is the followingSALI(t) ≡ min(d − , d + ) , (10)where d − and d + are the alignments indices. For computing the SALI we tracksimultaneously the time-evolution of the main orbit and the two deviation vectors.The nature of an orbit can be determined by the numerical value of SALI at theend of the numerical integration t max . In particular, if SALI > − the orbit issaid to be regular, whereas if SALI < − we can thus conclude that the orbit ischaotic. On the other hand, when the value of SALI lies in the interval [10 − , − ]we have the case of a sticky orbit a and further numerical integration is needed toreveal the true character of the orbit.In our numerical integrations the maximum time is 5000 dtu (dimensionless timeunits), corresponding to about 9445.93 Julian years. It should be mentioned, thatorbits which do not escape or collide to Jupiter after a numerical integration of 5000dtu are considered as bounded orbits (regular, sticky or chaotic). a With the term “sticky orbit” we refer to a special type of orbit which behaves as a regular onefor long integration times before it exhibits its true chaotic nature. arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036
Orbital dynamics of the Sun-Jupiter system Our next task is to define appropriate numerical criteria in order to distinguishbetween the above-mentioned types of motion. In order to consider a more realisticapproach, we assume that Jupiter is a finite body, taking into account its meanradius approximately by 66854 km (about 8 . × − dimensionless lengthunits). The numerical integration stops when an orbit reaches the surface of Jupiter,thus producing an orbit leaking in the configuration space. Furthermore, an escapingorbit to the Sun realm must satisfy the condition x < x ( L ) − δ , with δ = 0 . x >x ( L ) + δ , with δ = 0 .
04. At this point, we must clarify that the tolerances δ and δ have been included in the escape criteria in order to avoid the incorrectclassification of the unstable Lyapunov orbits as escaping orbits.The equations of motion (1-2) are numerically integrated by means of a doubleprecision Bulirsch-Stoer integrator, using a numerical routine written in standard FORTRAN 77 , with a fixed time step equal to 10 − . Here it should be noted that theBulirsch-Stoer integrator is both faster and more accurate than a double precisionRunge-Kutta-Fehlberg 7(8) algorithm with Cash/Karp coefficients. Throughout allour computations, the Jacobi constant of Eq. (3) was conserved with fractional ac-curacy of about 10 − , or even better. The Lemaitre’s global regularization methodis applied in the case of collision orbits, when the test particle moves aroundJupiter into a region of radius 10 − .For the numerical integration of the sets of initial conditions of orbits, in alltypes of two-dimensional colour-coded diagrams, that we will be presented in thefollowing section, we needed roughly between 17 hours and 3.5 days of CPU timeon an Intel (cid:114) Quad-Core i7 2.4 GHz PC. Moreover, all graphical illustrations pre-sented in this paper have been created using the latest version 11.2 of the softwareMathematica (cid:114) .
4. Orbit classification
In this section we will perform a thorough analysis of initial conditions of orbitsin all possible Hill region configurations. Parallel to the classification we shall alsorecord the time scale (or time period) of the collision and the time scale of theescapes.In the following subsections we shall present colour-coded diagrams, thus follow-ing the methods also used in Refs. 41, 42. In these diagrams, each pixel is assigneda specific colour according to the particular type of the nature of the orbit. Thesecolour-coded diagrams are, in a way, a modern version of the classical Poincar´e Sur-face of section, where the phase space is a complex mixture of basins b of boundedmotion (regular or chaotic), escape and collision. Our numerical calculations indi-cate that the vast majority of bounded basins correspond to regular orbits, where a b By the term “basin” we refer to a local set of initial conditions which lead to a certain final state(collision, escaping or bounded motion). arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe
Fig. 2. Colour-coded basin diagrams for Case I. (a-left): C = 1 . C = C . Thecolour code is as follows: non-escaping regular orbits (blue), trapped sticky orbits (magenta),trapped chaotic orbits (yellow), collision orbits (cyan). (Colour online only.)Fig. 3. Distribution of the collision time of the orbits for the values of the Jacobi constant ofFig. 2(a-b). Large collision times are identified by darker colour, while the initial conditions of alltypes of bounded orbits are shown in white. (Colour online only.) third integral of motion is present. This additional integral poses new restrictions tothe available phase space and therefore it prevents them from escaping to infinity. Case I:
C > C Our investigation begins with the scenario where the Hill region configurationsconsist of small disks around the secondary (Jupiter). In this case only two typesof motion are possible: collision and bounded motion. The basin diagrams for twovalues of the Jacobi constant are presented in Fig. 2(a-b). The outermost black solidline is the ZVC which is defined as J ( x, y, ˙ x = 0 , ˙ y = 0) = C .In panel (a) of Fig. 2, where C = 1 . P . Such periodic orbits have a reflectionarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system symmetry, with respect to the horizontal x axis. On the other hand, the stabilityisland on the right side of Jupiter is composed of initial conditions that correspondto quasi-periodic orbits around P , travelling in counter-clockwise (prograde) sense,with respect to the rotating frame of reference. It is known that the regular orbitson the left side of Jupiter are much more stable than those on the right side of thesecondary, in relation to the variation of the Jacobi constant C .Panel (b) of Fig. 2 illustrates the orbital structure of the configuration ( x, y )plane when C = C . It is seen that a third thin basin of initial conditions of non-escaping regular orbits emerges near the left boundary of the ZVC. There is nodoubt that the main difference with respect to what we seen in panel (a) of Fig. 2concerns the area at the right side of Jupiter. We observe that the correspondingstability island which was present for C = 1 .
525 has split into pieces, while a trappedchaotic domain appears into the right side of Jupiter. With a much closer look wemay identify, inside the trapped chaotic area, several isolated initial conditions whichcorrespond to sticky orbits with sticky period larger than 5000 dtu. It should beemphasized that in Ref. 42 similar results have been reported, regarding the Pluto-Charon system, for a Jacobi value very close to the critical value C . However inthat case there was no further classification of the bounded orbits into sub-categories(ordered, sticky and chaotic).In the following Fig. 3(a-b), using tones of blue, we show the distribution of colli-sion times on the configuration ( x, y ) space. Light colors correspond to fast collisionorbits, dark colors indicate large collision times, while white colour denote all typesof bounded motion. It is interesting to note that in both energy levels the initialconditions of collision orbits form complicated colour layers (zones). Furthermore,one may observe that the orbits with the highest values of collision times have initialconditions mainly located in the vicinity of the boundaries of bounded basins. Case II: C > C > C The next case considers the second Hill region configurations in which the testparticle is allowed to move between Sun and Jupiter, through the open channelaround L . The orbital structure of the configuration ( x, y ) space is unveiled in Fig.4(a-d), using colour-coded diagrams.Our numerical calculations indicate that in the interval 1 . < C < C ,despite the existence of a bottleneck channel around L (see Fig. 1), there is noevidence of escaping orbits to the Sun domain. This clearly means that, for this rangeof the Jacobi constant values, the channel is very narrow and therefore all chaoticorbits need much more than 5000 dtu to eventually escape from the scatteringregion. In panel (a) of Fig. 4, where C = 1 . mixture of sticky, chaotic and escaping orbits. At theboundaries of the chaotic region and along the symmetry axis, we notice a chain ofsix small stability islands (archipelago), mainly composed by quasi-periodic orbits.It should be noted that except for the existence of the stability archipelago at thearch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe
Fig. 4. Colour-coded basin diagrams for Case II. (a-upper left): C = 1 . C = 1 . C = 1 . C = C . The colour code is as follows:non-escaping regular orbits (blue), trapped sticky orbits (magenta), trapped chaotic orbits (yellow),collision orbits (cyan), escaping orbits to Sun realm (red). (Colour online only.) frontier of the chaotic zone, the structure of the basin diagram is almost identicalto the one presented earlier in panel (b) of Fig. 2. When C = 1 . C = 1 . C = 1 . C = C , one may observe in panel (d) of Fig. 4, that othertypes of secondary resonant quasi-periodic orbits appear, while again about one-third of the configuration ( x, y ) plane is occupied by initial conditions of escapingorbits. From the shapes of several types of basins, shown in Fig. 4(a-b), it becomesmore evident that the parametric evolution of the bounded, collision as well asescape basins in the Sun-Jupiter system, is very different, compared to previousstudies regarding celestial bodies in the Solar System. It deserves mentioning that in our numerical calculations we have followed theapproach used in Ref. 10 and of course in Refs. 41, 42, i.e., we consider that an orbitarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036
Orbital dynamics of the Sun-Jupiter system escapes to the Sun domain if the test particle passes through L , even if its finalstate is different for very long integration time. This assumption is supported bythe fact that just a very small portion of orbits that initially enter the Sun realmreturn inside the scattering region.In Fig. 5(a-d), using tones of blue for collision orbits and a rainbow pallet forescaping orbits, we display the corresponding distributions for collision and escapetimes of the orbits. Once more, light colors correspond to short collision/esacpetimes, dark colors indicate large collision/escape times, while white colour denotesall types of bounded motion. By inspecting the spatial distribution of the vari-ous different ranges of escape times, we are able to associate medium escape timewith the stable manifold of a non-attracting chaotic invariant set, which spreadout throughout the chaotic sea. On the other hand, the largest escape times areassociated to the sticky orbits, surrounding the stability islands.arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe
Fig. 6. Colour-coded basin diagrams for Case III. (a): C = 1 . C = 1 . C =1 . C = 1 . C = 1 . C = C . The colour code is as follows: non-escapingregular orbits (blue), trapped sticky orbits (magenta), trapped chaotic orbits (yellow), collisionorbits (cyan), escaping orbits to Sun realm (red), escaping orbits to the exterior realm (green).(Colour online only.) Case III: C > C > C From the point of view of planetary systems and celestial mechanics, the case C >C > C constitutes the Hill region configurations with the most dynamical interest(the reader can find a detailed discussion about this topic in Ref. 10). When C For C = 1 . C = 1 . x, y )plane. However, the portion of escaping orbits to the exterior realm is almost doublearch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe Fig. 8. Colour-coded basin diagrams for Case IV. (a-upper left): C = C , (b-upper right): C =1 . C = 1 . C = 1 . 4. The colour code is the same as in Fig.6. (Colour online only.) with respect to that of the escaping orbits to the Sun realm. When C = 1 . C , one may observe in panel (f) of Fig. 6 that the initialconditions of orbits that escape to the exterior region dominate the configurationplane by occupying more than 70% of its area. Furthermore, there is a weak presenceof collision and bounded regular orbits around the Jupiter, which are manifestedthrough the existence of small basins.In Fig. 7(a-f) we illustrate the corresponding distribution of the escape andcollision time of orbits on the configuration ( x, y ) space. A closer look at the timescale, given in the appended colour-bar, shows that in most of the cases, more than90% of the initial conditions of the orbits escape from the scattering region, in lessarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system than 10 dtu.When C > C > C the channel around L opens, thus allowing orbits to escapeto the exterior region also from the left side of the primary (Sun). However, sincewe decided to focus our study in the vicinity of the secondary (Jupiter), this escapechannel in the ZVC is not visible and therefore this energy case has limited physicalmeaning in our study. Case IV: C ≤ C The last case under consideration involves the scenario when the test particle canfreely travel all over the configuration ( x, y ) plane with no restrictions of energet-ically forbidden regions. Again, all the different aspects of the numerical approachremain exactly the same as in the previously studied cases.Fig. 8(a-d) reveals the parametric evolution of the orbital structure of the config-uration space, through the colour-coded diagrams. The most important phenomenawhich take place, as the value of the Jacobi constant decreases, are the following:arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe • The area on the ( x, y ) plane covered by initial conditions of orbits thatescape to the exterior region grows rapidly and for C < . • The area of the basin composed of initial conditions of orbits that escapeto the Sun realm constantly decreases. • The portion of the initial conditions that lead to collision with Jupiterseems to be almost unperturbed by the change on the value of the Jacobiconstant. In particular, a thin collision basin is always present, which seemsto survive even at very low values of C , or equivalently, at extremely highlevels of the total orbital energy. • We find no numerical evidence of non-escaping regular orbits for C < . C ≤ C ) the collision to the secondary is very fast, becausethe vast majority of the corresponding initial conditions of the orbits have collisionrates lower than about 0.5 dtu. A summarized analysis of the numerical results It would be very informative for the reader to observe a summarized monitoring ofthe percentage of orbits pertaining to the different categories, introduced in section3. In Fig. 10, we show a diagram with the parametric evolution of all types ofpercentages, as a function of the Jacobi constant C . From this figure, we can seethat at very high values of C , regular bounded motion is the most populated typeof motion, occupying more than 70% of the configuration space. However, as thevalue of the Jacobi constant decreases the percentage of non-escaping regular orbitsdiminishes until C = C , where at this point, the magnitude of the slope changesand the percentage reduction continues much slower. Bounded regular motion ispossible up to about C = 1 . C there is no numericalevidence of this kind of motion. For most of the investigated range of the Jacobivalues chaotic motion corresponds to extremely low percentages (less than 0.5%),the only observed peak occurs at exactly C = C and it is equal to about 22%. Asthe channel around L opens for C < C the percentage of escaping orbits to theSun realm increases until C = C , while for lower values of C the trend is reversed.For C < C the second channel around the Lagrange point L also opens andthe percentage of orbits escaping to the exterior realm displays a rapid increase.For C < . x, y ) plane,while for extremely low values of the Jacobi constant ( C < C ) they occupy morethan 90% of the configuration space. The rate of collision orbits gets reduced assoon as C < C and tends asymptotically to zero. However, our analysis indicatesarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system x, y ) plane, in terms of the Jacobi constant C . The vertical dashed black linesindicate the four critical values of C . (Colour online only.)Fig. 11. (a-left): Orbital structure of the ( x, C ) plane. (b-right): The distribution of the corre-sponding collision and escape time of the orbits. The colour codes are the same as in Figs. 6 and7, respectively. (Colour online only.) that collision motion is possible for all tested values of C , even at extremely lowpercentages (less than 1%).Taking into account the above-mentioned analysis, we may say that the evolutionof the percentages of each basin in the Sun-Jupiter system is very similar to thatobserved for the Pluto-Charon system in Ref. 42. The main difference concerns thearch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe Fig. 12. Evolution of (a-left): the percentages of all types of orbits and (b-right): the averagelogarithmic collision and escape time ( < log ( t ) > ) of orbits on the ( x, C )-plane, as a function ofthe value of the Jacobi constant C . (Colour online only.) interval C < C , where for the Pluto-Charon system we had used more complicatedescape criteria (more escape sectors), thus leading to different percentages.The colour-coded diagrams in the configuration ( x, y ) space provide enoughinformation about the phase space mixing. However, such analysis is performedonly for fixed values of the Jacobi constant and for orbits that traverse the surfaceof section, either progradely or retrogradely. In order to surmount these barriers,Ref. 19, introduced a new plane of representation that can supply useful informationabout the classification of the orbits by using the section y = ˙ x = 0, ˙ y > 0. Withthis new approach the Jacobi integral can be used as an independent variable, andtherefore, the orbital structure of the Sun-Jupiter system can be monitored using acontinuous spectrum of values of C . In panel (a) of Fig. 11 we present the orbitalstructure of the ( x, C ) plane when x ∈ [ x ( L ) , x ( L )] and C ∈ [1 . , . f ( x, C ) = J ( x, y = 0 , ˙ x = 0 , ˙ y = 0) = C, (11)while the critical values of C are pointed by horizontal dashed black lines.The two stability islands, corresponding to prograde and retrograde motionaround Jupiter are now visible. Once more it is manifest that choosing smallervalues of the Jacobi constant C , induces a larger number of orbits escaping to theexterior domain through L , while larger values of C lead to an overpopulation ofbounded orbits. Just below the right stability island, which ends at about C = C ,there is a fractal mixture of initial conditions corresponding to both types of es-arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system x, C ) plane. Initial conditions with the samefinal state in both classical Newtonian and PN dynamics are shown in green, while converginginitial conditions, with different final states, are shown in red. (Colour online only.) caping orbits. The colour-coded diagram shows us exactly how the fractality ofthe several basin boundaries strongly depends not only on the Jacobi constant butalso on the spatial variable. In particular, one can observe a very interesting phe-nomenon. It is seen that the fractality of the basin boundaries, which is related tothe unpredictability, migrates from the upper right side to the lower left side of thesecondary (Jupiter), for low values of the Jacobi constant (i.e. high values of thetotal orbital energy).The parametric evolution of the percentages of all types of orbits on the ( x, C )plane is (see panel (a) of Fig. 12), in general terms, very similar to that discussedearlier in Fig. 10, for the configuration ( x, y ) plane. The combined analysis of theorbit classification, in both types of planes, strongly suggests that at very highvalues of the Jacobi constant C bounded regular motion completely dominate alltypes of planes. At very low values of C on the other hand, escaping motion tothe exterior realm is, by far, the mots populated type of motion. It would be alsointeresting to shed some light to the collision as the well as to the escape time of theorbits. The evolution of the average logarithmic value of the collision and escapetime ( < log ( t ) > ) of the orbits on the ( x, C ) plane, as a function of the valueof the Jacobi constant C , is given in panel (b) of Fig. 12. It is seen there that inmost of the cases the average escape time of the orbits to the exterior realm is lowerthan the average escape time of the orbits to the Sun realm. It is interesting to notethat the peak of the escape time of the Sun realm is about 63 time units, when C = 1 . C = 1 . C = 1 . (cid:15) = 0) and the PNarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe dynamics ( (cid:15) = 1). Our numerical analysis strongly suggests that in general termsthe influence of the PN terms is rather weak, however it does affect the final stateof the orbits. To prove this we reclassified the initial conditions of the orbits on the( x, C ) plane (see Fig. 11), by setting (cid:15) = 0. The corresponding colour-coded diagram,with the classification of the orbits, looks almost identical to that shown in panel(a) of Fig. 11. More precisely, only about 12% of the total initial conditions displaysa different final state between CN and PN. In order to observe the differences on thefinal state of the orbits, due to the inclusion or not of the PN terms, we have to zoomin and examine local regions of the ( x, C ) plane. In Fig. 13 we see a local regionof the ( x, C ) plane, with the highest portion of converging initial conditions. Greencolour indicated initial conditions with the same final state in both CN and PN,while the divergent initial conditions are shown in red. It is seen that both typesof initial conditions of orbits are not randomly distributed on the ( x, C ) plane,as we can distinguish basins of converging and not converging initial conditions.Additional computations reveal that inside the local region, shown in Fig. 13, allpossible divergent initial conditions exist. In particular, since there are five mainfinal states (trapped chaotic, non-escaping regular, collision, escaping to Sun realmand escaping to the exterior region) and two different dynamical cases (CN and PN)the corresponding total number of permutations is N = 5! / (5 − 5. Discussion The orbital dynamics of a small body (e.g., a spacecraft, comet, meteor or aster-oid) in the presence of the Sun-Jupiter system has been numerically investigated.By using the PN equations of motion for the planar circular restricted three-bodyproblem, we have performed an orbital classification in a scattering region aroundJupiter. To do so, we have determined the position of the Lagrange points and thecorresponding values of the Jacobi constant at these points. After numerically inte-grating several sets of initial conditions, for all possible Hill region configurations,we managed to classify the orbits into four main categories: bounded orbits aroundJupiter, escaping orbits to Sun realm, escaping orbits to the exterior realm andorbits that lead to collisions with Jupiter. Furthermore, the SALI chaos indicatorhas been used in order to further classify bounded orbits into three sub-categories:regular orbits, sticky orbits, and chaotic orbits.Despite the fact that the PN correction terms are modulated by a factor of10 − , it is clear from our results that, in general, the differences on the Newtonianand PN Lagrange points are non-negligible, mainly for the triangular points L and L (approx. 1 km). Our result is in agreement with previous findings for the post-Newtonian collinear points reported in Ref. 40, where the authors have found thatthe correction of distance for the triangular points is 923 m. Consequently, the finalstates of the orbits in both approaches can be significantly different. Our numericalarch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system C x CN PN1.50708 0.939075 0 -11.50701 0.935000 -1 01.50702 0.935250 0 11.50702 0.937350 0 21.50705 0.935100 1 01.50703 0.937355 2 01.50708 0.937300 0 91.50709 0.937350 9 01.50811 0.938225 1 21.50804 0.938450 2 11.50768 0.937500 1 91.50811 0.939150 2 91.50700 0.937750 9 11.50700 0.935375 9 21.50839 0.936675 1 -11.50709 0.937775 2 -11.50881 0.938725 -1 11.50869 0.939200 -1 21.50711 0.937400 9 -11.50717 0.937800 -1 9analysis strongly suggests that a refined version of the Sun-Jupiter system, thatincludes general relativistic correction terms, could certainly update most of ourcurrent knowledge about the dynamics of this particular planetary system.It is important to note that along the paper we neglect the perturbations pro-duced by other planets. This approximation is based on the fact that despite therelativistic contribution, due to the Sun, is 6 orders of magnitude smaller than theone of the Newtonian terms, this contribution is 4 orders of magnitude larger thanthe effect due to the presence of the second most-massive planet in the Solar System,i.e. Saturn.As far as we know, this is the first time that the nature of motion, in the vicinityof Jupiter, is explored in such a thorough and systematic manner, through the orbitclassification using the two-dimensional colour coded diagrams. Therefore we mayclaim that our paper adds considerable new information to the field of planetaryand celestial dynamics.We hope that the current numerical results to be useful in the active field of theorbital dynamics of the PN version of the restricted three-body problem. Takinginto account that our present outcomes are positive, as well as encouraging, it is inour future plans to expand our investigation into three dimensions, thus revealingthe orbital content of the full six-dimensional phase space.arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 E.E. Zotos & F.L. Dubeibe Acknowledgments One of the authors (FLD) gratefully acknowledges the financial support provided byUniversidad de los Llanos and COLCIENCIAS, Colombia, under Grant No. 8840.We would like to express our warmest thanks to the anonymous referee for thecareful reading of the manuscript and for all the apt suggestions and commentswhich allowed us to improve both the quality and the clarity of our paper. Appendix A. Relationship between physical and dimensionlessunits The planar circular restricted three body problem has three natural scales associatedto the three fundamental mechanical quantities (time, distance, and mass): thedistance between the primaries a , the total mass of the system M = m + m , andthe angular velocity of their orbital motion ω . In the canonical system of units,introduced in section 2.1 all this three quantities are dimensionless and normalizedto 1, i.e. a = M = 1 and ω ≈ 1. Therefore it implies that in the resulting system ofunits the kinematic quantities are also dimensionless.Taking into account that the average distance between Jupiter and the Sunis 5.20336301 AU, and the orbital period of Jupiter is about 11.8701 years, theconversion of astronomical units and years to the canonical units are as follows: 1AU = 1/5.20336301 and 1 year = 2 π/ . c = 63232 . c = 22945 . 23, 33 As can be noted, we have assumed that the post-Newtonian contribution tothe angular velocity is too small, i.e ω = 1 + ω /c ≈ 1. In order to prove thisassumption, let us start by considering that unlike the classical Newtonian system,the post-Newtonian angular velocity depends on the mass parameter µ and thespeed of light c , i.e. ω = 1 + µ (1 − µ ) − c , (A.1)where the mass parameter, for the Sun-Jupiter system, is given by µ =0 . v/c = 0 . v = ωa ,such that ω/c = v/c and from Eq. (A.1) vc = 1 c (cid:18) µ (1 − µ ) − c (cid:19) , (A.2)replacing the numerical values for v/c and µ , and solving for c we get c =22945 . ω ≈ Appendix B. Coefficients for the Lagrange points The exact numerical coefficients entering equations (8) and (9) are the following:arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system C = 1 . C = − . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = − . C = 2 . C = − . C = 14 . C = − . C = 20 . C = − . C = 4 . C = − . C = 0 . C = − . C = 0 . C = − . C = 0 . C = 0 . C = 0 . C = − . C = 0 . C = 1 . C = − . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = − . C = 6 . C = − . C = 29 . C = − . C = 42 . C = − . C = 8 . C = 0 . C = 0 . C = − . C = − . C = − . C = 0 . C = 0 . C = 0 . C = − . C = 0 . E.E. Zotos & F.L. Dubeibe C = − . C = − . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . C = 0 . 0; ˜ C = − . × − ;˜ C = 8 . × − ; ˜ C = − . C = 2 . × − ; ˜ C = − . C = 0 . C = − . C = 0 . C = − . C = 0 . T = 0 . T = − . T = 2 . × − ; ˜ T = 0 . T = 8 . × − ; ˜ T = − . × − ; T = 0 . T = − . T = 2 . × − ; ˜ T = − . T = − . × − ; ˜ T = 8 . × − ; References 1. Aguirre J., Viana R.L., Sanju´an M.A.F., Phys. Rev. E (2001) 066208.2. Arnold V.I., Russian Mathematical Surveys (1963) 85.3. Batygin K., Laughlin G., ApJ (2008) 1207.4. Bhatnagar K.B., Hallan P.P., Celest. Mech. Dyn. Astron. (1997) 271.5. Binney J., Tremaine S., Galactic Dynamics (Princeton Univ. Press, Princeton, USA,2008).6. Brumberg V.A., Relativistic Celestial Mechanics (Moscow, 1972).7. Brumberg V.A., Essential Relativistic Celestial Mechanics (Hilger, Bristol, 1991).8. Contopoulos G., In Memoriam D. Eginitis , 159. D. Kotsakis, Ed., (Athens, 1976).9. Darriba L.A., Maffione N.P., Cincotta P.M., Giordano C.M., International Journal ofBifurcation and Chaos (2012) 1230033.10. de Assis S.C., Terra M.O., Celest. Mech. Dyn. Astron. (2014) 105.11. Douskos C.N., Perdios E.A., Celest. Mech. Dyn. Astron. (2002) 317.12. Dubeibe F.L., Lora-Clavijo F.D., Guillermo A.G., A&SS (2017) 97.13. Einstein A., Infeld L., Hoffmann B., Annals of Mathematics (1938) p. 65-100. arch 21, 2018 0:30 WSPC/INSTRUCTION FILE IJMPD˙1850036 Orbital dynamics of the Sun-Jupiter system 14. Fitzpatrick R., An introduction to celestial mechanics (Cambridge University Press,2012).15. G´omez G., Llibre J., Mart´ınez R., Sim´o, C., Dynamics and Mission Design Near Li-bration Points, Volume I: Fundamentals: The Case of Collinear Libration Points (WorldScientific, Singapore, 2001).16. G´omez G., Sim´o C., Llibre J., Mart´nez R., Dynamics and Mission Design NearLibration Points, Volume II: Fundamentals: The Case of Triangular Libration Points (World Scientific, Singapore, 2001).17. G´omez G., Jorba A., Sim´o C., Masdemont J., Dynamics and Mission Design NearLibration Points, Volume III: Advanced Methods for Collinear Points (World Scientific,Singapore, 2001).18. G´omez G., Jorba A., Sim´o C., Masdemont J., Dynamics and Mission Design NearLibration Points, Volume IV: Advanced Methods for Triangular Points (World Scientific,Singapore, 2001).19. H´enon M., A&A (1969) 223.20. Huang G., Wu X., Physical Review D (2014) 124034.21. Krefetz E., AJ (1967) 471.22. Laskar J., Large scale chaos and marginal stability in the solar system. In Chaos inGravitational N-Body Systems pp. 115-162, (Springer, Netherlands, 1996).23. Lhotka C., Celletti A., Icarus (2015) 249.24. Lyapunov A.M., Ann. Fac. Sci. Toulouse (1907) 203.25. Maindl T.I., Dvorak R., A&A (1994) 335.26. Musielak Z.E., Quarles B., Reports on Progress in Physics (2014) 065901.27. Nagler J., Phys. Rev. E (2004) 066218.28. Nagler J., Phys. Rev. E (2005) 026227.29. NASA Space Science Data Coordinated Archive: http://nssdc.gsfc.nasa.gov/ 30. Poincar´e H., New methods of celestial mechanics Vol 1, (Paris Gauthier-Villars, 1982).31. Poisson E., Will C.M., 2014, Gravity: Newtonian, post-Newtonian, Relativistic, Cam-bridge University Press32. Press H.P., Teukolsky S.A, Vetterling W.T., Flannery B.P., Numerical Recipes inFORTRAN 77 Nonlinear Analysis: Theory,Methods & Applications (2001) 3413.34. Sim´o C., Stuchi T., Physica D (2000) 1.35. Skokos C., Journal of Physics A (2001) 10029.36. Stuchi T.J., Yokohama T., Corrˆea A.A., Sol´orzano R.H., Sanchez D.M., WinterS.M.G., Winter O.C, Adv. Space Res (2008) 1715.37. Szebehely V., Theory of Orbits:the restricted problem of three bodies (Academic Press,New York, 1967)38. Will C.M., Proceedings of the National Academy of Sciences (2011) 5945.39. Wolfram S., The Mathematica Book Fifth Edition. (Wolfram Media, Champaign,2003).40. Yamada K., Asada H., Phys. Rev. D (2012) 124029.41. Zotos E.E., A&SS (2015) 4.42. Zotos E.E., A&SS (2015) 7.43. Zotos E.E., A&SS361