Revealing the Newton-Raphson basins of convergence in the circular pseudo-Newtonian Sitnikov problem
Euaggelos E. Zotos Md Sanam Suraj Mamta Jain Rajiv Aggarwal






eez]Euaggelos E. Zotoscor1
evzotos@physics.auth.gr

mss]Md Sanam Suraj

mss]Mamta Jain

ra]Rajiv Aggarwal

[cor1]Corresponding author

[eez]Department of Physics, School of Science,
Aristotle University of Thessaloniki, GR-541 24, Thessaloniki, Greece

[mss]Department of Mathematics, Sri Aurobindo College,
University of Delhi, 110017, Delhi, India

[ra]Department of Mathematics, Deshbandhu College,
University of Delhi, Kalkaji, New Delhi-110019, Delhi, India


In this paper we numerically explore the convergence properties of the pseudo-Newtonian circular restricted problem of three and four primary bodies. The classical Newton-Raphson iterative scheme is used for revealing the basins of convergence and their respective fractal basin boundaries on the complex plane. A thorough and systematic analysis is conducted in an attempt to determine the influence of the transition parameter on the convergence properties of the system. Additionally, the roots (numerical attractors) of the system and the basin entropy of the convergence diagrams are monitored as a function of the transition parameter, thus allowing us to extract useful conclusions. The probability distributions, as well as the distributions of the required number of iterations are also correlated with the corresponding basins of convergence.



Sitnikov problem – pseudo-Newtonian dynamics – Basins of convergence – Fractal basin boundaries






§ INTRODUCTION



The Sitnikov problem is a special formulation of the restricted three-body problem which describes the motion of a test particle which oscillates along the vertical z-axis, perpendicular to the configuration (x,y) plane, where two primaries, of equal masses, move in circular or elliptic orbits with common barycentre i.e. the axes origin. In particular, when the primaries move in circular orbits, the problem describes the so called MacMillan problem <cit.>. The history of the Sitnikov problem, when the primaries move in circular orbit, begins with Ref. <cit.>, where this dynamical model was discussed for the first time. Later on, the Sitnikov problem, of two primary bodies, has been investigated by many authors (see e.g. Refs. <cit.>).

Several types of perturbations, such as primaries with either prolate <cit.> or oblate shape <cit.>, as well as the radiation pressure <cit.>, have been added for making the system of three bodies more realistic. Another well-studied aspect of the Sitnikov three-body problem is the study of the families of periodic orbits and the corresponding bifurcations <cit.>. In addition, the stability of motion in the same system has been investigated in <cit.>, where it was found that in the case where the mass of the test particle is not negligible the energetically allowed regions of motion grow with increasing value of the third body, while at the same time the amplitude of the oscillation, along the vertical z-axis, gradually increases.

As the extension of the restricted three-body problem to the restricted four-body is natural, in the same vein the extension of the Sitnikov three-body problem to the Sitnikov four-body problem is quite obvious. In Ref. <cit.> the existence, the linear stability as well as the bifurcations in the Sitnikov family of straight line periodic orbits in the restricted four-body problem were studied. In this work, it was revealed that the Sitnikov family has only one stability interval, while only twelve 3-dimensional families of symmetric periodic orbits exist, which bifurcate from twelve corresponding critical Sitnikov periodic orbits. On the other hand, the Sitnikov family has infinitely many stability intervals in the restricted three-body Sitnikov problem, which results to infinitely many Sitnikov critical orbits.

Over the years, the Sitnikov four-body problem has been investigated by many authors (see e.g., Ref. <cit.>), including also the oblateness of the primaries (see e.g., Refs. <cit.>), the effects of the radiation pressure (see e.g., Ref. <cit.>), and families of periodic orbits and bifurcations (see e.g., Ref. <cit.>). Moreover, the region of motion in the Sitnikov four-body problem when the fourth mass is finite has also been determined <cit.>, by exploring the structure of the phase space with the use of proper selection of surfaces of section. It was observed that for low energy, the central fixed point is stable, while for higher values of the energy it splits into an unstable and two stable fixed points.

In Ref. <cit.> the stability of the vertical motion and its bifurcations into families of 3-dimensional periodic orbits in the Sitnikov restricted N-body problem were revealed. Additionally, it was found that for N ≥ 4 there exists only one interval of stable vertical solutions for every N. However, as N increases, the stability intervals increase in size. They also demonstrated that a sizable region of bounded orbits close to the z-axis also exists. More precisely, these orbits are small near the endpoints of the interval, while they increase in size at the middle of the same interval.

For modelling the motion of the infinitesimal mass (test particle) in the classical N-body problem, especially for the three and four body-problems, various, more realistic, modifications have been proposed, mainly by introducing additional perturbative terms to the classical Newtonian effective potential. By using the Einstein-Infeld-Hoffmann theory <cit.> the first-order post-Newtonian equations of motion, in the restricted three-body problem, have been deduced by several authors (see e.g., Refs. <cit.>). Similarly, the orbital dynamics of the planar circular restricted three-body problem, in the context of the pseudo-Newtonian approximation by using the Fodor-Hoenselaers-Perjé procedure, has been analyzed in Ref. <cit.>. Furthermore, in Ref. <cit.> the influence of the separation between the primary bodies was discussed and it was revealed that the post-Newtonian dynamics substantially differ from the corresponding classical Newtonian dynamics, provided that the distances between the primary bodies is sufficiently small.

