Revealing the Newton-Raphson basins of convergence in the circular pseudo-Newtonian Sitnikov problem

TL;DR

This study investigates how the Newton-Raphson method's convergence behavior and basin structures vary in the pseudo-Newtonian Sitnikov problem with three and four primaries, highlighting the influence of a transition parameter.

Contribution

It provides a systematic analysis of convergence basins, fractal boundaries, and basin entropy variations in the pseudo-Newtonian Sitnikov problem, which was not previously explored.

Findings

Convergence basins are significantly affected by the transition parameter.

Fractal basin boundaries are identified and characterized.

Basin entropy varies systematically with the transition parameter.

Abstract

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.

Tables

Table 1. Table 1: Aggregated results regarding the values of the most probable number of iterations N ∗ superscript 𝑁 N^{*} , the location parameter a 𝑎 a and the diversity b 𝑏 b , regarding the basin diagrams shown earlier.

Figure	$ϵ$	$N^{\ast }$	$a$	$b$
3a	0	6	$N^{\ast }$	1.15
4c	0.2	5	$N^{\ast }+1$	1.18
4f	0.5	6	$N^{\ast }$	1.80
4i	0.666	13	$N^{\ast }+1$	2.03
5c	2/3	50	$N^{\ast }+1$	1.97
5f	0.67	11	$N^{\ast }+1$	1.92
5i	1	6	$N^{\ast }+1$	1.65
3c	0	6	$N^{\ast }$	1.15
6c	0.4	6	$N^{\ast }$	1.11
6f	0.8	5	$N^{\ast }+1$	1.27
6i	1	5	$N^{\ast }+1$	1.42

Equations

$${\Omega }_3(x,y,z)=\sum _{i=1}^2\frac{m_i}{r_{3i}}-\frac{ϵ}{2c^4}\sum _{i=1}^2\frac{m_i^3}{r_{3i}^3}+\frac{1}{2}\left(x^2+y^2\right),$$

$$r_{31}$$

$$r_{32}$$

$${\Omega }_4(x,y,z)=\sum _{i=1}^3\frac{m_i}{r_{4i}}-\frac{ϵ}{2c^4}\sum _{i=1}^3\frac{m_i^3}{r_{4i}^3}+\frac{1}{2}\left(x^2+y^2\right),$$

$$r_{41}$$

$$r_{42}$$

$$r_{43}$$

$$\ddot{x}-2\dot{y}=\frac{\partial {\Omega }_i}{\partial x},$$

$$\ddot{y}+2\dot{x}=\frac{\partial {\Omega }_i}{\partial y},$$

$$\ddot{z}=\frac{\partial {\Omega }_i}{\partial z},$$

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

$${\Omega }_3(z)=\frac{1}{r_3}\left(1-\frac{ϵ}{8r_3^2}\right),$$

$${\Omega }_4(z)=\frac{1}{r_4}\left(1-\frac{ϵ}{18r_4^2}\right),$$

$$\ddot{z}$$

$$J_z(z,\dot{z})=2{\Omega }_i(z)-{\dot{z}}^2=C_z,$$

$$\displaystyle f_{3}(\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr};\epsilon)$$

$$\displaystyle f_{4}(\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr};\epsilon)$$

$$\displaystyle\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}\left(8\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-3\epsilon+2\right)=0,$$

$$\displaystyle\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}\left(6\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-\epsilon+2\right)=0.$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{i}=\pm\frac{\sqrt{3\epsilon-2}}{2\sqrt{2}},\ \ \ i=1,2,$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{i}=\pm\frac{\sqrt{\epsilon-2}}{\sqrt{6}},\ \ \ i=1,2,$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n}-\frac{f_{i}(\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n};\epsilon)}{f_{i}^{\prime}(\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n};\epsilon)},\ \ \ i=3,4,$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\frac{12\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{3}\left(8\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-5\epsilon+2\right)}{8\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}\left(8\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-6\epsilon+1\right)+3\epsilon-2},$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\frac{3\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{3}\left(18\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-5\epsilon+6\right)}{6\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}\left(6\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-2\epsilon+1\right)+\epsilon-2},$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\frac{12\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{3}}{8\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-1},$$

$$\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n+1}=\frac{9\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{3}}{6\ooalign{$z$\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}^{2}-1}.$$

$$P(N|a,b)=\frac{1}{2b}\{\begin{array}{ll}\exp \left(-\frac{a-N}{b}\right),& \text{if }N<a\\ \exp \left(-\frac{N-a}{b}\right),& \text{if }N\ge a\end{array}.$$

$$S_i=\sum _{j=1}^{m_i}{P}_{i,j}{\log }_{10}\left(\frac{1}{P_{i,j}}\right),$$

$$S=\sum _{i=1}^N{S}_i.$$

$$S_b=\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{m_i}{P}_{i,j}{\log }_{10}\left(\frac{1}{P_{i,j}}\right).$$

$$S_{bb}=\frac{S}{N_b}$$

$$\frac{{\left(m_1+m_2+m_3\right)}^2}{m_1{m}_2+m_1{m}_3+m_2{m}_3}>3{\left(\frac{1+\kappa }{3-\kappa }\right)}^2,$$

Full text

††journal: International Journal of Non-Linear Mechanics

1 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 [1]. The history of the Sitnikov problem, when the primaries move in circular orbit, begins with Ref. [2], 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. [3, 4, 5, 6, 7, 8, 9, 10, 11]).

