Gravitationally quantized orbits in the solar system: computations based on the global polytropic model

The so-called"global polytropic model"is based on the assumption of hydrostatic equilibrium for the solar system, or for a planet's system of statellites (like the jovian system), described by the Lane-Emden differential equation. A polytropic sphere of polytropic index $n$ and radius $R_1$ represents the central component $S_1$ (Sun or planet) of a polytropic configuration with further components the polytropic spherical shells $S_2$, $S_3$, ..., defined by the pairs of radii $(R_1,\,R_2)$, $(R_2,\,R_3)$, ..., respectively. $R_1,\,R_2,\,R_3,\, ...$, are the roots of the real part $\mathrm{Re}(\theta)$ of the complex Lane-Emden function $\theta$. Each polytropic shell is assumed to be an appropriate place for a planet, or a planet's satellite, to be"born"and"Live". This scenario has been studied numerically for the cases of the solar and the jovian systems. In the present paper, the Lane-Emden differential equation is solved numerically in the complex plane by using the Fortran code DCRKF54 (modified Runge-Kutta-Fehlberg code of fourth and fifth order for solving initial value problems in the complex plane along complex paths). We include in our numerical study some trans-Neptunian objects.


Introduction
In this study, we approach the issue on "gravitational quantization of orbits" in the solar system, or in systems of satellites of planets, by exclusively considering laws of classical mechanics. In particular, we take as basis of our treatment the equations of hydrostatic equilibrium for a nondistorted star or planet. These equations yield the well-known Lane-Emden differential equation, which is solved in the complex plane by using the so-called "complex plane strategy" (Sec. 2), developed by the first author for numerical treatment of certain astrophysical problems (see e.g. [1]; see also [2]).
Classical mechanics is also used by some authors for treating this issue (see e.g. in [3] the so-called "vibrating membrane model"). Mostly, several investigators use ideas arising within the framework(s) of scale relativity as an extension of Einstein's principle of relativity [4]; relativity theory regarding the finite propagation speed of gravitational interaction [5]; and quantum mechanics, like appropriate Bohr-Sommerfeld discretization ( [6], [7]), or Schrödinger-type equations [8].

The Lane-Emden Equation in the Framework of the Complex-plane Strategy
The so-called "complex-plane strategy" (CPS) proposes and applies numerical integration of "ordinary differential equations" (ODE, ODEs) in the complex plane, either along an interval I r ⊂ R when the independent variable r is real, or along a contour C ⊂ C when r is complex. Integrating in C is necessary when the "initial value problem" (IVP, IVPs) under consideration is defined on ODEs: (i) suffering from singularities and/or indeterminate forms in R, and/or (ii) involving terms that become undefined in R when the independent variable r exceeds a particular value. A detailed review of CPS is given in [9] (Sec. 3.1). The equations of hydrostatic equilibrium for a nondistorted star are written as where P (r) is the pressure, m(r) the mass inside a sphere of radius r, and ρ(r) the density. For the polytropic models we use the "polytropic EOS" ( [10], Chapter IV, Eq. (1)) and the "normalization equations" ( [10], Chapter IV, Eqs. (8a), (10a), respectively) where K is the polytropic constant, Γ the adiabatic index defined by Γ = 1 + (1/n), n ∈ [0, 5) ⊂ R the polytropic index, λ the polytropic unit of density set equal to the central density of the star, λ = ρ(r = 0), and α the polytropic unit of length set equal to ( [10], Chapter IV, Eq. (10b)) Thus θ n is the density measured in such units, so-called "classical polytropic units" (abbreviated "cpu"), and ξ the length measured also in cpu.
By introducing Eqs. (2)-(3) into Eqs. (1a, b), we obtain the so-called "Lane-Emden equation" (cf. [10], Chapter IV, Eq. (11)) which, integrated along a specified integration interval with initial conditions yields as solution the "Lane-Emden function" θ[I ξ ⊂ R] ⊂ R. However, the Lane-Emden equation (5) involves (i) the indeterminate form θ ′ /ξ at the origin, and (ii) the "raised-to-real-power" term θ n which becomes undefined for θ(ξ) < 0 and, also, suffers from a "non-monodromy syndrome" in the sense that multiple-valued logarithmic functions are involved in the representation of θ n (see e.g. [11], Secs. 26-28). To avoid such syndroms, we apply to the IVP established on the equations (5) and (7) the complex-plane strategy. In particular, we assume that the independent variable ξ is a "complex distance", ξ =ξ + iξ ∈ C, and the integration proceeds along a complex path parallel to the real axis and at a relatively small imaginary distance from it. Hence, we perform numerical integration along a contour C ⊂ C, being parallel to the real axis R and distancing iξ 0 from it, i.e. along the straight line-segment joining the points ξ 0 and ξ end in C. For the constant imaginary partξ 0 of the complex distance ξ, we usually take it to lie in the interval 10 −9 , 10 −3 . Thus the Lane-Emden function θ becomes complex-valued function in one complex variable, In the franework of CPS, the initial conditions (7) are written as whereθ 0 = 1 ( cf. Eq. (7a)). The initial value for the imaginary partθ 0 is selected to be small compared to the initial value for the real part, lying usually in the interval [10 −9 , 10 −3 ]. In certain cases, this initial value can be set equal to zero, though systematic numerical experiments show that the presence of a nonzero initial value stabilizes and accelerates the complex-plane integration procedure.
Readers interested on issues of this section can find full details in [9] (Sec. 3.2).