Undoubtedly, in a dynamical system an issue of great importance is knowing the locations of the equilibrium points. However, in most of the cases, such as in the restricted N-body problem (with N ≥ 3), the position of the libration points cannot be computed by means of analytical methods. Consequently, the only viable alternative is the use of numerical methods. Unfortunately, the outcomes of any numerical technique are directly related to the initial conditions, used to the iterative procedure. More precisely, for initial conditions inside the basins of convergence the iterative procedure lead very quickly to a solution. On the other hand, for initial conditions located in the vicinity of the fractal basin boundaries the numerical scheme usually requires a considerable amount of iterations for reaching to a libration point (root). Therefore, the convergence properties of a system is a highly important issue, because this information reveals the optimal (regarding fast convergence) initial conditions, for numerically locating an equilibrium point.

Over the past decades, the basins of convergence, associated with the libration points, have been investigated mainly by using the multivariate version of the Newton-Raphson method. At the same time, various types of perturbing terms have been added the effective potential of the restricted problem of three and four bodies (see e.g., Refs. <cit.>). Very recently, in Ref. <cit.> we discussed the basins of convergence, associated with the equilibrium points in the pseudo-Newtonian planar circular restricted three-body problem, by using the multivariate version of the Newton-Raphson iterative scheme. Our analysis revealed that the transition parameter strongly influences the existence and the stability of the libration points as well as the topology of the basins of convergence.

For almost all the N-body systems (with N ≥ 3) we have no analytic equations, regarding the position of the equilibrium points of the systems. This directly implies that we have to use numerical methods for determining the position of the libration points. However, the outcomes of the numerical algorithms are affected by the specific initial conditions, which we use as starting points of the iterative schemes. More precisely, there are initial conditions which display a very fast convergence, while there are also starting points for which the iterative procedure requires a substantial number of iterations before reaching to a solution. The first type of initial conditions compose the so-called “basins of convergence", while the latter type usually appears in the vicinity of the fractal basin boundaries. Therefore, it is undeniable that knowing the convergence regions of a dynamical system is an issue of paramount importance, because they indicate the optimal initial conditions for the iterative schemes. At the same time, the basins of convergence show us all the fractal regions in which all the corresponding initial conditions should be avoided as starting points.

In Ref. <cit.> it was proved that some pseudo-Newtonian systems are more dynamically stable, compared to their Newtonian equivalent. On this basis, it is very interesting to explore the pseudo-Newtonian version of the Sitnikov problem of three and four primary bodies. The main aim of this work is to unveil the convergence properties of the system by determining how the transition parameter influences the geometry as well as the shape of the Newton-Raphson basins of convergence.

It should be emphasized that the configuration of the circular restricted four-body problem, with equally massed primaries, is dynamically unstable (the corresponding proof is given in <ref>). Nevertheless, we decided to conduct numerical calculations in this system, mainly for comparison reasons of the applied numerical method's behavior to the computation of the equilibrium points, as well as the associated basins of convergence. Additionally, several previous works have been devoted on the four-body problem with three equal masses (see e.g., Refs. <cit.>).

The present paper has the following structure: the most important properties of the dynamical system are presented in Section <ref>. The parametric evolution of the position of the equilibrium points is investigated in Section <ref>. The following Section contains the main numerical results, regarding the evolution of the Newton-Raphson basins of convergence. In Section <ref> we demonstrate how the transition parameter influences the basin entropy. Our paper ends with Section <ref>, where we emphasize the main conclusions of this work.




§ PROPERTIES OF THE DYNAMICAL SYSTEM



Let us briefly recall the main properties of the dynamical systems under consideration. A dimensionless, barycentric rotating system of coordinates Oxyz, is considered for both cases. For the description of the planar motion of the test particle we choose a rotating reference frame, where the center of mass of the primaries coincides with its origin.


  
  * The pseudo-Newtonian circular restricted three-body problem: Two primary bodies, P_1 and P_2, are located on the Ox axis with masses m_1 = μ and m_2 = 1 - μ, respectively, where μ = m_2/(m_1 + m_2) ≤ 1/2 is the mass parameter <cit.>. The centers of both primaries are located at (x_1, 0, 0) and (x_2, 0, 0), where x_1 = - μ and x_2 = 1 - μ.

      According to Ref. <cit.> the effective potential function of the pseudo-Newtonian circular restricted three-body problem, with only the first correction terms, is given by
      
    Ω_3(x,y,z) = ∑_i=1^2m_i/r_3i - ϵ/2 c^4∑_i=1^2m_i^3/r_3i^3 + 1/2(x^2 + y^2),

      where of course
      
    r_31   = √((x - x_1)^2 + y^2 + z^2),
    
          r_32   = √((x - x_2)^2 + y^2 + z^2),

      are the distances of the test particle from the two primary bodies.
  
  * The pseudo-Newtonian circular restricted four-body problem: Three primary bodies, P_1, P_2, and P_3, with masses m_1, m_2, and m_3, respectively are located on the configuration (x,y) plane, while their mutual distances remain constant and form an equilateral triangle. In the special case, where the three primaries have equal masses (m_1 = m_2 = m_3 = 1/3) the coordinates of the centers are (x_1,0,0), (x_2,y_2,0), and (x_3,y_3,0), where x_1 = √(3)/3, x_2 = - √(3)/6, x_3 = x_2, y_2 = 1/2, and y_3 = - y_2.

      In the same vein of the pseudo-Newtonian three-body problem, the effective potential function of the pseudo-Newtonian four-body problem (where again we consider only the first correction terms) is
      
    Ω_4(x,y,z) = ∑_i=1^3m_i/r_4i - ϵ/2 c^4∑_i=1^3m_i^3/r_4i^3 + 1/2(x^2 + y^2),

      where once more
      
    r_41   = √((x - x_1)^2 + y^2 + z^2),
    
          r_42   = √((x - x_2)^2 + (y - y_2)^2 + z^2),
    
          r_43   = √((x - x_3)^2 + (y - y_3)^2 + z^2),

      are the distances of the third body from the three primaries. In <ref> we prove that the equilateral triangular configuration of the system of three primary bodies is always unstable.
      

