Conic-Helical Motion in the Three-Body Problem: Star-Planet-Moon Systems and Relativistic Effects in Binary-Star-Planet Systems

The subject of the study is the following variation of the general three-body problem in astrophysics. There are two relatively heavy bodies rotating around their barycentre, one of the bodies being significantly heavier than the other one. In addition, there is a third, much lighter body whose orbit is relatively close to the lighter of the two heavy bodies, the orbit being not in the plane of the rotation of the two heavy bodies, so that the motion of the entire three-body system is really three-dimensional. One particular kind of such systems are star-planetmoon systems. Another particular kind of such systems are binary-star-planet systems (i.e., a planet around a binary star). For brevity, we call such 3-dimensional 3-body systems “3D3BS”.


Introduction
The subject of the study is the following variation of the general three-body problem in astrophysics. There are two relatively heavy bodies rotating around their barycentre, one of the bodies being significantly heavier than the other one. In addition, there is a third, much lighter body whose orbit is relatively close to the lighter of the two heavy bodies, the orbit being not in the plane of the rotation of the two heavy bodies, so that the motion of the entire three-body system is really three-dimensional. One particular kind of such systems are star-planet-moon systems. Another particular kind of such systems are binary-star-planet systems (i.e., a planet around a binary star). For brevity, we call such 3-dimensional 3-body systems "3D3BS".
Conic-helical orbits in 3D3BS-the astrophysical analogue of 1-electron Rydberg quasi molecules studied previously by Oks [1,2] are of interest for several reasons. First, the possibility of stable (or metastable) conic-helical orbits in 3D3BS is a fundamental problem in its own right, as a relatively new chapter of the centuries-old classical three-body problem. Second, one of its applications-namely, to binarystar-planet systems-is especially significant for the search of an extraterrestrial life. Indeed, according to estimates [3][4][5], approximately 50% of binary stars could support habitable terrestrial planets within stable orbital ranges.
The result in paper I was obtained by applying the standard general analytical method for systems that can be separated into rapid and slow subsystems. The method was applicable because there are ranges of parameters, specified in paper I, where the primary frequency Ω of the conic-helical motion of the planet in binary-star-planet systems is much greater than the Kepler frequency ω of the star's revolution around their barycentre. Even with the allowance for the star's rotation, the trajectory of the planet is still conic-helical. The plane of the quasicircular planetary orbit undergoes relatively small oscillations along the rotating interstellar axis and the radius of the planetary orbit also undergoes relatively small oscillations. The last, but not least: in paper I, there were also obtained positive results concerning the transitability (and thus, detectability) of such planets.
In the present paper, first, we extend that study to star-planetmoon systems-also in frames of the nonrelativistic classical mechanics. We complement analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, by exact simulations showing that the moon can have practically stable conic helical orbits around the planet, the average plane of the orbits being perpendicular to the axis connecting the planet and the star. Second, we extend that study to the relativistic classical mechanics. We show that relativistic effects can become significant in conic-helical orbits of a planet around a binary star for the situations where the mass of the planet is relatively small (such planets are socalled planetoids). Again, we complement analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, by exact simulations showing that the planet can have relatively stable conic helical orbits around the lighter star, the average plane of the orbits being perpendicular to the interstellar axis.
The paper is organized as follows. In Section 2 we apply those nonrelativistic results to star-planet-moon systems and demonstrate the stability of the moon in conic-helical orbits. In Section 3 we present the relativistic study of OBSS and of the stability of planetoids in conic helical orbits. In Section 4 we summarize the results into conclusions.

