Revealing the basins of convergence in the planar equilateral restricted four-body problem

TL;DR

This study investigates how the mass parameter affects the shape and fractality of convergence basins in the planar equilateral restricted four-body problem, revealing the influence of system parameters on dynamical behavior.

Contribution

It provides a systematic numerical analysis of the Newton-Raphson basins of convergence and their dependence on the mass parameter in this specific four-body problem.

Findings

Mass parameter significantly influences basin shape and fractality.

Convergence regions vary with changes in $m_3$.

Attracting domains correlate with iteration counts.

Abstract

The planar equilateral restricted four-body problem where two of the primaries have equal masses is used in order to determine the Newton-Raphson basins of convergence associated with the equilibrium points. The parametric variation of the position of the libration points is monitored when the value of the mass parameter $m_3$ varies in predefined intervals. The regions on the configuration $(x,y)$ plane occupied by the basins of attraction are revealed using the multivariate version of the Newton-Raphson iterative scheme. The correlations between the attracting domains of the equilibrium points and the corresponding number of iterations needed for obtaining the desired accuracy are also illustrated. We perform a thorough and systematic numerical investigation by demonstrating how the dynamical parameter $m_3$ influences the shape, the geometry and the degree of fractality of the converging regions. Our numerical outcomes strongly indicate that the mass parameter is indeed one of the most influential factors in this dynamical system.

Figures

Tables

Table 1. Table 1: The values of the location parameter μ 𝜇 \mu and the diversity b 𝑏 b , related to the most probable number N ∗ superscript 𝑁 N^{*} of iterations, for all the studied cases shown in Figs. 8 , 11 , and 14 .

Figure	$m_3$	$N^{\ast }$	$\mu$	$b$
8a	0.0001	33	$N^{\ast }$	5.09905946
8b	0.01	17	$N^{\ast }+1$	3.22370099
8c	0.05	13	$N^{\ast }$	2.86554309
8d	0.22	10	$N^{\ast }+2$	3.58862734
8e	0.25	12	$N^{\ast }+2$	7.46437951
8f	0.258	13	$N^{\ast }+2$	12.84964306
8g	0.26	16	$N^{\ast }+5$	11.56009567
8h	0.275	10	$N^{\ast }+6$	17.09014571
8i	0.2882761	7	-	-
11a	0.2882762	7	$N^{\ast }+4$	2.22864018
11b	0.29	8	$N^{\ast }+1$	1.55241568
11c	0.30	8	$N^{\ast }$	1.70743943
11d	0.35	7	$N^{\ast }+1$	1.68094743
11e	0.37	11	$N^{\ast }$	1.66194465
11f	0.38	12	$N^{\ast }$	2.28836276
11g	0.40	9	$N^{\ast }+1$	1.94075381
11h	0.42	7	$N^{\ast }$	1.93663675
11i	0.4402	6	-	-
14a	0.4403	6	$N^{\ast }+1$	3.41485683
14b	0.46	6	$N^{\ast }+1$	2.44498156
14c	0.48	6	$N^{\ast }+1$	2.55876653
14d	0.4999	5	$N^{\ast }+2$	3.18349760

Equations

$$\Omega (x,y)=\sum _{i=1}^3\frac{m_i}{r_i}+\frac{1}{2}\left(x^2+y^2\right),$$

$$r_i=\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2},i=1,2,3,$$

$$x_1=\frac{|K|M}{K},$$

$$y_1=0,$$

$$x_2=-\frac{|K|\left[\left(m_2-m_3\right)m_3+m_1\left(2m_2+m_3\right)\right]}{2KM},$$

$$y_2=\frac{\sqrt{3}}{2}\frac{m_3}{m_2^{3/2}}\frac{\sqrt{m_2^3}}{M},$$

$$x_3=-\frac{|K|}{2\sqrt{M}},$$

$$y_3=-\frac{\sqrt{3}}{2}\frac{1}{m_2^{1/2}}\frac{\sqrt{m_2^3}}{M},$$

$$|K|=m_2\left(m_3-m_2\right)+m_1\left(m_2+2m_3\right),$$

$$M=\sqrt{m_2^2+m_2{m}_3+m_3^2}.$$

$${\Omega }_x(x,y)=\frac{\partial \Omega }{\partial x}=\ddot{x}-2\dot{y}=x-\sum _{i=1}^3\frac{m_i\left(x-x_i\right)}{r_i^3},$$

$${\Omega }_y(x,y)=\frac{\partial \Omega }{\partial y}=\ddot{y}+2\dot{x}=y-\sum _{i=1}^3\frac{m_i\left(y-y_i\right)}{r_i^3},$$

$${\Omega }_{xx}(x,y)=\frac{{\partial }^2\Omega }{\partial x^2}=1+\sum _{i=1}^3\frac{m_i\left[2{\left(x-x_i\right)}^2-{\left(y-y_i\right)}^2\right]}{r_i^5},$$

$${\Omega }_{xy}(x,y)=\frac{{\partial }^2\Omega }{\partial x\partial y}=3\sum _{i=1}^3\frac{m_i\left(x-x_i\right)\left(y-y_i\right)}{r_i^5},$$

$${\Omega }_{yx}(x,y)=\frac{{\partial }^2\Omega }{\partial y\partial x}={\Omega }_{xy}(x,y),$$

$${\Omega }_{yy}(x,y)=\frac{{\partial }^2\Omega }{\partial y^2}=1-\sum _{i=1}^3\frac{m_i\left[{\left(x-x_i\right)}^2-2{\left(y-y_i\right)}^2\right]}{r_i^5}.$$

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