In both cases, the parameter ϵ acts as a transition parameter from classical Newtonian to pseudo-Newtonian dynamics, while its values lie in the interval ϵ∈ [0, 1]. Indeed, when ϵ = 0 both effective potential functions, Ω_3(x,y,z) and Ω_4(x,y,z), are reduced to the corresponding versions of classical Newtonian dynamics. On the other hand, when ϵ = 1 we have the case of the full pseudo-Newtonian problem.

We adopt a system of units in which the sum of the masses of the primaries, the angular velocity (ω), the speed of light c, the distance R between the primaries as well as the gravitational constant G are equal to unity.

The equations of motion describing a test particle (a body of a negligible mass m, with respect to the masses of the primaries) moving under the mutual gravitational attraction of the primaries read

    ẍ - 2 ẏ = ∂Ω_i/∂ x, 
       ÿ + 2 ẋ = ∂Ω_i/∂ y, 
       z̈ = ∂Ω_i/∂ z,

where i = 3,4.

The above system of differential equations admits only one integral of motion (also known as the Jacobi integral), which is
given by the following Hamiltonian

    J(x,y,z,ẋ,ẏ,ż) = 2Ω_i(x,y,z) - (ẋ^2 + ẏ^2 + ż^2 ) = C,

where ẋ, ẏ, and ż are the velocities, while C is the numerical value of the Jacobi constant which is conserved.



The potential function of the circular Sitnikov problem can be obtained if we set equal masses to all primaries and x = y = 0. For the restricted three-body problem with μ = 1/2 we have

    Ω_3(z) = 1/r_3(1 - ϵ/8r_3^2),

where r_3 = √(z^2 + 1/4). Similarly, the potential function of the circular Sitnikov pseudo-Newtonian four-body problem reads

    Ω_4(z) = 1/r_4(1 - ϵ/18r_4^2),

where r_4 = √(z^2 + 1/3). It is evident that Eqs. (<ref>) and (<ref>) describe the motion of a massless test particle which oscillates along a straight line which is perpendicular to the orbital configuration (x,y) plane of the primary bodies with equal masses. In Fig. <ref> we present the geometry of the Sitnikov problem, related to the pseudo-Newtonian circular three and four-body problems.

Consequently, the equations, regarding the motion of the test particle along the vertical z axis, have the form

    z̈   = d/d z(Ω_3(z)) = - z/r_3^3(1 - 3ϵ/8r_3^2), 
    z̈   = d/d z(Ω_4(z)) = - z/r_4^3(1 - ϵ/6r_4^2),

while the corresponding Jacobi integral, for the vertical motion, becomes

    J_z(z,ż) = 2 Ω_i(z) - ż^2 = C_z,

with i=3,4.




§ PARAMETRIC VARIATION OF THE EQUILIBRIUM POINTS



From now on, the z coordinate is considered as a complex variable and it is denoted by , thus following the approach which was successfully used in Refs. <cit.>. At this point, it should be explained that the transition to complex variables is imperative. This is because, as it was presented in Ref. <cit.>, all the impressive and beautiful fractal basin structures can be observed only on the complex plane.

In order to obtain the equilibrium points (roots) of the vertical motion, the right hand side of Eq. (<ref>) should be taken equal to zero as

    f_3(;ϵ)    = - 4(8^2 - 3ϵ + 2)/(1 + 4^2)^5/2 = 0, 
    
    f_4(;ϵ)    = - 3 √(3)(6^2 - ϵ + 2)/2(1 + 3^2)^5/2 = 0,

which are reduced to

    (8^2 - 3ϵ + 2) = 0, 
       (6^2 - ϵ + 2) = 0.




Looking at Eqs. (<ref>) we observe that the root = 0 is always present, regardless the value ϵ of the transition parameter. This root corresponds to the inner collinear equilibrium point L_1 of the circular restricted three and four-body problems. However since the left hand side of Eqs. (<ref>) is a third order polynomial it means that there are two additional roots, given by

    _i = ±√(3ϵ - 2)/2√(2),    i = 1,2,

in the case of two primary bodies and

    _i = ±√(ϵ - 2)/√(6),    i = 1,2,

when three primaries are present.

The nature of these two roots strongly depends on the numerical value ϵ of the transition parameter. Our numerical analysis reveals that for the circular Sitnikov pseudo-Newtonian three-body problem, along with the = 0 root

  
  * When ϵ < -2/3 there are two imaginary roots.
  
  * When ϵ = 0 or ϵ = 2/3 only the root = 0 exists.
  
  * When ϵ > 2/3 there are two real roots.