Application of the Non-Relativistic Results to Star-Planet-Moon Systems
We consider, as an example, an Earth-like planet around a Sun-like star, so that the ratio of the masses / ′ of the planet and the star is 3 × 10 -6 . The separation between the planet and the star is R=1 AU. The planet has a moon. Figure 1 shows the ratio Ω/ω of the frequency Ω of the moon revolution around the planet-star axis to the frequency ω of the planet rotation around the star versus the scaled projection w=z/R of the average plane of the moon's orbit on the planet-star axis. We use the cylindrical coordinates with the z-axis along the planet-star axis in the direction from the planet to the star. It is seen that in the range of 0<w<10 -4 presented in Figure 1, the ratio Ω/ω >> 1, so that the separation of rapid and slow subsystems is justified. Figure 2 presents the scaled amplitude δw 0 of the oscillations of the moon's orbit along the planet-star axis versus the scaled coordinate w=z/R. Figure 3 similarly shows the scaled amplitude δv 0 of the oscillations of the scaled radius v=ρ/R of the moon's orbit in the plane perpendicular to the planet-star axis versus the scaled coordinate w. It is seen that these oscillations are small.
In addition to the above analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, we also performed exact simulations of the moon's motion, as the moon revolves around the planet-star axis while this axis rotates with the Kepler frequency corresponding to the star-planet twobody problem. Figure 4 shows the dependence of m, the scaled (by R) y-coordinate in the Cartesian righthanded xyz-system with the z-axis being the planet-star axis, on the scaled time τ defined in Equation (A.4) from Appendix, for the initial value of the scaled projection of the average plane of the moon's orbit on the planet-star axis w 0 =10 -4 , i.e., 150000 km, which is of the same order of magnitude as the distance between the Earth and the Earth's Moon.       It is clearly seen from Figures 4-7 that the conic-helical orbit of the moon in the star-planet-moon system is stable. Thus, under the influence of a relatively distant star, the plane of the moon's orbit could orient itself to be perpendicular to the planet-star axis. The trajectory of the moon becomes conic-helical.