Several types of perturbations, such as primaries with either prolate [12] or oblate shape [13], as well as the radiation pressure [14], 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 [15, 16, 7]. In addition, the stability of motion in the same system has been investigated in [17], 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. [18] 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. [19]), including also the oblateness of the primaries (see e.g., Refs. [20, 21]), the effects of the radiation pressure (see e.g., Ref. [22]), and families of periodic orbits and bifurcations (see e.g., Ref. [18]). Moreover, the region of motion in the Sitnikov four-body problem when the fourth mass is finite has also been determined [23], 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. [24] 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\geq 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 [25] the first-order post-Newtonian equations of motion, in the restricted three-body problem, have been deduced by several authors (see e.g., Refs. [26, 27, 28]). 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. [29]. Furthermore, in Ref. [30] 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\geq 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. [12, 31, 32, 33, 34, 35, 36]). Very recently, in Ref. [34] 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\geq 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. [37] 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 A). 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. [38, 39, 40]).

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

2 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}=\mu$ and $m_{2}=1-\mu$, respectively, where $\mu=m_{2}/(m_{1}+m_{2})\leq 1/2$ is the mass parameter [41]. The centers of both primaries are located at $(x_{1},0,0)$ and $(x_{2},0,0)$, where $x_{1}=-\mu$ and $x_{2}=1-\mu$.

According to Ref. [29] the effective potential function of the pseudo-Newtonian circular restricted three-body problem, with only the first correction terms, is given by

[TABLE]

where of course

[TABLE]

are the distances of the test particle from the two primary bodies. 2. 2.

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}=\sqrt{3}/3$, $x_{2}=-\sqrt{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

[TABLE]

where once more

[TABLE]

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

In both cases, the parameter $\epsilon$ acts as a transition parameter from classical Newtonian to pseudo-Newtonian dynamics, while its values lie in the interval $\epsilon\in[0,1]$. Indeed, when $\epsilon=0$ both effective potential functions, $\Omega_{3}(x,y,z)$ and $\Omega_{4}(x,y,z)$, are reduced to the corresponding versions of classical Newtonian dynamics. On the other hand, when $\epsilon=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 $(\omega)$, 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

[TABLE]

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

[TABLE]

where $\dot{x}$, $\dot{y}$, and $\dot{z}$ 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 $\mu=1/2$ we have

[TABLE]

where $r_{3}=\sqrt{z^{2}+1/4}$. Similarly, the potential function of the circular Sitnikov pseudo-Newtonian four-body problem reads

[TABLE]

where $r_{4}=\sqrt{z^{2}+1/3}$. It is evident that Eqs. (7) and (8) 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. 1 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

[TABLE]

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

[TABLE]

with $i=3,4$.

3 Parametric variation of the equilibrium points

From now on, the $z$ coordinate is considered as a complex variable and it is denoted by

$z$

  , thus following the approach which was successfully used in Refs. [12, 42, 35, 43]. At this point, it should be explained that the transition to complex variables is imperative. This is because, as it was presented in Ref. [44], 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. (9) should be taken equal to zero as

[TABLE]

which are reduced to

[TABLE]

Looking at Eqs. (12) we observe that the root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 is always present, regardless the value $\epsilon$ 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. (12) is a third order polynomial it means that there are two additional roots, given by

[TABLE]

in the case of two primary bodies and

[TABLE]

when three primaries are present.

The nature of these two roots strongly depends on the numerical value $\epsilon$ of the transition parameter. Our numerical analysis reveals that for the circular Sitnikov pseudo-Newtonian three-body problem, along with the \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 root

When $\epsilon<-2/3$ there are two imaginary roots.

• 2.

When $\epsilon=0$ or $\epsilon=2/3$ only the root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 exists.

• 3.

When $\epsilon>2/3$ there are two real roots.

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

When $\epsilon<-2$ there are two imaginary roots.

• 2.

When $\epsilon=0$ or $\epsilon=2$ only the root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 exists.

• 3.

When $\epsilon>2$ there are two real roots.

However since the numerical value of the transition parameter $\epsilon$ cannot exceed 1, real roots are impossible for the circular Sitnikov pseudo-Newtonian four-body problem. It should be noted that $\epsilon=1$ is the maximum allowed value of the transition parameter, corresponding to full pseudo-Newtonian dynamics, while all higher values $(\epsilon>1)$ have no physical meaning. It is seen, that the value $\epsilon=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. 2(a-b) shows the parametric evolution of the roots on the complex plane, when $\epsilon\in[0,1]$, with \mathcal{R}=Re[\ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}] and \mathcal{I}=Im[\ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}]. For the case where two primaries are present we see in panel (a) of Fig. 2 that as long as $\epsilon>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 $\epsilon=2/3$ they collide with the central root $R_{1}$ and only the root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 survives thus having multiplicity 3. When $\epsilon>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 $\epsilon$. Finally, when $\epsilon=1$ we the two real roots $\pm\sqrt{2}/4$ reach their maximum distance from the origin. When three primaries are present it is seen in panel (b) of Fig. 2 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 $\epsilon=1$ we have the imaginary roots $\pm i\sqrt{6}/6$ and their evolution is terminated, thus avoiding the collision with the origin.