while for the circular Sitnikov pseudo-Newtonian four-body problem applies that

  
  * When ϵ < -2 there are two imaginary roots.
  
  * When ϵ = 0 or ϵ = 2 only the root = 0 exists.
  
  * When ϵ > 2 there are two real roots.

However since the numerical value of the transition parameter ϵ cannot exceed 1, real roots are impossible for the circular Sitnikov pseudo-Newtonian four-body problem. It should be noted that ϵ = 1 is the maximum allowed value of the transition parameter, corresponding to full pseudo-Newtonian dynamics, while all higher values (ϵ > 1) have no physical meaning. It is seen, that the value ϵ = 2/3 is in fact a critical value of the transition parameter, since it determines the change on the nature of the two roots, in the case of the circular Sitnikov pseudo-Newtonian three-body problem.



Useful information could be extracted from the parametric variation of the roots on the complex plane, as a function of the transition parameter. Fig. <ref>(a-b) shows the parametric evolution of the roots on the complex plane, when ϵ∈ [0,1], with ℛ = Re[] and ℐ = Im[]. For the case where two primaries are present we see in panel (a) of Fig. <ref> that as long as ϵ > 0 two imaginary roots R_12 and R_13 appear. As the value of the transition parameter increases both imaginary roots tend to the origin and when ϵ = 2/3 they collide with the central root R_1 and only the root = 0 survives thus having multiplicity 3. When ϵ > 2/3 two real roots R_22 and R_23 emerge from the origin O and they move away from each other, with increasing value of ϵ. Finally, when ϵ = 1 we the two real roots ±√(2)/4 reach their maximum distance from the origin. When three primaries are present it is seen in panel (b) of Fig. <ref> that the evolution of the imaginary roots R_2 and R_3 follow a similar evolution with respect to the first case (with roots R_12 and R_13) of two primaries. When ϵ = 1 we have the imaginary roots ± i √(6)/6 and their evolution is terminated, thus avoiding the collision with the origin.




§ THE NEWTON-RAPHSON BASINS OF CONVERGENCE











For obtaining the basins of convergence on the complex plane, associated with the roots of the system, we have to numerically solve the corresponding mono-parametric equation. The easiest way is to use the classical Newton-Raphson method of second order, through the iterative scheme

    _n+1 = _n - f_i(_n;ϵ)/f_i'(_n;ϵ),    i=3,4,

where _n is the value of the  at the n-th step of the iterative process. By plugging formulae (<ref>) in Eq. (<ref>) we obtain

    _n+1 = 12^3 (8^2 - 5 ϵ + 2 )/8^2 (8^2 - 6 ϵ + 1) + 3ϵ - 2,

for the case of two primary bodies and

    _n+1 = 3^3 (18^2 - 5 ϵ + 6 )/6^2 (6^2 - 2 ϵ + 1) + ϵ - 2,

when three primaries exist on the configuration plane.

The numerical algorithm for revealing the basins of convergence on the complex plane works as follows: The iterative procedure is activated by an initial condition of the form of a complex number = a + ib, with ℛ = a and ℐ = b and continues until a root is reached, with the predefined numerical accuracy. In our case the numerical accuracy is set to 10^-15, for both real and imaginary parts. Using a dense uniform grid of 1024 × 1024 (ℛ,ℐ) nodes, as initial conditions, we perform a double scan of the complex plane. During the iterative procedure we monitor the required number N of iterations, along with the classification of the nodes. For our experiments the maximum allowed number of iterations is set to N_ max = 500.

If the iterative procedure leads to a root it means that the Newton-Raphson method converges for the particular initial condition. Here it should be clarified that the numerical method, in general terms, does not converge equally well for all the initial conditions on the complex plane. All the initial conditions which lead to the same root (final state of the iterative procedure) compose the Newton-Raphson basin of convergence/attraction (we will also use the terms attractive regions/domains). Furthermore, we would like to emphasize that the Newton-Raphson basins of convergence should not be mistaken with the basins of attraction which are encountered in dissipative systems.

The Newton-Raphson basins of convergence when ϵ = 0 are presented in panels (a), for the case of two primaries, and (c), for the case of three primaries, of Fig. <ref>. We see that the converging initial conditions are mainly located near the center, while they form a rhomboidal fractal shape. On the other hand, the vast majority of the complex plane is covered by a unified sea (yellow region) of initial conditions for which the Newton-Raphson iterative scheme tends very quickly to extremely large complex numbers. This numerical behaviour implies that for these initial conditions tend, theoretically, to infinity. In both cases (two and three primaries) the overall geometry of the converging regions are qualitatively very similar. This can be explained, in a way, if we examine the form of the corresponding iterative schemes. For the case of two primary bodies the Newton-Raphson iterative scheme with ϵ = 0 takes the form

    _n+1 = 12^3/8^2 - 1,

while for the case of three primaries we have that

    _n+1 = 9^3/6^2 - 1.

We observe that the complex functions entering both iterative schemes are completely identical and only the numerical coefficients are different. This should be the main reason of the similar patterns observed in panels (a) and (c) of Fig. <ref>. In panels (b) and (d) of the same figure the distribution of the corresponding number (N) of iterations required for obtaining the desired accuracy is given using tones of blue.

In the following subsections we will determine how the transition parameter ϵ affects the structure of the Newton-Raphson basins of convergence in the circular Sitnikov problem with two and three primary bodies. For the classification of the nodes on the complex plane we will use color-coded basin diagrams, in which each pixel is assigned a different color, according to the final state (root) of the corresponding initial condition.



 §.§ Case I: Two primary bodies present and 0 < ϵ < 2/3