Relativistic Effects in Binary-Star-Planet Systems
The relativistic force acting on the planet (scaled by its mass, i.e., the force divided by the mass of the planet) is given by the formula where and V are the acceleration and the velocity of the planet, respectively, in the inertial reference frame, and ( ) 2 The additional acceleration given by Equation (A.1) from Appendix can be substituted into (1) as the last factor, with V as the time derivative of the radius-vector in the rotating reference frame.
where the z-axis is along the interstellar axis in the direction from the lighter star of the mass to the heavier star of the mass ′. After the calculation, we obtain the following values for the components of the force: Since ω << Ω, the dominant term is the z-projection of the force (given by the third equation in (3)). By analogy with the non-relativistic case from paper I, we derive the small oscillations of the scaled axial w and radial v coordinates about their equilibrium values, due to the dominant force term (compare to Equation (A.3)): where v 0 is the equilibrium value of ρ/R and the tilde above ω and Ω means the scaling by multiplying by (R 3 /Z) 1/2 . We can now find the conditions, under which the amplitude of the oscillations is small while the condition ω << Ω stays valid. In the relativistic case, the Hamiltonian, divided by the mass m of the planet, h=H/m for the 3D3BS case is given by where G is the gravitational constant. In a circular state, the components of the linear momentum in the cylindrical coordinates p z =p ρ =0 and |p φ |/m=const=L. Using the scaling we write the scaled energy for the circular state: For the relativistic motion one has: Using the scaled radius v=ρ/R and the first formula in (6), we can find the speed of the planet in the relativistic case in the units of the speed of light: The equilibrium values for p and r can be obtained from differentiating (7) with respect to w and v. The first differentiation yields the same relation as in the non-relativistic case, so the equilibrium value of p is the squared right-hand side of the second line in (4). The second differentiation (with respect to v), with the later substitution of the equilibrium value of v (or p), yields the equilibrium value of ℓ, which is related to r by ℓ 2 =α/r 2 , as mentioned above. Using the substitution [13] 1/3 (11) which significantly simplifies the formulas in the two-Coulombcenter problem, we find the speed of the planet in the circular state (in the units of the speed of light): From (8) with the substituted value of V=Ωρ, by using the first relation in (6), the last relation in (9), and the formula ℓ 2 =α/r 2 , we find Ω = β 2 c α 1 -R pr (13) and its scaled counterpart α Ω β  2 R = c 1 -Z pr (14) As given in (5), and measuring now the mass of the star μ in the units of the mass of the Sun (in distinction to the nonrelativistic case in Section 2, where the star masses were measured in units of the mass of the planet), we get where ⊙ =1.989 × 10 33 g. Substituting (15) into (14), we obtain where the quantity s=GM ⊙ /c 2 (which is one half of the Schwarzschild radius of the Sun and is approximately equal to 147700 cm). As defined above, α 1/2 =Z/(cL), and substituting the scaling relation for L from (6), ℓ 2 =α/r 2 and (15), we have We substitute (17) into (12), and then substitute the resulting equation for β and the solution for α into (16), obtaining the equation for the scaled frequency of the revolution of the planet that only depends on the ratio R/(μs) and γ (or w-see (11)). Finally, from (4), and because the scaled Kepler frequency is ῶ=(1 + b) 1/2 the amplitude of the small oscillations about the equilibrium on the w-axis and v-axis are For deriving more explicit expressions for the amplitude of these small oscillations in the relativistic case, we use Equation (16) for the scaled frequency of the revolution of the planet, and Equation (12) for β, and Equation (A.4) for the scaled eigenfrequencies ω + and ω -, and the equilibrium value for p=v 0 2 (v 0 being given in Equation (4)), and the following equilibrium value of r: As a result, we arrive to the following final formula for the amplitudes of the small oscillations of the planet on the w-axis and the v-axis in the relativistic case: In turns out that the amplitude in the relativistic case is the same as in the non-relativistic case.
We checked that the ratio of the frequencies of the revolution of the planet Ω and the Kepler frequency ω of the revolution of the stars about their barycentre is much greater than 1. By using the previously found values of Ω and ω and substituting the equilibrium values of p, r, α and β, we derive the formula for the ratio k of the frequencies depending on the axial coordinate for the given interstellar distance R and the star masses μ and μ': From Equation (21) it can be found out that k in the relativistic case is equal to the nonrelativistic k NR multiplied by (1-β 2 ) 1/2 . Figure 8 shows the ratio k/k NR versus the scaled radius of the orbit v, for the masses of the stars μ=1 and μ'=100 (in the units of the mass of the Sun) and the interstellar distance R=100 a.u. It is seen that the relativistic effects become significant when v ~ 10 -9 or smaller.
The range of v=(1 ÷ 2) × 10 -9 for R=100 a.u. corresponds to the range of the radius of the planetary orbit ρ=(15 ÷ 30) km. Since the radius of the planet should be smaller than ρ, in this case it should be a planetoid. Figure 9 shows the dependence k(v) for the example where the mass of the lighter star μ=1 and the heavier star μ'=100 (in the units of the mass of the Sun) and for the interstellar distance R=100 a.u. (wide binary system). It is seen that the ratio of the frequencies is greater than 10 12 for v < 2 × 10 -9 , which means that such a system would be stable for a very long time.
In Figure 10 we present the plot of the scaled amplitude δ 0 of the axial oscillations of the planetary orbit from Equation 20 versus the scaled radial coordinate v for the mass of the lighter star μ=1 and the heavier star μ'=100, in the units of the mass of the Sun. It is seen that the amplitude of the oscillations is really much smaller-about quadrillion (!) times smaller-than the interstellar distance.
In addition to the above analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, we also performed exact simulations of the planetary motion, as the planet revolves around the interstellar axis while this axis rotates with the Kepler frequency corresponding to the binary star twobody problem. Figure 11 shows the dependence of m, the scaled (by R) y-coordinate of the planet in the Cartesian xyz-system with z-axis being the interstellar axis, on the scaled time τ defined in Equation (A.4) from Appendix, for the initial value of the scaled projection of the average plane of the planetary orbit on the interstellar axis w 0 =10 -25 , i.e., practically coinciding with the position of the lighter star in the binary system. Figure 12 similarly presents the dependence of n, the scaled (by R) x-coordinate of the planet in the above-mentioned Cartesian xyzsystem, on the scaled time τ defined in Equation (A.4), for the same initial conditions. Figure 13 shows the dependence of w, the scaled z-coordinate of the planet on scaled time τ defined in Equation (A.4), for the same initial conditions. We remind that that the z-axis is along the interstellar axis.
It is seen from Figures 11-13 that the conic-helical orbit of the planet in the binary-star-planet system is relatively stable. Thus, under the influence of the heavier star in the binary system, the plane of the planetary orbit, being at the lighter star, could orient itself to be perpendicular to the interstellar axis. The trajectory of the planet becomes conic-helical.