The Global Polytropic Model: Application to the Solar and Jovian Systems
In the so-called "global polytropic model" for the solar system ( [12], Sec. 1), the primary assumption is hydrostatic equilibrium (Eq. (1)). Due to the fact that θ[C ⊂ C] ⊂ C (Eq. (9)), the real partθ has a first root at ξ 1 =ξ 1 + iξ 0 , as well as further roots: a second root at ξ 2 =ξ 2 + iξ 0 , withξ 2 >ξ 1 , a third root at ξ 3 =ξ 3 + iξ 0 , withξ 3 >ξ 2 , etc. Thus the polytropic sphere of polytropic index n and radiusξ 1 is the central component or central body S 1 of a "resultant polytropic configuration" of which further components are the polytropic spherical shells S 2 , S 3 , S 4 , . . . , defined by the pairs of radii (ξ 1 ,ξ 2 ), (ξ 2 ,ξ 3 ), (ξ 3 ,ξ 4 ), . . . , respectively. Each polytropic shell can be considered as an appropriate place for a planet, or a satellite, to be born and live. We speak for a planet when the central body S 1 simulates the Sun [12]; in this case, the resultant polytropic configuration represents the solar system. On the other hand, we speak for a satellite when S 1 simulates a planet, say the Jupiter (for the Jupiter's system of satellites see [13]); then the resultant polytropic configutation represents the jovian system.
The most appropriate location for a planet, or for a planet's satellite, to settle inside a polytropic shell S j is the placeΞ j at which |θ| takes its maximum value inside S j , max|θ[S j ]| = |θ(Ξ j + iξ 0 )|. Now, concerning the solar system, the distance A E of the Earth (the most massive inner planet) from the Sun is the distance A J of the Jupiter (the most massive outer planet) from the Sun is and the distance A N of the Neptune (the most massive among the most distant outer planets) from the Sun is The triplet A E , A J , A N seems to be a representative triplet of distances for the solar system. To compute an "optimum polytropic index" n ⊙ for the Sun, we search for a particular value of the polytropic index which generates a sequence of maximum values max|θ[S j ]|, j = 2, 3, . . . , L, with the integer L taken sufficiently large, occuring at distancesΞ j , j = 2, 3, . . . , L among which there being as close to the given distances A E , A J , A N , as possible; in agreement with Eq. (3b), the unit of length for the Sun is given by α ⊙ = R ⊙ /ξ 1 . Note that, in such computations, the astronomical unit remains invariant, equal to 214.9487 R ⊙ (Eq. (11)) irrespective of the particular valueξ 1 . Computations presented in [12] (Sec. 2, Table I) give as optimum polytropic index for the Sun the value n ⊙ = 3.23.
In [13] (Sec. 1) we have developed a more general algorithm, called A[n], which can compute the optimum polytropic index n opt for a star possessing a planetary system, or for a planet possessing a system of satellites. In detail, to compute an optimum polytropic index n opt in order for a triplet of planets, or satellites, with distances A P 1 < A P 2 < A P 3 from the central body, to be accomodated inside the resultant polytropic configuration, we work as follows. A[n]-1. For a sequence of values n i , i = 1, 2, . . . , N n , we compute the corresponding sequence of distances α pj (n i ) =Ξ j (n i ), j = 2, 3, . . . , L, at which planets/satellites can be accomodated, with the integer L taken sufficiently large. A[n]-2. For each sequence {α pj (n i )}, we compute the two-dimensional "array of distance ratios" A[n]-3. We scan the N n arrays D(n i ; j, k) in order to find a value n opt generating a "maximum number of proper levels" related to the given ratios R 12 = A P 1 /A P 2 and R 13 = A P 1 /A P 3 . Note that two particular elements D(n i ; q, s) and D(n i ; q, t), t > s, constitute a "proper level" (and thus increase by one the number of the proper levels in favor of the polytropic index n i ) if and 100 × |R 13 − D(n i ; q, t)| R 13 ≤ τ, verified within a percent tolerance τ specified by the user. Applying A[n] to the Jupiter's system of satellites ( [13], Sec. 2, Tables I-III), we have computed an optimal value n J = 2.45 for this planet.

The Computations
To compile our programs, we use the gfortran compiler, licensed under the GNU General Public License (GPL; http://www.gnu.org/licenses/gpl.html). gfortran is the name of the GNU Fortran compiler belonging to the GNU Compiler Collection (GCC; http://gcc.gnu.org/). In our computer, it has been installed by the TDM-GCC "Compiler Suite for Windows" (http://tdm-gcc.tdragon.net/), which is free software distributed under the terms of the GPL.
Subroutines required for standard numerical procedures (e.g. interpolations of functions, rootfinding of algebraic equations, localizing extrema of functions, etc.) are taken from the SLATEC Common Mathematical Library, which is an extensive public-domain Fortran Source Code Library, incorporating several public-domain packages. The full SLATEC release is available in http://netlib.cs.utk. edu/.
To solve the complex IVPs involved in this investigation, we use the code DCRKF54 included in the Fortran package dcrkf54.f95 [2]. DCRKF54 is a Runge-Kutta-Fehlberg code of fourth and fifth order modified for the purpose of solving complex IVPs, which are allowed to have high complexity in the definition of their ODEs, along contours (not necessarily simple and/or closed) prescribed as continuous chains of straight-line segments; interested readers can find full details on dcrkf54.f95 in [2].
The header of DCRKF54 is given in [2] (Sec. 2.1, Part #[000]). On entry to DCRKF54, the input arguments are assigned the values A = r start , B = r end , N = n, Y = y start ; DEQS is the subroutine computing the vector derivative(s) function f ([2], Sec. 2.1, Eq. (1a), and discussion preceding this equation), HIN an initial stepsize, HMIN a minimum stepsize, HMAX a maximum stepsize. In this work, the input parameters are assigned the values H = 10 −3 , HMIN = 10 −6 , HMAX = 10 −1 . Furthermore, the input values next to HMAX (discussed in [2], Appendix A) are assigned the values ATOL = 10 −24 , RTOL = 10 −14 for double precision (KD=8, where KD is the overall "kind type parameter" explained in [2], Sec 2.2, discussion following Part #[060]), RTOL = 10 −15 for high precision (KD=10), and QLBD = 7.5 × 10 −1 . In fact, we almost use a pure relative error control, since ATOL is ∼ 10 orders of magnitude less than RTOL. Concerning the equality tolerance XTOL ([2], Sec. 2.2, Part #[050]), it takes the value XTOL = 32 × EMR, where EMR is the well-known "machine roundoff error" (i.e. the larger real number which does not change unity when added to it). Alternatively, the user can add XTOL to the call sequence by modifying the header of DCRKF54 and specify its input value in the calling program. On exit to DCRKF54, A has been hopefully advanced to B, Y is the solution vector at B, and H is the stepsize adapted so far, situation verified by the return value NFLAG = 1. If A needs further steps to arrive at B, then the return value is NFLAG = 2, whence we call again DCRKF54 leaving all of its arguments unchanged.
In this study, integrations proceed along the following members of the con- for satellites of the jovian system, planets of the solar system, and "trans-Neptunian objects" (TNO, TNOs), respectively. The contour class C 2Form represents forward-and-then-backward straight-line routes parallel and close to R obeying the special form (8). The endpoint of any such contour coincides with its start point, ξ end = ξ 0 . Hence, the true valueθ end at the endpoint ξ end coincides with its initial valueθ 0 (Eq. (10)) at the start point ξ 0 . From the numerical analysis point of view, this is an important fact, since the global percentage errors %E(θ) owing to DCRKF54 can be readily calculated. A discussion on various contours and their characteristics can be found in [2] (Sec. 5).

Numerical Results and Discussion
The following two subsections contain numerical results for the purpose of testing the code DCRKF54 and of comparing with previous (published) corresponding results. The third subsection contains results regarding some trans-Neptunian objects. Some earlier (unpublished) corresponding numerical results will not be quoted here, since the computations of the present study are more accurate and reliable.

Satellites of the jovian system
The jovian system of satellites constitutes a short-distance integration problem, since the related complex IVP is solved along the contour C 2B (Eq. (17)). Integrating by the code DCRKF54 has been proved very accurate, since the global percentage error has been found to be %E(θ) ≤ 2 × 10 −9 ; whence, the numerical results in Table 1 can be safely quoted (as they do) with seven decimal digits. As said in Sec. 3, the optimum polytropic index for the Jupiter is n opt (J) = n J = 2.45. All symbols involved in Table 1 are explained in [13] (see especially Sec. 2; comparisons can be made with respective results of Table III).
In both studies, Europa, Ganymede, and Callisto occupy Shells No 4, No 5, and No 7, respectively.

Planets of the solar system
Treating planets of the solar system is a long-distance integration problem, since the so-defined complex IVP is solved along the contour C 2C (Eq. (18)). Numerical integration by the code DCRKF54 gives very accurate results, since %E(θ) ≤ 2×10 −9 ; accordingly, the numerical results in Table 2 are quoted again with seven decimal figures. As said in Sec. 3, the optimum polytropic index for the Sun is n ⊙ = 3.23. All symbols involved in Table 2 are explained in [13] (especially in Sec. 2; comparisons can be made with respective results of Tables  I and II). In both studies, Earth, Jupiter, and Neptune occupy Shells No 7, No 11, and No 18, respectively.

Trans-Neptunian Objects
Computing quantities related to TNOs constitutes a very-long-distance integration problem, since the corresponding complex IVP is solved along the contour C 2D (Eq. (19)). A global percentage error %E(θ) ≤ 9 × 10 −8 has been verified for the code DCRKF54, which is quite satisfactory for this case. Our model reproduces the sharp division between "plutinos" and "classical Kuiper belt objects" (also called "cubewanos"); in particular, the border between Shell No 19 and Shell No 20 is at ∼ 40 AU, and the "Kuiper cliff" appears at ∼ 47 AU.
It seems that Quaoar and Varuna, due to the very small eccentricities of their orbits, are ideal candidates for obeying the global polytropic model (for similar comments regarding jovian satellites, see [13], Sec. 2; for comments regarding planets, see [15], Sec. 2)). In fact, an object in almost circular orbit implies that it has evolved under mild processes, which, in turn, maintain the sensitive global polytropic character of such a system.