4 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

[TABLE]

where \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}_{n} is the value of the

$z$

  at the $n$-th step of the iterative process. By plugging formulae (11) in Eq. (15) we obtain

[TABLE]

for the case of two primary bodies and

[TABLE]

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 \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=a+ib, with $\mathcal{R}=a$ and $\mathcal{I}=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\times 1024$ $(\mathcal{R},\mathcal{I})$ 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_{\rm 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 $\epsilon=0$ are presented in panels (a), for the case of two primaries, and (c), for the case of three primaries, of Fig. 3. 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 $\epsilon=0$ takes the form

[TABLE]

while for the case of three primaries we have that

[TABLE]

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. 3. 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 $\epsilon$ 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.

4.1 Case I: Two primary bodies present and $0<\epsilon<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}(\ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr};\epsilon)=0 has, apart from the \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=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. 4. 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. 4 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 $(\mathcal{R},\mathcal{I})$ 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. 4 indicate the best fit to the right-hand side $N>N^{*}$ of them (more details are given in subsection 4.4).

We see that as soon as $\epsilon>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.

• 2.

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 $\epsilon=0.2$ to 13 when $\epsilon=0.666$.

• 3.

In panel (h), where $\epsilon=0.666$ it is seen that throughout the convergence regions the corresponding iterations are significantly higher with respect to panel (b), where $\epsilon=0.2$. We suspect that this behaviour is directly related with the fact that we approaching the critical value $\epsilon=2/3$.

4.2 Case II: Two primary bodies present and $2/3\leq\epsilon\leq 1$

The next case under consideration involves the scenario where there are two real roots, along with the \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0 root. In the first column of Fig. 5 we present the Newton-Raphson basins of convergence for three values of the transition parameter. In panel (a) of Fig. 5, which corresponds to the critical value of $\epsilon$ 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 $\epsilon>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.

• 2.

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 $\epsilon=2/3$ to 6 when $\epsilon=1$.

• 3.

The basin boundaries become more noisy, which means that the fractality of the several basin boundaries on the complex plane increases.

4.3 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}(\ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr};\epsilon)=0 has two imaginary roots, along with the $R_{1}$ root \ooalign{z\cr\hfil\rule[2.15277pt]{1.99997pt}{0.25832pt}\hfil\cr}=0. The evolution of the geometry of the Newton-Raphson basins of convergence is depicted in the first column of Fig. 6, 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 $\epsilon$ 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.

• 2.

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}$.

• 3.

The most probable number $N^{*}$ of iterations displays a minor reduce from 6, when $\epsilon=0.4$, to 5, when $\epsilon=1$.

4.4 An overview analysis

Since the basin diagrams, presented earlier in subsections 4.1, 4.2, and 4.3, 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 $\epsilon$. Such diagrams are presented in Fig. 7(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. 7 we see that the evolution of all types of the initial conditions is very smooth, except near the critical value $\epsilon=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. [45, 46, 47]). 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

[TABLE]

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\geq 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 1. 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.

5 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. [48, 49]. In Fig. 8 we present a magnification of a local area corresponding to the basin diagram shown earlier in panel (g) of Fig. 4, where $\epsilon=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. 8. 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. 9 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 [50] 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]\times[-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

[TABLE]

where $m_{i}\in[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

[TABLE]

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

[TABLE]

Using the value $\varepsilon=0.005$, as suggested in Ref. [50], we calculated the numerical value of the basin entropy $S_{b}$ of the complex plane, for several values of the transition parameter $\epsilon$. 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. 10 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. 4, 5, and 6, but also from additional levels of the transition parameter $\epsilon$.

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. 10a 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 $\epsilon>2/3$ it climbs up again and (ii) at the full pseudo-Newtonian case, where $\epsilon=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. 8). In order to obtain quantitative results regarding the fractal nature of the basins we computed the boundary basin entropy [50]

[TABLE]

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. 10 we demonstrate the evolution of the boundary basin entropy $S_{bb}$, as a function of the transition parameter $\epsilon$. 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.

6 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 $\epsilon$ 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. 2. 2.

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. 3. 3.

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_{\rm max}$) to final states which are not roots of the system, thus displaying a false convergence. 4. 4.

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. 5. 5.

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 $\epsilon\to 1$ (post-Newtonian dynamics), while on the other hand the basin entropy tends to zero when $\epsilon\to 0$ (classical Newtonian dynamics).

For classifying the nodes on the complex plane we used a double precision FORTRAN 77 code [51]. 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*®* [52] 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.

Appendix A Stability of the equilateral triangular configuration of the four-body problem

According to [53] 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

[TABLE]

where $\kappa$ 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 $\kappa=2$) but also in the pseudo-Newtonian problem (where $\kappa=4$).