Conclusion
We extended the results of paper I to star-planet-moon systems-in     frames of the nonrelativistic classical mechanics. We complemented analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, by exact simulations. The exact simulations demonstrated that the moon can have practically stable conic-helical orbits around the planet, the average plane of the orbits being perpendicular to the axis connecting the planet and the star. It is counterintuitive that under the influence of a relatively distant star, the plane of the moon's orbit could orient itself to be perpendicular to the planet-star axis and the orbit becomes conic-helical.
We also extended the nonrelativistic study from paper I to the relativistic classical mechanics (and also corrected some printing errors from paper I). We showed that relativistic effects in binarystar-planet systems can become significant for the situations where the mass of the planet is relatively small (such planets are so-called planetoids). In this case, the conic helical trajectory of the planetoid encloses the star of the lighter mass, the average plane of the planetary orbit practically coinciding with the position of the star of the lighter mass. It is counterintuitive that a distant heavier star can cause such a modification of the planetoid orbit around the lighter star. We showed that in this case the frequency of the rapid subsystem (i.e., the frequency of the planetoid revolution around the lighter star) exceeds the frequency of the slow subsystem (i.e., the Kepler frequency of the star's rotation around their barycentre) by more than a trillion times. This means that the standard analytical method of separating rapid and slow subsystems, which we used, is well justified and that this configuration will remain stable for a very long time [14,15].
Finally, we complemented analytical results, obtained by the standard accurate approximate method based on the separation of rapid and slow subsystems, by exact simulations. The simulations confirmed that the planet can have relatively stable conic-helical orbits around the lighter star, the average plane of the orbits being perpendicular to the interstellar axis.

Appendix: corrections of misprints/errors in some formulas from paper I and new results based on the corrected formulas
In this Appendix we present corrections to three formulas from paper I. These are three formulas that are used in the main text of the present paper. Because of misprints/errors in these formulas in paper I, it is impossible to simply refer in the main text of the present paper to the corresponding formulas from paper I, which is why it is necessary to present here the corrected formulas numbered below as (A.1), (A.2), and (A.3).
The analysis in paper I was performed in a rotating frame. The frame was rotating at the same frequency as the Kepler frequency ω of the stars rotation. In this frame, there is an additional force (see, e.g., Landau & Lifshitz 1960;Goldstein 1980): [14,15] ( ) The order of magnitude of the velocity v in Equation (A.1) is related to the average radius ρ of the planetary orbit as v ~ Ωρ. Here Ω is the main frequency of the planetary motion. Therefore, Equation (A.1) can be approximated as where is the unit vector along the z-axis, which is the interstellar axis. The motion in the zρ-plane corresponded to a 2D-oscillator driven by the force given by Equation (A.2). In the scaled coordinates w = z/R, v = ρ/R, the solution for the oscillation amplitudes is where, are the scaled (dimensionless) time, as well as the scaled frequencies ω and Ω, respectively. The similarly scaled eigenfrequencies ω + and ωof the above two-dimensional oscillator are The equilibrium value of the scaled radial coordinate 0 ( , ) in Equations (A.3, A.5) was defined in Equation (4).
Based on the corrected formula (A.3) for the small oscillations of this driven two-dimensional oscillator, below we present new illustrations demonstrating the ranges of the stability of the conic-helical orbits for various parameters of this nonrelativistic 3-body problem. Figure 14 shows the scaled amplitude δw 0 =2v 0 (w, b) ω s f p (cos 2α) / |ω + 2 -ω -2 | of the oscillations of the planetary orbit along the interstellar axis versus the scaled projection w = z/R of the average plane of the planetary orbit on  the interstellar axis, for three values of the ratio b= ′/ of the stellar masses. It is seen that δw 0 << 1 (i.e., the amplitude of the oscillations is much smaller than the interstellar distance) for b greater or of the order of 10. For the validity condition Ω >> ω of this result to be satisfied with a large margin of "safety", the average plane of the planetary orbit should be very close to the star of the smaller mass (and the closer it is in the range of w < 0.05, the smaller is the amplitude δw 0 , as seen from Figure  14). As long as the primary frequency Ω of the planet revolution about the interstellar axis exceeds by many orders of magnitude the Kepler frequency ω of stars rotation, the conic-helical planetary orbit would remain stable for a very long time. (Rigorously speaking, the planet is in a metastable state, which is an analogy with atomic/molecular systems: they have metastable states living by many orders of magnitude longer than other states of the system.) Figure 15 shows the scaled amplitude δv = 2v 0 (w, b) ω s f p (sin 2α) / |ω + 2 -ω -2 | of the oscillations of the radius of the planetary orbit versus the scaled projection w = z/R of the average plane of the planetary orbit on the interstellar axis, for three values of the ratio b= ′/ of the stellar masses. It is seen that for the range of b greater or of the order of 10 (required for having δw 0 << 1), we have also δv << 1.