We begin with the first case where two primary bodies are present on the configuration (x,y) plane, while the equation f_3(;ϵ) = 0 has, apart from the = 0 root, two imaginary roots. The Newton-Raphson basins of convergence on the complex plane, for three values of the transition parameter, are illustrated in the first column of Fig. <ref>. It is seen that in all cases the area of all the types of the basins of convergence is finite. On the contrary, outside the convergence regions the vast majority of the complex plane is covered by initial conditions which do not converge to any of the three roots (yellow regions) but they tend asymptotically to infinity.

In the second column of Fig. <ref> we present the corresponding number N of iterations, using tones of blue, while the corresponding probability distribution of the required iterations is given in the third column of the same figure. The definition of the probability P is the following: if N_0 complex initial conditions (ℛ,ℐ) converge, after N iterations, to one of the roots then P = N_0/N_t, where N_t is the total number of nodes in every basin diagram. Moreover, in all plots the tails of the histograms extend so as to cover 97% of the corresponding distributions of iterations. The vertical, red, dashed lines in the probability histograms denote the most probable number N^* of iterations. The blue lines in the histograms of Fig. <ref> indicate the best fit to the right-hand side N > N^* of them (more details are given in subsection <ref>).

We see that as soon as ϵ > 0 several lobes emerge at the boundaries of the convergence region corresponding to the central root R_1. With increasing value of the transition parameter the most important changes which occur on the complex plane are the following:

  
  * The area of the basins of convergence corresponding to all three roots increase, while at the same time the overall shape of the convergence region changes from rhomboidal to spherical.
  
  * The probability distribution of the required iterations moves to the right, which implies that additional iterations are required. Indeed, we observe that the most probable number N^* of iterations increases from 6, when ϵ = 0.2 to 13 when ϵ = 0.666.
  
  * In panel (h), where ϵ = 0.666 it is seen that throughout the convergence regions the corresponding iterations are significantly higher with respect to panel (b), where ϵ = 0.2. We suspect that this behaviour is directly related with the fact that we approaching the critical value ϵ = 2/3.




 §.§ Case II: Two primary bodies present and 2/3 ≤ϵ≤ 1



The next case under consideration involves the scenario where there are two real roots, along with the = 0 root. In the first column of Fig. <ref> we present the Newton-Raphson basins of convergence for three values of the transition parameter. In panel (a) of Fig. <ref>, which corresponds to the critical value of ϵ one can see that the convergence region corresponding to root R_1 is a perfect circle with radius 0.5. However, in panel (b) of the same figure one can observe that the initial conditions which form the circular region produce interesting inner structures, depending on the number of the required iterations.

As we proceed to higher values of the transition parameter the main changes, regarding the geometry of the convergence areas, are the following:

  
  * The extent of the basins of convergence of all three roots increases, while the geometry of the overall converging region alters from fractal spherical to elliptic (oval). As long as ϵ > 2/3 two main lobes appear, corresponding to basins of convergence of the central roots, while these are lobes gradually unified thus producing a singe region.
  
  * The probability distribution of the required iterations moves to the left, which implies that less iterations are required. Indeed, it is observed that the most probable number N^* of iterations decreases from 50, when ϵ = 2/3 to 6 when ϵ = 1.
  
  * The basin boundaries become more noisy, which means that the fractality of the several basin boundaries on the complex plane increases.




 §.§ Case III: Three primary bodies present



We close with the last case where three primary bodies circulate on the configuration (x,y) plane, while the equation f_4(;ϵ) = 0 has two imaginary roots, along with the R_1 root = 0. The evolution of the geometry of the Newton-Raphson basins of convergence is depicted in the first column of Fig. <ref>, where we present three basin diagrams for three values of the transition parameter. It becomes evident that this is the least interesting case, from the dynamical point of view.

Indeed, with increasing value of ϵ the overall structure of the complex plane hardly changes and the most important phenomena which take place are the following:

  
  * The area of the basin of convergence, corresponding to the central root R_1, increases, while at the same time a smaller increase on the extent of the other types of converging regions (corresponding to roots R_2 and R_3) is observed.
  
  * Apart from the central lobes, corresponding to roots R_2 and R_3, several smaller convergence regions emerge at the boundaries of the main basin of convergence, corresponding to root R_1.
  
  * The most probable number N^* of iterations displays a minor reduce from 6, when ϵ = 0.4, to 5, when ϵ = 1.




 §.§ An overview analysis



Since the basin diagrams, presented earlier in subsections <ref>, <ref>, and <ref>, have a fixed and equal size we could monitor the evolution of the percentages of the different types of initial conditions, as a function of the transition parameter ϵ. Such diagrams are presented in Fig. <ref>(a-b). It is seen that in both cases, regarding the number (two or three) of the primary bodies, the percentage of initial conditions for which the Newton-Raphson iterative scheme tends to infinity is constantly reduced. On the other hand, the rates of all the types of converging initial conditions are increased, following an almost linear growth. In Panel (a) of Fig. <ref> we see that the evolution of all types of the initial conditions is very smooth, except near the critical value ϵ = 2/3, where the dynamical properties of the system change.