$$\dot{x}=\dot{y}=\ddot{x}=\ddot{y}=0.$$

$$\{\begin{array}{l}{\Omega }_x(x,y)=0\\ {\Omega }_y(x,y)=0\end{array}.$$

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n-J^{-1}f({\mathbf{x}}_n),$$

$$x_{n+1}$$

$$y_{n+1}$$

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

Full text

††journal: Astrophysics and Space Science

1 Introduction

The topic of dynamical systems of few-bodies has always been one of the most fascinating fields in celestial mechanics and dynamical astronomy. Especially these days, with the detection of more than 3500 extra-solar planetary systems (see http://exoplanets.eu, update: November 3, 2016), the few-problem problem strongly attracts the scientific interest.

There is no doubt that one of the most well investigated versions of the few-body problem is the circular or elliptic restricted (or not) three-body problem (Szebehely, 1967). In the same vein, the planar restricted four-body problem describes the motion of a test particle with infinitesimal mass (with respect to the masses of the primaries) moving inside the gravitational field of three primary bodies. There are two main configurations regarding the position of the three primary bodies: (i) the Eulerian configuration, where all three primaries lie on the same axis and (ii) the Lagrangian or triangular configuration, where the three primaries always lie at the vertices of an equilateral triangle. For the latter configuration we have the case of the planar equilateral restricted four-body problem (PERFBP). Usually, for the corresponding configurations we use the term “central configurations” due to the fact that the accelerations of the three primary bodies are proportional to the corresponding radius-vectors, while they are directed toward the common center of gravity (Marchal, 1990).

Over the years, the four-body problem has been used, by many researchers, for several practical applications, such as describing real celestial systems. For example: Van Hamme & Wilson (1986); Kloppenborg (2010) for the epsilon Aurigae system, Robutel & Gabern (2006); Melita et al. (2008) for the Sun-Jupiter-Saturn system, Scharz et al. (2009a) for a system of a star, two massive planets and a massless Trojan, Scharz et al. (2009b) and references therein for a system of a star, a brown dwarf, a gas giant and a massless Trojan, Ceccaroni & Biggs (2010) and references therein and Baltagiannis & Papadakis (2013) for the Sun-Jupiter-Trojan Asteroid, Spacecraft system, Simó et al. (1995); Jorba (2000); de Almeida Prado (2005); Machuy et al. (2007) for the Sun-Earth-Moon system, where the fourth body can be a space vehicle.

A very interesting topic of the four-body problem is the location of the periodic orbits (see e.g., Soulis et al., 2008; Baltagiannis & Papadakis, 2011b; Burgos-García & Delgado, 2013; Álvarez-Ramírez et al., 2015a). Similarly, the determination of the position and the stability of the equilibrium points of the four-body problem is another issue of great importance (see e.g., Simó, 1995; Hadjidemetriou, 1980; Michalodimitrakis, 1981; Leandro, 2006; Papadakis, 2007; Baltagiannis & Papadakis, 2011a; Papadouris & Papadakis, 2013; Álvarez-Ramírez et al., 2015b; Singh & Vincent, 2015, 2016).

In dynamical systems knowing the basins of attraction associated with the equilibrium points is very important since this knowledge reveals some of the most inartistic properties of the system. For obtaining the basins of convergence we use an iterative scheme and we perform a scan of the configuration $(x,y)$ plane in order to determine from which of the equilibrium points (attractors) each initial condition is attracted by. The attracting domains in several types of dynamical systems have been numerically investigated. The Newton-Raphson iterative method was used in Douskos (2010) to explore the basins of attraction in the Hill’s problem with oblateness and radiation pressure, while in Zotos (2016a) the multivariate version of the same iterative scheme has been used to unveil the basins of convergence in the restricted three-body problem with oblateness and radiation pressure. Furthermore, the Newton-Raphson converging domains for the photogravitational Copenhagen problem (see e.g., Kalvouridis, 2008), the electromagnetic Copenhagen problem (see e.g., Kalvouridis & Gousidou-Koutita, 2012), the four-body problem (see e.g., Baltagiannis & Papadakis, 2011a; Kumari & Kushvah, 2014), the ring problem of $N+1$ bodies (see e.g., Croustalloudi & Kalvouridis, 2007; Gousidou-Koutita & Kalvouridis, 2009), or even the restricted 2+2 body problem (see e.g., Croustalloudi & Kalvouridis, 2013) have been studied.

In this paper we shall work as in Zotos (2016a), thus following the same numerical techniques and methodology, and we will try to reveal the Newton-Raphson basins of attraction on the configuration $(x,y)$ plane for special case of the PERFBP where two of the three primary bodies have equal masses.

The present paper is organized as follows: In Section 2 we present the basic properties of the considered mathematical model. In section 3 the parametric evolution of the position of the equilibrium points is investigated as the value of the mass parameter $m_{3}$ varies in predefined intervals. In the following Section, we conduct a thorough and systematic numerical exploration by revealing the Newton-Raphson basins of attraction of the PERFBP with two equal masses and how they are affected by the value of the mass parameter. Our paper ends with Section 5, where the main conclusions of this work are presented.

2 Presentation of the mathematical model

We consider a system of units where the units of length, mass and time are taken in such a way so that the sum of the masses $(m_{1}+m_{2}+m_{3})$, the distance between the primaries $(R)$ and of course the angular velocity $(\omega)$ to be equal to unity. Consequently, the gravitational constant is $G=1$.

Without the loss of generality, we may assume that one of the primary bodies is located on the positive $x$ axis, at the origin of time. For describing the motion of the four-body system, we use axes rotating with uniform angular velocity. Furthermore, the three primaries move always on the $(x,y)$ plane, while their mutual distances remain constant with respect to the time (see Fig. 1). In this study, we consider only circular orbits of the primaries around their common center of gravity.

The effective potential function in a synodic system of coordinates is defined as

[TABLE]

where

[TABLE]

are the distances to the respective primaries.

Under the assumption that the three primaries move always on the $(x,y)$ plane and the axes are chosen in such a way so the $m_{1}$ is on the $x$ axis the coordinates of the primaries are given by the following relations (see also Baltagiannis & Papadakis, 2011b)

[TABLE]

where

[TABLE]

The above-mentioned relations (3) apply for the general case where $m_{1}\neq m_{2}\neq m_{3}$. In the special case where $m_{2}=m_{3}<1/2$ and of course $m_{1}=1-2m_{3}$ the coordinates of the centers of the three primaries are: $P_{1}(x_{1},y_{1})=(m_{3}\sqrt{3},0)$, $P_{2}(x_{2},y_{2})=(\left(2m_{3}-1\right)\sqrt{3}/2,1/2)$, and $P_{3}(x_{3},y_{3})=(x_{2},-y_{2})$.

The equations describing the motion of an infinitesimal mass (test particle) in the usual dimensionless rectangular rotating coordinate system read (Moulton, 1900)

[TABLE]

where dots denote the time derivatives.

In the same vein, the second order derivatives (which will be needed for the multivariate Newton-Raphson iterative scheme) are written as

[TABLE]

The system of differential equations (5) admits the integral of the total orbital energy (also known as the Jacobi integral of motion)

[TABLE]

where $\dot{x}$ and $\dot{y}$ are the velocities, while $C$ is the Jacobi constant which is conserved.

3 Parametric variation of the equilibrium points

It is well known that the necessary and sufficient conditions for the existence of every equilibrium point are

[TABLE]

Therefore, the coordinates of the positions of all the coplanar equilibrium points of the PERFBP can be numerically derived by solving the following system of partial differential equations

[TABLE]

At this point, we would like to emphasize that for the rest of the paper we shall deal only with the case where two of the three primaries have equal masses $(m_{2}=m_{3})$. When $y_{0}=0$ the second of equations (5) is fully satisfied because there are only terms of $x_{2}$, $y_{2}$, $x_{3}$, and $y_{3}$ which cancel each other due to the symmetry of the system. Therefore, for the PERFBP with two equal masses collinear equilibrium points exist for every possible value of the mass parameter $m_{3}$. The coordinates $(x_{0},0)$ of the collinear equilibrium points are obtained from the first of equations (5) for $y_{0}=0$.

In Baltagiannis & Papadakis (2011a) it was shown that the total number of the equilibrium points in the PERFBP with two equal masses in not constant but it depends on the value of the mass parameter $m_{3}$. In particular

• •

When $m_{3}\in(0,0.2882761]$ there are eight equilibrium points: two collinear and six non-collinear points (see panel (a) of Fig. 2).

• •

When $m_{3}\in[0.2882762,0.4402]$ there are ten equilibrium points: four collinear and six non-collinear points (see panel (b) of Fig. 2).

• •

When $m_{3}\in[0.4403,0.5)$ there are eight equilibrium points: four collinear and four non-collinear points (see panel (c) of Fig. 2).

In Fig. 2 we see how the intersections of equations $\Omega_{x}=0,\ \Omega_{y}=0$ define on the configuration $(x,y)$ plane the positions of the equilibrium points when (a): $m_{3}=0.15$, (b): $m_{3}=0.36$, and (c): $m_{3}=0.47$.

In this investigation we shall reveal how the mass parameter $m_{3}$ influences the positions of the equilibrium points, when it varies in the interval $m_{3}\in(0,1/2)$. Our results are illustrated in Fig. 3, where we present the space-evolution of all the equilibrium points on the configuration $(x,y)$ plane. One may observe that as value of the mass parameter $m_{3}$ tends to zero the equilibrium points $L_{7}$ and $L_{9}$, as well as $L_{8}$ and $L_{10}$ move towards the centers $P_{2}$ and $P_{3}$, respectively. It is seen that the libration points $L_{2}$ and $L_{3}$ emerge only when $m_{3}\geq 0.2882762$ ($L_{2}$ and $L_{3}$ completely coincide for $m_{3}$ = 0.2882762). Furthermore, at the special case where $m_{3}\to 0.4402$ it was found that the equilibrium points $L_{9}$ and $L_{10}$ collide with each other on the $x$ axis and they disappear for higher values of the mass parameter. Furthermore, it it interesting to note that as the value of the mass parameter tends to 1/2 the libration points $L_{3}$, $L_{4}$, $L_{5}$, and $L_{6}$ tend to collide with each other on the $x$-axis. In fact, when $m_{3}=1/2$ only five equilibrium points exist since the PERFBP degenerates to the restricted three-body problem with equal masses (also known as the Copenhagen problem). In addition, for $m_{3}\to 1/2$ the equilibrium points $L_{7}$ and $L_{8}$ as well as the centers of the primary bodies 2 and 3 tend to the vertical $y$ axis.

The center of the primary 1 $(P_{1})$ moves always on the $x$ axis in the interval $(0,\sqrt{3}/2)$, while the centers $P_{2}$ and $P_{3}$ of the primaries with equal masses move on the parallel lines $y=1/2$ and $y=-1/2$, respectively in the interval $(-\sqrt{3}/2,0)$, when the value of the mass parameter varies in the interval $(0,1/2)$. In the PERFBP with two equal masses $(m_{2}=m_{3})$ the $x$ axis $(y=0)$ is the only axis of symmetry (observe in Fig. 3 the symmetry of the parametric evolution of all points with respect to the $x$ axis).

The stability of all the equilibrium points of the PERFBP with two equal masses has been numerically investigated in Baltagiannis & Papadakis (2011a) (see Table 1).

4 The basins of attraction

There is no doubt that the most famous numerical method for solving systems of equations is the Newton-Raphson method. This method is also applicable to systems of multivariate functions $f({\bf{x}})=0$, through the iterative scheme

[TABLE]

where $J^{-1}$ is the inverse Jacobian matrix of the system of differential equations $f({\bf{x_{n}}})$, where in our case it is described in Eqs. (9).

With trivial matrix calculations (see e.g., Appendix in Zotos, 2016b) we can obtain the following iterative formulae for each coordinate

[TABLE]

where $x_{n}$, $y_{n}$ are the values of the $x$ and $y$ coordinates at the $n$-th step of the iterative process, while the subscripts denote the corresponding partial derivatives of first and second order of the effective potential function $\Omega(x,y)$.

The Newton-Raphson algorithm works as follows: an initial condition $(x_{0},y_{0})$ on the configuration plane activates the code and the iterative process continues until one of the equilibrium points of the system is reached, with some predefined accuracy. In most of the cases the successive approximation points create a crooked path line (see Fig. 4). The initial condition may or may not converge to one of the libration points which act as attractors. If the crooked path leads to one of the equilibrium point then the iterative method converges for the particular initial condition. A Newton-Raphson basin of attraction111It should be clarified and clearly emphasized that the Newton-Raphson basins of convergence should not be mistaken, by no means, with the classical basins of attraction which exist in dissipative systems. The difference between the Newton-Raphson basins of convergence and the basins of attraction in dissipative systems is huge. This is true because the attraction in the first case is just a numerical artifact of the Newton-Raphson iterative method, while in dissipative systems the attraction is a real dynamical phenomenon, observed through the numerical integration of the initial conditions. or convergence (also known as attracting region or domain) is composed of all the initial conditions that lead to a specific attractor (equilibrium point).

One may claim that knowing the basins of attraction of a dynamical system is an issue of paramount importance because these attracting regions may reflect some of the most important qualitative properties of the system in question. This can be justified by taking into account the fact that the derivatives of both first and second order of the effective potential function $\Omega(x,y)$ are included in the iterative formulae (11).

For revealing the structures of the basins of attraction on the configuration $(x,y)$ plane we define a dense uniform grid222Needless to say the the initial conditions corresponding to the three centers $(P_{1},P_{2},P_{3})$ of the primaries are excluded from the grid because for these values the distances $r_{i}$, $i=1,2,3$ to the primaries are zero and therefore several terms of the formulae (11) become singular. of $1024\times 1024$ initial conditions (nodes), which will be used as the initial values of the numerical algorithm. The iterative procedure begins and stops only when an accuracy of $10^{-15}$ regarding the position of the attractors has been achieved. A double scanning of the configuration plane is performed in order to classify all the available initial conditions that lead to a specific equilibrium point (or attractor). While classifying the initial conditions we also record the number $N$ of required iterations in order to obtain the aforementioned accuracy. It is evident that there is a strong correlation between the required number of iterations and the desired accuracy; the better the accuracy the higher the required iterations. In this study we set the maximum number of iterations $N_{\rm max}$ to be equal to 500.

In panel (a) of Fig. 5 we present the Newton-Raphson basins of attraction when $m_{3}=1/3$, which means that all three primaries have equal masses. For each basin of convergence we use different color, while the positions of all the attractors (equilibrium points) are pinpointed by small black dots. All non-converging points are shown in white. We observe that in the case of three equal masses there are two additional axes of symmetry, $y=\sqrt{3}x$ and $y=-\sqrt{3}x$, along with the $y=0$, which exists in the case of two equal masses. In panel (b) 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. Looking the color-coded plot in Fig. 5a we may say that the shape of the basins of convergence corresponding to equilibrium points $L_{4}$, $L_{7}$, and $L_{8}$ look like exotic bugs with many legs and many antennas, while the shape of the basins of attraction corresponding to all other libration points look like butterfly wings.

In the following we shall try to determine how the mass parameter $m_{3}$ influences the structure of the Newton-Raphson basins of attraction, considering three cases regarding the total number and also the type (collinear and non-collinear) of the equilibrium points.

For the classification of the initial conditions on the $(x,y)$ plane we will use modern color-coded diagrams. In these diagrams, each pixel is assigned a specific color according to the particular attractor (equilibrium point). The size of the two-dimensional grids, or in other words the minimum and the maximum values of $x$ and $y$, are chosen differently in each case so as to have a complete view of the basin structures created by the attractors.

4.1 Case I: Two collinear points and six non-collinear points

Our investigation begins with the case where two collinear and six non-collinear equilibrium points are present, that is when $0<m_{3}\leq 0.2882761$. In Fig. 6 we present a large collection of color-coded plots illustrating the Newton-Raphson basins of convergence for several values of the mass parameter $m_{3}$. We observe that the existence of one very large primary body and two small ones substantially influences the structure of the attracting basins, with respect to what we seen earlier in Fig. 5a where all primaries had equal masses. It is seen that well-formed basins of convergence cover the majority of the configuration $(x,y)$ plane. However, the boundaries of all these basins exhibit a highly fractal333When we state that a domain displays fractal structure we simply mean that it has a fractal-like geometry however, without conducting, at least for now, any specific calculations for computing the fractal dimensions as in Aguirre et al. (2009). structure and we may say that they behave as a “chaotic sea”. The meaning of chaos is justified taking into account that if we choose a starting point $(x_{0},y_{0})$ inside these fractal areas we will observe that the choice is highly sensitive. In particular, even a slight change in the initial conditions leads to a completely different final destination (different attractor). This implies that in these areas it is almost impossible to predict from which of the libration points each initial condition is attracted by.

As the value of the mass parameter $m_{3}$ increases the structure of the configuration $(x,y)$ plane changes drastically. We found that the evolution of the structure of the $(x,y)$ plane, with respect to the mass parameter $m_{3}$, does not follow a specific pattern. On the contrary, we observe sudden and completely unpredicted changes which appear even with a slight change of the value of the mass parameter (see e.g. panels (f) and (g), where $m_{3}=0.258$ and $m_{3}=0.26$, respectively).

Looking at Fig. 6(a-i) one may easily observe a very interesting phenomenon related with the extent of the basins of convergence. Indeed, the extent of the basins of convergence corresponding to non-collinear points is always finite, while on the other hand the extent of the attracting domains of the two collinear points, $L_{1}$ and $L_{4}$, extends to infinity. In most of the examined cases it is the $L_{1}$ basins which is infinite. Additional numerical calculations indicate that for $0<m_{3}\leq 0.2195$, $0.248\leq m_{3}\leq 0.259$, and $0.2745\leq m_{3}\leq 0.276$ the attracting domain with infinite area corresponds to libration point $L_{1}$, while for $0.2196\leq m_{3}\leq 0.228$ and $0.2725\leq m_{3}\leq 0.274$ the basins of convergence corresponding to collinear point $L_{4}$ are infinite. For all the other values of $m_{3}$, always in the range $0<m_{3}<0.2882761$, it is nearly impossible to know beforehand which of the two collinear equilibrium points dominates with infinite basins of attraction.

In panel (i) of Fig. 6, where $m_{3}=0.2882761$ (or in other words, equal to the first critical value of $m_{3}$), the basins of convergence are surrounded by a highly chaotic mixture. This mixture is composed of initial conditions of three types: (i) initial conditions attracted by $L_{1}$; (ii) initial conditions attracted by $L_{4}$; (iii) initial conditions for which the multivariate Newton-Raphson method does not converge to any of the attractors. These non-converging initial conditions are in reality extremely slow converging points for which the multivariate Newton-Raphson method requires more than 500 iterations. Indeed, if we increase the maximum allowed number of iteration from 500 to 5000 we will see that all non-converging points finally converge, sooner or later, to either $L_{1}$ or $L_{4}$.

The distribution of the corresponding number $(N)$ of iterations required for obtaining the desired accuracy is provided in Fig. 7(a-i), using tones of blue. It is more than evident that initial conditions inside the basins of attraction converge relatively fast $(N<20)$, while the slowest converging points $(N>50)$ are those in the vicinity of the basin boundaries. In the same vein, in Fig. 8(a-i) the corresponding probability distribution of iterations is presented. The definition of the probability $P$ is the following: if we assume that $N_{0}$ initial conditions $(x_{0},y_{0})$ converge to one of the available attractors, after $N$ iterations, then $P=N_{0}/N_{t}$, where $N_{t}$ is the total number of initial conditions in every color-coded diagram. The blue lines in the histograms of Fig. 8 indicate the best fit to the right-hand side $N>N^{*}$ of them (more details are given in subsection 4.4).

With increasing value of $m_{3}$ the the most probable number $(N^{*})$ of iterations is reduced from 33 when $m_{3}=0.0001$ to 11 when $m_{3}=0.22$. For $m_{3}>0.22$ the most probable number of iteration slightly increases up to $m_{3}=0.26$, while for higher values the tendency is reversed and finally for $m_{3}=0.2882761$ we have $N^{*}=7$. Moreover, it was found that for $m_{3}>0.26$ the range of the probability distribution of the required iterations significantly increases. Indeed, we observe in panel (i) of Fig. 8 that for $m_{3}=0.2882761$ almost the entire range of available iterations, $N\in[0,500]$, is occupied.

4.2 Case II: Four collinear points and six non-collinear points

We continue our exploration with the case where four collinear and six non-collinear libration points exist, that is when $0.2882762\leq m_{3}\leq 0.4402$. The Newton-Raphson basins of attraction for nine values of the mass parameter $m_{3}$ are presented in Fig. 9(a-i). We observe that for values of $m_{3}$ just above the first critical value $(m_{3}=0.2882761)$ the extent of the basins of attraction corresponding to equilibrium points $L_{9}$ and $L_{10}$ are much higher than all the other basins except of course for the attracting domain of libration point $L_{2}$ which extends to infinity. However, as we proceed to higher values of $m_{3}$ the area of the basins of attraction of equilibrium points $L_{9}$ and $L_{10}$ slowly decreases and other attractors take over the configuration space.

In this case $(m_{3}\in[0.2882762,0.4402])$ all four collinear equilibrium points can have infinite basins of attraction, while those corresponding to non-collinear points always have finite area. Our calculations reveal that each attractor with infinite attracting domain dominate in different range of values of the mass parameter $m_{3}$. Being more precisely, we found that

• •

Attractor $L_{1}$ has infinite domains when $0.382\leq m_{3}\leq 0.3945$.

• •

Attractor $L_{2}$ has infinite domains when $0.2882762\leq m_{3}\leq 0.296$, $0.3806\leq m_{3}\leq 0.3819$, and $0.404\leq m_{3}\leq 0.4402$.

• •

Attractor $L_{3}$ has infinite domains when $0.2961\leq m_{3}\leq 0.372$, $0.377\leq m_{3}\leq 0.3805$, and $0.396\leq m_{3}\leq 0.403$.

• •

Attractor $L_{4}$ has infinite domains when $0.3721\leq m_{3}\leq 0.3769$ and $0.3946\leq m_{3}\leq 0.3959$.

When the value of $m_{3}$ approaches the second critical value (0.4402) the area of all finite basins of convergence is reduced and therefore we can easily distinguish, through the color-code diagrams, all the different attracting domains. Indeed, the bug-like structures of the domains corresponding to equilibrium points $L_{4}$, $L_{7}$, and $L_{8}$ and the butterfly-wing shapes of the other basins are again visible. It is interesting to note that in panel (i) of Fig. 9, where $m_{3}=0.4402$, the basins of attraction corresponding to libration points $L_{9}$ and $L_{10}$ heavily suppress the attracting domains of the collinear point $L_{3}$.

Our computations suggest that in this range of values of $m_{3}$ non-converging points are present only when $m_{3}=0.4402$. Even in for this value of the mass parameter however, the relative fraction of non-converging points is so small $(<0.01\%)$ that their presence is almost negligible. Once more, these non-converging initial conditions are in reality extremely slow converging points with $N>500$.

In Fig. 10(a-i) we illustrate the distribution of the corresponding number $(N)$ of iterations required for obtaining the desired accuracy. Looking at panel (i) of Fig. 10 we may observe a very strange phenomenon. So far, the lowest required number of iterations have been identified for initial conditions inside the basins of attraction. However for $m_{3}=0.4402$ (see panel (i) of Fig. 10) it is seen all the initial conditions which are attracted either by $L_{9}$ or $L_{10}$ possess the highest numbers of iteration $(N>50)$ observed in this case. At the moment, we cannot give a logical explanation for this phenomenon (the required number of iterations of initial conditions inside the basins of convergence to be higher than the iterations of the initial conditions in the vicinity of the basin boundaries) which still remains a mystery. Perhaps the fact that this case corresponds to the second critical value of the mass parameter $m_{3}=0.4402$ gives an explanation to this strange and bizarre behavior.

The corresponding probability distribution of iterations is given in Fig. 11(a-i). The the most probable number $(N^{*})$ of iteration starts at 7 for $m_{3}=0.2882762$, then it increases up to $m_{3}=0.38$, while for $m_{3}>0.38$ it decreases until it reaches the lowest observed value, $N^{*}=6$, for $m_{3}=0.4402$.

4.3 Case III: Four collinear points and four non-collinear points

The last case under consideration concerns the scenario according to which the PERFBP with two equal masses has four collinear and four non-collinear equilibrium points $(0.4403\leq m_{3}<1/2)$. In Fig. 12(a-d) we provide the Newton-Raphson basins of convergence for four values of the mass parameter $m_{3}$. In this range of values of $m_{3}$ the changes on the configuration $(x,y)$ plane, due to the variation of the mass parameter, are not so prominent as in the previous two cases. Thus we decided to present in Fig. 12 only four (instead of nine) color-coded convergence diagrams.

As it is seen in Fig. 3 for $m_{3}>0.4402$ the pair of the non-collinear libration points $L_{9}$ and $L_{10}$ disappears. If we compare the color-coded diagrams shown in panel (i) of Fig. 9 and in panel (a) of Fig. 12 we will see that their structure is almost identical with only one major difference: for $m_{3}=0.4403$ the $L_{3}$ attractor has completely assimilated the basins of convergence of equilibrium points $L_{9}$ and $L_{10}$, which are present for $m_{3}=0.4402$. As the value of the mass parameter tends asymptotically to 1/2 the following phenomena take place in the configuration $(x,y)$ plane

• •

The area of the attracting domains corresponding to libration points $L_{1}$, $L_{7}$, and $L_{8}$ increases.

• •

The area of the basins of converging of the attractors $L_{5}$ and $L_{6}$ seems almost unperturbed.

• •

The area of the basins of attractions corresponding to $L_{3}$ and $L_{4}$ decreases.

• •

The convergence domain of collinear point $L_{2}$ is infinite. In fact, in this case $(0.4403\leq m_{3}<1/2)$ only the attractor $L_{2}$ has infinite attracting domain.

• •

The shape (bug-like or butterfly wing-like) of all the basins of attraction remains almost unperturbed.

• •

In this case, there is no numerical evidence of non-converging initial conditions, whatsoever.

Looking at panel (d) of Fig. 12, where $m_{3}=0.4999$, we can numerically verify a statement we made earlier in Section 3. We see that the equilibrium points $L_{3}$, $L_{4}$, $L_{5}$, and $L_{6}$ are very close to one another, while they tend to collide to the limit $m_{3}\to 1/2$. Now if we plot the corresponding four basins of convergence with the same color we will obtain a color-code diagram almost identical to that presented in panel (f) of Fig. 3 in Zotos (2016a), while the only difference will be its orientation (a rotation by 90*∘*).

The corresponding number $(N)$ of required iterations for the desired accuracy is shown in Fig. 13(a-d), while the probability distribution of iterations is presented in Fig. 14(a-d). It was found, that for more than 95% of the examined initial conditions on the configuration $(x,y)$ plane the iterative scheme (11) needs no more than 30 iterations for obtaining the desired accuracy. Furthermore, the average value of required number $(N)$ of iterations remains almost constant, $N^{*}=6$, throughout the range of values of $m_{3}$. Only for $m_{3}=0.4999$ it was measured a little bit lower at $N^{*}=5$.

4.4 An overview analysis

The color-coded convergence diagrams on the configuration $(x,y)$ space, presented earlier in Figs. 6, 9, and 12, provide sufficient information regarding the attracting domains however for only a fixed value of the mass parameter $m_{3}$. In order to overcome this handicap we can define another type of initial conditions which will allow us to scan a continuous spectrum of $m_{3}$ values rather than few discrete levels. The most easy configuration is to set on the two coordinates $(x,y)$ equal to zero, while the initial value of the other coordinate will vary in the interval $[-100,100]$. This technique allows us to construct, once more, a two-dimensional plane in which the $x$ or the $y$ coordinate is the abscissa, while the value of $m_{3}$ is always the ordinate. In panel (a) of Fig. 15 we present the attracting domains of the $(x,m_{3})$ plane when $m_{3}\in(0,1/2)$, while in panel (b) of the same figure the distribution of the corresponding number $(N)$ of required iterations for obtaining the Newton-Raphson basins of attraction is shown. In panels (c) and (d) of Fig. 15 the corresponding results of the $(y,m_{3})$ plane are shown. In both cases, the two critical values of $m_{3}$, which delimit the three cases, are indicated using horizontal black dashed lines.

It is interesting to note that when $0.235<m_{3}<0.2882$ a highly chaotic layer is present in both types of diagrams. Therefore, in this interval the task of knowing beforehand which of the collinear equilibrium points has infinite basins of attraction is next to impossible. Looking at panels (b) and (d) of Fig. 15 it becomes more than evident that the vast majority of initial conditions in this range of values of $m_{3}$ do converge to one of the attractors however, after a relatively high number of iterations $(N>50)$. Finally, it should be noted that for about $m_{3}=0.258$ we detected a small portion of non-converging points. Additional numerical calculations indicated that these non-converging points are in reality extremely slow converging points which need more than 500 of iterations in order to reach to one the attractors. A similar behaviour was also observed for the non-converging points in subsections 4.1 and 4.2. In general terms we may say that in all other regions apart from the interval $m_{3}\in(0.235,28882)$ the multivariate Newton-Raphson method requires no more than 30 iterations in order to converge to one of the attractors.

So far we have discussed the fractality of the several two-dimensional planes only in a qualitative way. More precisely, we seen that the highly fractal areas are those in which we cannot predict from which attractor (equilibrium point) each initial condition is attracted. On the other hand, inside the basins of convergence the degree of fractality is zero and the final state of the initial conditions is well known and of course predictable. At this point we shall provide a quantitative analysis regarding the degree of fractality for the $(x,m_{3})$ and $(y,m_{3})$ planes shown earlier in panels (a) and (c) of Fig. 15, respectively. In order to measure the degree of fractality we have computed the uncertainty dimension (Ott, 1993) for different values of the mass parameter $m_{3}$, thus following the computational method introduced in Aguirre et al. (2001). Obviously, this degree of fractality is completely independent of the initial conditions we used to compute it.

The evolution of the uncertainty dimension $D_{0}$ for both $(x,m_{3})$ and $(y,m_{3})$ planes, as a function of the mass parameter $m_{3}$, is shown in Fig. 16(a-b). The computation of the uncertainty dimension was done for only a “1D slice” of initial conditions of Figs. 15, and for that reason $D_{0}\in(0,1)$. It is interesting to note that in both types of planes for the first critical value of $m_{3}$, that is $m_{3}=0.2882761$, the uncertainty dimension tends to one. This means that for that critical value there is a total fractalization of the corresponding planes and the chaotic set becomes “dense” in the limit. On the contrary, for the second critical value $(m_{3}=0.4402)$ the uncertainty dimension is almost zero. Looking carefully both panels of Fig. 15 we may conclude that the highest degree of fractality is observed when $0.16<m_{3}<0.18$ and especially when $0.2<m_{3}<0.2882761$. On the other hand the lowest degree of fractality is measured in the third case $(m_{3}>0.4402)$, where four collinear and four non-collinear equilibrium points exist.

Before closing this numerical investigation we would like to shed some light to the probability distributions of iterations presented in Figs. 8, 11, and 14. In particular, it would be very interesting to try to obtain the best fit of the tails444By the term “tails” of the distributions we refer to the right-hand side of the histograms, that is, for $N>N^{*}$. of the distributions. For finding the best fit of the tails we tried several single types of distributions (Laplace, Maxwell-Boltzmann, Rayleigh, Pascal, Poisson, etc). Our calculations strongly indicate that in the vast majority of the cases the Laplace distribution is the best fit to our data. The only two cases where the Laplace distribution fails to fit the corresponding numerical data are the cases corresponding to the two critical values of the mass parameter $m_{3}$ (see panels (i) in Figs. 8 and 11).

The probability density function (PDF) of the Laplace distribution is given by

[TABLE]

where $\mu$ is the location parameter, while $b>0$, is the diversity. In our case we are interested only for the $x\geq\mu$ part of the distribution function.

In Table 1 we present the values of the location parameter $\mu$ and the diversity $b$, as they have obtained through the best fit, for all cases discussed in Figs. 8, 11, and 14. One may observe that for most of the cases the location parameter $\mu$ is very close to the most probable number $N^{*}$ of iterations, while in some cases these two quantities coincide. Here we would like to emphasize that the Laplace distribution is only a first good approximation to our data. Additional numerical calculations indicate that if we use a mixture of several types of distributions, instead of a single type of distribution (i.e., the Laplace distribution), the fit is much better. However we feel that this task is out of the scope and the spirit of this paper and therefore we did not pursue it.

5 Concluding remarks

The scope of this research paper was to numerically determine the basins of convergence associated with the equilibrium points. In the PERFBP with two equal masses the number, the position and the type of the libration points strongly depends on the value of the mass parameter $m_{3}$. With the help of the multivariate version of the Newton-Raphson iterative scheme we managed to unveil the extraordinary and magnificent structures of the basins of attraction corresponding to the equilibrium points of the dynamical system. These basins play an important role as they describe how each point on the configuration $(x,y)$ plane is attracted by the libration points which act as attractors. Our numerical exploration revealed how the position of the equilibrium points and of course the structure of the attracting areas are influenced by the mass parameter $m_{3}$. Furthermore, we related the several basins of attraction with the corresponding distribution of the required number of iterations. To our knowledge, this is the first time that such a thorough and systematic numerical investigation, regarding the basins of attraction, takes place in the PERFBP and this is exactly the novelty as well as the importance of the current work.

The main results of our numerical research are the following:

We observed that the change on the value of the mass parameter $m_{3}$ mostly influences the shape and the geometry of the basins of attraction when $0<m_{3}<0.4402$. Indeed, in this range it is almost impossible to know beforehand the structure of the attracting regions, as even a tiny change on the value of $m_{3}$ leads to a complete different structure on the $(x,y)$ plane. For larger values of $m_{3}$ the influence of the same parameter on the attracting regions is much more milder. 2. 2.

In all examined cases the area of the basins of convergence corresponding to the non-collinear libration points is finite. On the other hand, all four collinear equilibrium points can have infinite attracting domains. In particular, it was found that when $0.2882762\leq m_{3}\leq 0.4402$ each collinear attractor dominate in different ranges of the values of $m_{3}$. 3. 3.

The iterative method was found to converge very fast $(0\leq N<15)$ for initial conditions around each equilibrium point, fast $(15\leq N<25)$ and slow $(25\leq N<50)$ for initial conditions that complement the central regions of the very fast convergence, and very slow $(N\geq 50)$ for initial conditions of dispersed points lying either in the vicinity of the basin boundaries, or between the dense regions of the equilibrium points. 4. 4.

In general terms we concluded that the average value of required iterations $(N^{*})$ for obtaining the desired accuracy decreases with increasing value of the mass parameter. In almost all cases, the Newton-Raphson method, for more than 90% of the initial conditions, requires less than 60 iterations to converge to one of the available attractors. 5. 5.

Our calculations strongly suggest that non-converging points on the configuration $(x,y)$ plane exist only when $0<m_{3}<0.4402$. Being more precisely, they appear mostly for values of $m_{3}$ just before the two critical values ($m_{3}=0.2882761$ and $m_{3}=0.4402$). A deeper analysis on these points revealed the fact that these points are not true non-converging points. In reality they are extremely slow converging points which require a huge number of iterations $(N>500)$ in order to reach to one of the attractors. 6. 6.

In the case where eight equilibrium points exist (two collinear and six non-collinear) we observed the highest degree of fractality, especially when $0.2<m_{3}<0.2882761$. On the other hand, in the third case $(m_{3}>0.4402)$, where four collinear and four non-collinear equilibrium points exist, we measured the lowest degree of fractality. 7. 7.

Our tests indicate that our numerical data, corresponding to the histograms with the probability distributions of the required iterations, are best fitted by the Laplace probability distribution function (PDF). Only the cases just before the two critical values of the mass parameter (which have long tails) cannot be fitted well by a Laplace PDF.

A double precision code, written in standard FORTRAN 77, has been deployed for performing all the required numerical calculations regarding the basins of convergence. For the graphical illustration of the paper, we used the latest version 11.0 of Mathematica*®* Wolfram (2003). For the classification of each set of the initial conditions on the several types of two-dimensional planes, we needed about 5 minutes of CPU time using a Quad-Core i7 2.4 GHz PC, depending of course on the required number of iterations. When an initial condition had converged to one of the attractors with the predefined accuracy the iterative procedure was effectively ended and proceeded to the next available initial condition.

Judging by the novel results revealed through the detailed and systematic numerical exploration we believe that we successfully completed our computational task. We hope that our investigation and the corresponding outcomes to be useful in the field of attracting domains in the PERFBP. Taking into account that the current analysis of the case with two equal masses was encouraging it is our future plans to study the general case with three unequal masses. In addition, it would be highly interesting to try and use other types of iterative formulae (of higher order with respect to the classical Newton-Raphson method) and compare the similarities as well as the differences of the structures of the basins of attraction.

Acknowledgments

I would like to express my warmest thanks to the anonymous referee 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.