Additional interesting information could be extracted from the probability distributions of iterations presented in the third row of the basin diagrams. More precisely, we could try to obtain the best fit for the tails of the probability histograms. Using the term “tails" we refer to the right-hand side of the histograms, where N > N^*, where of course N^* is the most probable value. The ideal choice for fitting the tails of the histograms is the Laplace distribution, due to the fact that this particular type of distribution is commonly used in systems where transient chaos is present (see e.g., Refs. <cit.>). Our analysis indicates that the Laplace distribution is, in the vast majority of the cases, the best fit for the corresponding data.

For the Laplace distribution the corresponding probability density function (PDF) is defined as

    P(N | a,b) = 1/2bexp(- a - N/b),    if  N < a 
    exp(- N - a/b),    if  N ≥ a
     .

It is seen that two quantities enter the function P(N). The first one a is widely known as the location parameter, while b is of course the diversity. Note that in our case, only the N ≥ a part of the distribution function is needed for fitting the tails of the histograms.

The aggregated results regarding the numerical values of the location parameter as well as the diversity, for all the studied cases, are presented in Table <ref>. It is very interesting to note that the most probable number of iterations is, in most of the cases, very close to the location parameter a, while in some of them both these quantities completely coincide, thus indicating the ideal choice of the Laplace distribution.






§ PARAMETRIC EVOLUTION OF THE BASIN ENTROPY







In the basin diagrams of the previous section we observed the presence of highly fractal regions, mainly located near the vicinity of the basin boundaries. By using the term fractal we simply imply that the particular areas display a fractal-like geometry, however without computing the corresponding fractal dimension as in Refs. <cit.>. In Fig. <ref> we present a magnification of a local area corresponding to the basin diagram shown earlier in panel (g) of Fig. <ref>, where ϵ = 0.666. Here the fractal basin boundaries are much more clear and distinct. It is known that the final state (root) of initial conditions inside these fractal areas is highly sensitive. Specifically, even the slightest change of the initial conditions automatically leads to a completely different root, which is a classical indication of chaos. Therefore, for the initial conditions located close to the basin boundaries it is almost impossible to predict their final states (roots).

To demonstrate the sensitivity of the fractal regions we chose three initial conditions inside a fractal region of Fig. <ref>. All three initial conditions have the same imaginary part (0.2416), while their real parts differ only at the fourth decimal figure, which implies that the initial conditions are extremely close with each other. In panel (a) of Fig. <ref> we depict the crooked path line created by the successive approximations-points that are followed by the Newton-Raphson iterator. We see that for the first eleven iterations the Newton-Raphson iterator follows almost identical paths for all three initial conditions. However after the eleventh iteration the three paths start to diverge, thus leading to three different final states (roots R_1, R_12, and R_13). Similar behavior applies for all the nodes with initial conditions inside the fractal areas of the basin diagrams.

In the previous section for describing the degree of fractality of the basin diagrams on the complex plane we used only qualitative arguments. In order to enrich our analysis we must also present quantitative information, regarding the fractal geometry of the basins of convergence. For this purpose we decided to use the basin entropy which was recently introduced in <cit.> and measures the degree of the basin fractality (or unpredictability). This new tool examines the topological properties of the convergence regions and provides quantitative results about their fractality.



The algorithm which describes how the basin entropy works is the following: Let us assume that we define a certain region R = [-1,1] × [-1,1] on the complex plane which we later divide into a rectangular grid of N squares boxes (cells). Inside each of these cells there might exist between 1 and N(A) attractors (roots in our case). Then we denote as P_i,j the probability that inside a square box i the attractor is j. Obviously, the initial conditions (nodes), inside each box, are completely independent. Therefore, the Gibbs entropy of each square cell is given by

    S_i = ∑_j=1^m_iP_i,jlog_10(1/P_i,j),

where m_i∈ [1,N_A] is the total number of the attractors (roots) that exist inside each square box i.

Adding all the individual entropies of the N square cells of the rectangular grid on the complex plane we obtain the total entropy of the region R as

    S = ∑_i=1^N S_i.

Consequently, the total entropy, associated with the entire amount of N cells, is called basin entropy and it can be calculated as

    S_b = 1/N∑_i=1^N∑_j=1^m_iP_i,jlog_10(1/P_i,j).


Using the value ε = 0.005, as suggested in Ref. <cit.>, we calculated the numerical value of the basin entropy S_b of the complex plane, for several values of the transition parameter ϵ. At this point, it should be emphasized that the initial conditions, on the complex plane, for which the Newton-Raphson iterative scheme fails to converge, were counted as an additional type of basin, which coexist along with the regular basins of convergence, associated with the equilibrium points of the system. Panel (a) of Fig. <ref> illustrates the parametric evolution of the basins entropy, as a function of the transition parameter. Here it should be noted that for this diagram we used results not only from the cases, of Figs. <ref>, <ref>, and <ref>, but also from additional levels of the transition parameter ϵ.

One may observe that as the value of the transition parameter increases the basin entropy also increases for both cases (two and three primary bodies). Looking the diagram of Fig. <ref>a we can identify two phenomena which should be emphasized: (i) at the critical value of the transition parameter the basin entropy displays a sudden as well as huge drop, almost to zero however, as soon as ϵ > 2/3 it climbs up again and (ii) at the full pseudo-Newtonian case, where ϵ = 1, the basin entropy, corresponding to the Sitnikov problem with two primary bodies, has almost twice the value of the case of the Sitnikov problem with three primaries.

Despite the thorough study of the basins of convergence, their fractal boundaries are only visible in magnifications of local areas of the basin diagrams (see e.g., Fig. <ref>). In order to obtain quantitative results regarding the fractal nature of the basins we computed the boundary basin entropy <cit.>

    S_bb = S/N_b

where N_b is the number of boxes containing more than one attractor.

According to the so-called “log 2 criterion", the basin boundaries are not smooth (fractal) if S_bb > log 2 (the converse is not true). In panel (b) of Fig. <ref> we demonstrate the evolution of the boundary basin entropy S_bb, as a function of the transition parameter ϵ. It is evident that in both cases (systems of two and three primary bodies) it is always S_bb > log 2 which implies that the basin boundaries are always fractal in both systems under consideration.




§ CONCLUDING REMARKS



The convergence properties of the pseudo-Newtonian circular Sitnikov problem of three and four primary bodies have been numerically investigated. Using the classical, yet very effective, Newton-Raphson iterative scheme, we managed to reveal the basins of convergence on the complex plane. Additionally, we successfully determined the influence of the transition parameter ϵ on the roots as well as on several important quantities of the system.

As far as we know, there are no previous related numerical studies on the convergence properties of this particular dynamical system. Therefore, all the contained outcomes are novel and add to our existing knowledge on the basins of convergence in dynamical systems.

The following list contains the most important results of our numerical investigation:

  
  * Imaginary roots are possible for both cases, regarding the number of the primary bodies (two and three). Real roots on the other hand, are possible only in the case of the Sitnikov problem corresponding to the pseudo-Newtonian circular restricted three-body problem.
  
  * It was found that all the basins of convergence, corresponding to all three roots, have finite area, regardless the value of the transition parameter. Our numerical analysis indicates that the vast majority of the complex plane is covered by initial conditions which do not converge to any of the three roots. Furthermore, additional computations revealed that for all these initial conditions the Newton-Raphson iterator lead very fast to extremely large complex numbers (either real or imaginary), which implies that for these initial conditions the numerical method converges to the infinity.
  
  * It should be emphasized that our classification of the initial conditions on the complex plane did not report any false-converging nodes. It should be explained that by the term “false-converging" nodes we refer to initial conditions for which the iterative scheme leads (for N < N_ max) to final states which are not roots of the system, thus displaying a false convergence.
  
  * Near the critical value of the transition parameter we identified several types of converging areas for which the corresponding number of required iterations is relatively high, with respect to near by basins of other roots. We suspect that this phenomenon is inextricably linked with the fact that near the critical point the dynamics of the system, such as the total number of the equilibrium points (roots), drastically changes.
  
  * It was observed that the basin entropy of the complex plane is highly influenced by the transition parameter. More precisely, the highest value of S_b is exhibited when ϵ→ 1 (post-Newtonian dynamics), while on the other hand the basin entropy tends to zero when ϵ→ 0 (classical Newtonian dynamics).


For classifying the nodes on the complex plane we used a double precision  code <cit.>. The required CPU time, per grid of initial conditions, was less then 5 minutes, using an Intel^ Quad-CoreTM i7 2.4 GHz PC. The latest version 11.3 of the software Mathematica^ <cit.> was used for developing all the graphical illustration of the paper.

In this article we used the simplest iterative method (the Newton-Raphson) for reveling the basins of convergence on the complex plane. An undeniably challenging task for a future work would be to use other iterative schemes (of order n > 2) and determine the corresponding similarities and differences on the convergence properties.




§ ACKNOWLEDGMENTS


The authors would like to express their warmest thanks to the two anonymous referees for the careful reading of the manuscript and for all the apt suggestions and comments which allowed us to improve both the quality and the clarity of the paper.






§ STABILITY OF THE EQUILATERAL TRIANGULAR CONFIGURATION OF THE FOUR-BODY PROBLEM



According to <cit.> the stability of the equilateral triangular configuration of the four-body system strongly depends on the masses of the three primaries, as well as on the type of forces between the three main bodies. In particular, the system is stable only if

    (m_1 + m_2 + m_3 )^2/m_1 m_2 + m_1 m_3 + m_2 m_3 > 3 (1 + κ/3 - κ)^2,

where κ is the power of the attraction law between the primaries.

For the case where all three primary bodies have equal masses (that is for m_1 = m_2 = m_3 = 1/3), it can be easily derived, from the above criterion, that the equilateral triangular configuration of the tree primaries is unstable not only in the classical Newtonian problem (where κ = 2) but also in the pseudo-Newtonian problem (where κ = 4).








McM11 W.D. McMillan, Astronomical Journal 27 (1911) 11.

P07 P. Pavanini, Sopra una nuova categoria di soluzioni periodiche nel problema dei tre corpi Ann. Math. SerieIII, Tomo XIII 1907.

DS97 R. Dvorak, Y., Sui Sun, Celest. Mech. Dyn. Astron. 67 (1997) 87.

H09 J. Hagel, Celest. Mech. Dyn. Astron. 103 (2009) 251.

JE01 L. Jiménez-Lara, A, Escalona-Buendía, Celest. Mech. Dyn. Astron. 79 (2001) 97.

JP97 M.A., Jalali, S.H. Pourtakdoust, Celest. Mech. Dyn. Astron. 68 (1997) 151.

Per07 E.A. Perdios, Celest. Mech. Dyn. Astron. 99 (2007) 85.

PM87 E.A. Perdios, V.V. Markellos, Celest. Mech. 42 (1987) 187.

PK12 E.A. Perdios, V.S. Kalantonis, Celest. Mech. Dyn. Astron. 113 (2012) 377.

S65 K. Stumpff, Himmelsmechanik, Band II, VEB Deutscher Verlag der Wissenschaften, Berlin, 1965.

Z17b E.E. Zotos, Chaos, Solitons & Fractals 101 (2017) 8.

DKMP12 C.N. Douskos, V.S. Kalantonis, P. Markellos, E. Perdios, Astrophys. Space Sci. 337 (2012) 99.

RGH15 Rahman, M.A., Garain, D.N., Hassan, M.R., Astrophys. Space Sci. 357 (2015) 64.

PK06 Perdios, E.A., Kalantonis, V. S., In Recent Advances in Astronomy and Astrophysics, 848 (2006) 743.

BLO94 Belbruno, E., Llibre, J., Ollé, M., Celest. Mech. Dyn. Astr. 60 (1994) 99.

KPP08 Kalantonis, V.S., Perdios, E.A., Perdiou, A. E., Astrophys. Space Sci. 315 (2008) 323.

SBV07 Soulis, P.S., Bountis, T., Dvorak, R., Celest. Mech. Dyn. Astron. 99 (2007) 129.

SPB08 P.S. Soulis, K.E. Papadakis, T. Bountis, Celest. Mech. Dyn. Astron. 100 (2008) 251.

SH11 M.S. Suraj, M.R. Hassan, Proceedings of the Pakistan academy of sciences, 48, (2011) 117.

PA13a L.P. Pandey, I. Ahmad, Astrophys. Space Sci. 345 (2013) 73.

SH13 M.S. Suraj, M.R. Hassan, Proceedings of the Pakistan academy of sciences, 50 (2013) 61.

SH14 M.S. Suraj, M.R. Hassan, Astrophys. Space Sci. 349 (2014) 705.

PA13b L.P. Pandey, I. Ahmad, Astrophys. Space Sci. 348 (2013) 377.

BP09 T. Bountis, K.E. Papadakis, Celest. Mech. Dyn. Astron. 104 (2009) 205.

EIH38 A. Einstein, L. Infeld, B., Hoffmann, Annals of Mathematics 39 (1938) 65.

B72 V.A. Brumberg, Relativistic Celestial Mechanics Nauka, Moscow, 1972.

C76 G. Contopoulos, In: In Memoriam D. Eginitis, Kotsakis, D. (eds.) Athens, 1976, p. 159.

K67 E. Krefetz, Astronomical Journal 72 (1967) 471.

DLCG17 F.L. Dubeibe, F.D. Lora-Clavijo, A.G. Guillermo, Physics Let. A 381 (2017) 563.

HW14 G. Huang, X. Wu, Physical Review D 89 (2014) 124034.

SAA17 M.S. Suraj, R. Aggarwal, M. Arora, Astrophys. Space Sci. 362 (2017) 159.

SAP17 M.S. Suraj, M.C., Asique, U. Prasad, M.R. Hassan, K. Shalini, Astrophys. Space Sci. 362 (2017) 211.

Z16 E.E. Zotos, Astrophys. Space Sci. 361 (2016) 181.

Z17a E.E. Zotos, Astrophys. Space Sci. 362 (2017) 195.

ZSAS18 E.E. Zotos, Md S. Suraj, R. Aggarwal, S.K. Satya, Few Body Systems 59 (2018) 69.

Z18b E.E. Zotos, International Journal of Non-Linear Mechanics 103 (2018) 93-103.

SL06 A.F. Steklain, P.S. Letelier, Physics Letters A 352 (2006) 398.

BP11a A.N. Baltagiannis K.E. Papadakis, International Journal of Bifurcation and Chaos 21 (2011) 2179.

BP11b A.N. Baltagiannis K.E. Papadakis, Astrophys. Space Sci. 336 (2011) 357.

PP13 J.P. Papadouris K.E. Papadakis, Astrophys. Space Sci. 344 (2013) 21.

S67 V. Szebehely, Theory of Orbits, Academic Press, New York, 1967.

Z17 E.E. Zotos, Astrophys. Space Sci. 362 (2017) 190.

Z18a E.E. Zotos, Astrophys. Space Sci. 363 (2018) 113.

D10 C.N. Douskos, Astrophys. Space Sci. 326 (2010) 263.

ML01 A.E. Motter, Y.C. Lai, Physical Review E 65 (2001) 015205.

SASL06 J.M. Seoane, J. Aguirre, M.A.F. Sanjuán, Y.C. Lai, Chaos 16 (2006) 023101.

SS08 J.M. Seoane, M.A.F. Sanjuán, Physics Letters A 372 (2008) 110.

AVS01 J. Aguirre, J.C. Vallejo, M.A.F. Sanjuán, Physical Review E 64 (2001) 066208.

AVS09 J. Aguirre, R.L. Viana, M.A.F. Sanjuán, Reviews of Modern Physics 81 (2009) 333.

DWGGS16 A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin, M.A.F. Sanjuán, Scientific Reports 6 (2016) 31416.

PTVF92 H.P. Press, S.A., Teukolsky, W.T., Vetterling, B.P., Flannery, Numerical Recipes in FORTRAN 77, 2nd Ed., Cambridge University Press, Cambridge, USA, 1992.

W03 S. Wolfram, The Mathematica Book, Fifth Edition, Wolfram Media, Champaign, 2003.

R75 E.J. Routh, Proc. Lond. Math. Soc. 6, (1875) 86.




