{'fraction_non_alphanumeric': 0.035447431899747014, 'fraction_numerical': 0.015149732305701005, 'mean_word_length': 4.468951093951094, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Basins of attraction of equilibrium points in the planar circular restricted five-body problem

We numerically explore the Newton-Raphson basins of convergence, related to the libration points (which act as attractors), in the planar circular restricted five-body problem (CR5BP). The evolution of the position and the linear stability of the equilibrium points is determined, as a function of the value of the mass parameter. The attracting regions, on several types of two dimensional planes, are revealed by using the multivariate version of the classical Newton-Raphson iterative method. We perform a systematic investigation in an attempt to understand how the mass parameter affects the geometry as well as the degree of fractality of the basins of attraction. The regions of convergence are also related with the required number of iterations and also with the corresponding probability distributions.

are used to determine the dynamical properties found in the Solar System. These problems have plenty of applications in scientific research, not only in planetary physics, Solar System and galactic dynamics but also in various fields of astrodynamics and astrophysics as well. Also, several modifications are proposed by introducing perturbing terms to the effective potential to understand the motion of the infinitesimal mass in the problem of three and four bodies.
The restricted few-body problem is, beyond any doubt, the best tool for modeling the motion of test particles in several types of planetary systems in our Solar System. In particular, knowing the equilibrium points (also known as Lagrange points) of a system is an issue of paramount importance since they are very important especially for plotting the trajectories of spacecrafts. For example, a large number of NASA missions have been, or will be, sent to saddle point L 1 because this position is ideal for making observations of the Sun-Earth system, since there the gravitational attraction of the Earth partially cancels the gravitational attraction of the Sun. Moreover the L 2 saddle point can be thought as a gateway for bodies (e.g., dangerous asteroids that might collide with the Earth) entering the region in the Earth-Moon system. In addition, as it is well known, the Trojan asteroids can be found near the triangular points L 4 and L 5 . On this basis, knowing the equilibrium points of a system, and of course the corresponding basins of convergence, give us very important information regarding the most intrinsic properties of the dynamical system.
A number of articles have been published on planar central configurations of N-bodies mainly for N = 4, 5 and 7. The study of the motion of N point masses moving under their mutual gravitational attraction governed by Newton's gravitational law is always referred as classical N-body problem. The central configurations are authoritative in the N-body problem since they allow us to obtain the homographic solutions, such as the configuration of the N-bodies at any particular point of time, with respect to the inertial barycentric system, which stay similar to itself as the time changes. It is well known that the first three collinear homographic solutions and two another homographic solutions for N = 3, known as the equilateral triangle solutions, were found by Euler in 1767 and by Lagrange in 1772, respectively.
It was Ollöngren (1988) who introduced the gravitational five-body problem to discuss the motion of the fifth body of negligible mass, in comparison to other bodies. He has supposed that the three bodies with equal masses revolve on the same plane, around their gravitational center in circular orbit under their mutual gravitational pull. In addition, a mass of β > 0 times the mass of one of the three primary bodies is situated at the centre of mass. For β = 0, this particular case of the restricted five-body problem is reduced to the restricted four-body problem. He has obtained nine libration points in total, where three of them become stable when β > 43.18, while for smaller values all the libration points are linearly unstable. Kulesza et al. (2011) discussed the restricted rhomboidal five-body problem and showed that the total number of equilibrium points depends on the ratio of the semi-diagonals and there can be eleven, thirteen or even fifteen libration points, which are all unstable. In the same vein, Marchesin & Vidal (2013) studied the spatial restricted rhomboidal five-body problem and also studied the horizontal stability of its periodic solutions. In this coplanar problem of five bodies, the four primaries move two by two in circular or-bits, while the center of mass is taken as the origin, so that the primaries always maintain a rhombus configuration. Papadakis & Kavanos (2007) numerically explored the restricted five-body problem, when some or all of the primary bodies are sources of radiation, thus extending the work of Ollöngren. They revealed that the number of collinear libration points of the system depends on the radiation factors and of course on the mass parameter. Recently, Gao et al. (2017) studied the motion of the fifth body with infinitesimal mass in the axisymmetric restricted five-body problem, when the four primaries are maintaining the axisymmetric central configurations discussed inÉrdi & Czirják (2016). They revealed that there exist at most fifteen equilibrium points for convex configuration, while thirteen libration points exist for concave configuration, depending upon the angle coordinates. Moreover, they showed that all the libration points are linearly unstable for all the combinations of the angle coordinates, in three axisymmetric configurations.
In the present work we shall use the mathematical model of the restricted five-body problem introduced in Ollöngren (1988). Our aim is to reveal how the geometry as well as the shape of the Newton-Raphson basins of attraction are influenced by the mass parameter. Our paper has the following structure: the most important properties of the dynamical system are presented in Section 2. The parametric evolution of the position as well as of the stability 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 basin entropy evolves, as a function of the mass parameter, while in the next section we provide the main concluding remarks. Our paper ends with Section 7, where we emphasize related aspects which will be incorporated in future works.

Description of the mathematical model
In the circular restricted five-body problem (CR5BP) the four primaries, P 0 , P 1 , P 2 , and P 3 , move on coplanar circular orbits around their common center of gravity. We assume that the fifth body has a significantly smaller mass related to the masses of the primary bodies (m * m 0 , m 1 , m 2 , m 3 ). On this basis, the fifth body acts as an infinitesimal test particle and therefore it does not influence the circular motion of the four primaries.
For modeling the planar motion of the fifth body we choose a rotating frame of reference in which the origin O(0, 0) coincides with the center of the mass of the primary bodies. The dimensionless masses of the primaries are m 0 = β m, m 1 = m 2 = m 3 = m = 1, while the positions of their centers are: (x 0 , y 0 ) = (0, 0), (x 1 , y 1 ) = (1/ √ 3, 0), (x 2 , y 2 ) = (−x 1 /2, 1/2), and (x 3 , y 3 ) = (x 2 , −y 2 ). In  we present the planar configuration of the CR5BP. It is seen that the three bodies of mass m are located at the vertices of an equilateral triangle with side a = 1, while the fourth primary, with mass β m, is located at the center of the equilateral triangle.
Looking at Fig. 1 it becomes evident that the CR5BP admits a 2π/3 symmetry. Indeed, if we rotate the primary symmetry x axis (y = 0), successively through an angle of 2π/3 we obtain the two additional lines of symmetry of the system, that is y = √ 3 and y = − √ 3. According to Ollöngren (1988) the time-independent effective potential function of the CR5BP is where while are the distances of the fifth body from the four primaries. It is interesting to note that when β = 0 (which means that the central primary P 0 is absent) the system is automatically reduced to the circular equilateral restricted four-body problem.
The equations of motion, in synodic coordinates, of the test particle read where In the same vein, the partial derivatives of the second order, which will be needed later for the multivariate Newton-Raphson iterative scheme, are the following , The Jacobi integral of motion, corresponding to the system of differential equations (4) is given by whereẋ andẏ are the velocities, corresponding to coordinates x and y, respectively, while C is the numerical value of the Hamiltonian which is conserved.

Parametric evolution of the equilibrium points
For comparison reasons we can easily define a mass parameter µ = 1/(1 + β), similar to the mass parameter of the classical restricted three-body problem. Therefore we have that µ(0, 1] when β ∈ [0, ∞). The necessary and sufficient conditions, which must be fulfilled for the existence of coplanar equilibrium points, arė x =ẏ =ẍ =ÿ = 0.
(8) The corresponding coordinates (x, y) of the libration points can be determined by solving numerically the system of the first order derivatives The total number of libration points in the CR5BP is a function of the mass parameter µ. In particular -When µ ∈ (0, 0.98617275] there exist nine equilibrium points: three collinear and six non-collinear points (see panel (a) of Fig. 2). -When µ ∈ [0.98617276, 1) there exist fifteen equilibrium points: five collinear and ten non-collinear points (see panel (b) of Fig. 2). -When µ = 1 there are ten libration points, as in the restricted equilateral four-body problem with three equal masses.
The value µ * = 0.98617276 is a critical value of the mass parameter, since it delimits the point where the number of the equilibrium points changes. The intersections of the nonlinear equations Ω x = 0, and Ω y = 0 define the positions of the equilibrium points. Fig.  2(a-b) illustrates how these equations pinpoint, in each case, the location of the libration points, when (a): µ = 0.5 and (b): µ = 0.995. In the same diagram we explain the numbering, L i , i = 1, ..., 15, of all the equilibrium points.
In Fig. 3 we present the parametric evolution of the positions of the equilibrium points, on the configuration (x, y) plane, when µ ∈ (0, 1]. It is seen that as soon as µ is just Fig. 3 The parametric evolution of the positions of the equilibrium points, L i , i = 1, ..., 15, in the CR5BP, when µ ∈ (0, 1]. The arrows indicate the movement direction of the equilibrium points as the value of the mass parameter increases. The big black dots pinpoint the fixed centers of the primaries, while the small black dots (points A, B, and C) correspond to µ → 0, µ = µ * , and µ = 1, respectively. (Color figure online).
above zero three pairs of equilibrium points emerge from the centers P 1 (L 1 and L 2 ), P 2 (L 6 and L 8 ), and P 3 (L 7 and L 9 ) of the primaries. When µ = µ * three additional pairs of libration points appear. As the value of mass parameter grows the equilibrium points evolve following two different patterns. More precisely, L 1 , L 8 , L 9 , L 11 , L 12 , and L 14 move towards the center (0, 0), while all the other points move away from the center. Our analysis suggests that L 1 , L 8 , L 9 , L 11 , L 12 , and L 14 collide with the origin when µ = 1. We observe that all the equilibrium points evolve along the axes of symmetry y = 0, y = √ 3 and y = − √ 3. At this point, we would like to note that the centers of all the primary bodies are completely unaffected by the shift of the mass parameter.
Knowing the exact positions (x 0 , y 0 ) of the equilibrium points, we can easily determine their linear stability, through the nature of the four roots of the characteristic equation.
Our computations indicate that the vast majority of the libration points are always unstable, when the mass parameter µ varies in the interval (0, 1]. Only L 3 , L 4 , and L 5 can be stable however only for extremely small values of the mass parameter and especially when µ ∈ (0, 0.00226341].

The basins of attraction
There is no doubt that the most well-known numerical method for solving systems of nonlinear equations is the famous Newton-Raphson method. This method is applicable to systems of multivariate functions f (x) = 0 through the iterative scheme where f (x n ) denotes the system of equations, while J −1 is the corresponding inverse Jacobian matrix. In our case Eqs.
(9) describe the system of the differential equations. Decomposing the above-mentioned iterative scheme we obtain the following iterative formulae for each coordinate where x n , y n are the values of the x and y coordinates at the n-th step of the iterative process. The philosophy behind the Newton-Raphson method is the following: The numerical code is activated when an initial condition (x 0 , y 0 ), on the configuration plane, is given, while the iterative procedure continues until an equilibrium point (attractor) is reached, with the desired predefined accuracy. If the particular initial condition leads to one of the libration points of the system it means that the numerical method converges for that particular initial condition. At this point, it should be emphasized that in general terms the method does not converge equally well for all the available initial conditions. The sets of the initial conditions which lead to the same attractor compose the so-called Newton-Raphson basins of attraction or basins of convergence or even attracting domains/regions. Nevertheless, it should be clarified that the Newton-Raphson basins of convergence should not be mistaken, by no means, with the basins of attractions which are present in systems with dissipation.
From the iterative formulae of Eqs. (11) it becomes evident that they should reflect some of the most intrinsic properties of the Hamiltonian system. This is true if we take into account that they contain the derivatives of both first and second order of the effective potential function Ω(x, y).
A double scan of the configuration (x, y) plane is performed for revealing the structures of the basins of attraction. In particular, a dense uniform grid of 1024 × 1024 (x 0 , y 0 ) nodes is defined which shall be used as initial conditions of the iterative scheme. Evidently, the initial conditions of the centers of the primary bodies are of course excluded from all the grids because for these initial conditions the distances r i , i = 0, ..., 3 to the respective primaries are equal to zero and therefore several terms, entering formulae (11), become singular. The number N of the iterations, required for obtaining the desired accuracy, is also monitored during the classification of the nodes. For our computations, the maximum allowed number of iterations is N max = 500, while the iterations stop only when an attractor is reached, with accuracy of 10 −15 .
In the following subsections we will determine how the mass parameter µ affects the structure of the Newton-raphson basins of convergence in the CR5BP, by considering two cases regarding the total number of the attractors (equilibrium points). For the classification of the nodes on the configuration (x, y) plane we will use color-coded diagrams (CCDs), in which each pixel is assigned a different color, according to the final state (attractor) of the corresponding initial condition.

Case I: Nine equilibrium points
Our first case under investigation considers the scenario where nine equilibrium points exist, that is when 0 < µ ≤ 0.98617275. The evolution of the geometry of the basins of attraction, for nine values of the mass parameter, is illustrated in Fig. 4(ai). It is seen that in all cases the configuration (x, y) plane is covered by several well-defined basins of convergence which extend to infinity. Furthermore, in the vicinity of the basin boundaries we observe a highly fractal 1 mixture of initial conditions. This means that the basin boundaries are highly chaotic and therefore the final state (attractor) of an The color code, denoting the nine attractors, is as follows: L 1 (green); L 2 (red); L 3 (blue); L 4 (magenta); L 5 (orange); L 6 (indigo); L 7 (brown); L 8 (cyan); L 9 (yellow); non-converging points (white). (Color figure online). initial condition inside this area is highly sensitive. In particular, even the slightest change of the initial conditions automatically leads to a completely different attractor. Therefore for the initial conditions in the basin boundaries it is almost impossible to predict their final state (equilibrium point).
With increasing value of the mass parameter µ the structure of the configuration (x, y) plane changes drastically. The most noticeable changes can be summarized as follows: -The extent of the attracting domains, corresponding to equilibrium points L 3 , L 4 and L 5 , decreases. -The area of the basins of convergence, corresponding to libration points L 2 , L 6 and L 7 , increases. -The extent of the basins of attraction, corresponding to attractors L 1 , L 8 and L 9 , initially increases, while for µ > 0.9 the tendency is reversed.
In panel (i) of Fig. 4, where µ = 0.986172, we identify a considerable amount of non-converging initial conditions, which mainly surround the central region. However, additional numerical calculations revealed that these initial conditions are in fact extremely slow converging points, which do converge to one of the nine attractors, only after a considerable amount of iterations (N 500). In Fig. 5(a-i) we present the corresponding number N of iterations, using tones of blue. It is seen that initial conditions inside the basins of attraction converge relatively fast (N < 20), while the slowest converging points (N > 20) are those in the vicinity of the basin boundaries. The corresponding probability distribution of the required iterations is given in Fig. 6(a-i). In all plots the tails of the histograms extend so as to cover 95% of the corresponding distribution of iterations. The definition of the probability P is the following: if N 0 initial conditions (x 0 , y 0 ) converge, after N iterations, to one of the libration points then P = N 0 /N t , where N t is the total number of nodes in every CCD. Our analysis suggest that the most probable number N * of iterations (see the red vertical dashed lines in Fig. 6) constantly reduces from N * = 23, when µ = 0.001 to N * = 7, when µ = 0.986172.

Case II: Fifteen equilibrium points
In this case, where µ * ≤ µ < 1, there are fifteen equilibrium points: five on the x axis and ten with y 0. The Newton-Raphson basins of attraction for six values of the mass parameter are given in Fig. 7(a-f). It is seen that the basins of attraction, corresponding to libration points L i , i = 10, ..., 15, have the shape of butterfly wings and they also extend to in-finity, as those of the equilibrium points L i , i = 1, ..., 9, discussed in the previous subsection. Furthermore, we observe that the entire pattern (the overall structure on the configuration plane composed of all the different attracting domains) grows rapidly as the mass parameter tends to 1.
of initial conditions fails to obtain the desired accuracy. In particular, for these initial conditions the iterative scheme reaches an attractor (equilibrium point) with an accuracy of 10 −14 and then the convergence stops, even if we allow the iterative procedure to continue for more than 500 iterations. We believe that this numerical malfunction of the code is strongly related with the fact that we are just above the criti-cal value of the mass parameter, where the basin boundaries are very fractalized. Nevertheless, we count all these initial conditions as regular converging nodes, taking into account that an accuracy of 10 −14 is also sufficient and therefore acceptable.  The distribution of the corresponding number N of iterations, required for obtaining the desired accuracy, is illustrated in Fig. 8(a-f). We see that in most cases the vast majority of the initial conditions converge to one of the attractors (equilibrium points) within the first 25 iterations. The corresponding probability distributions are given in Fig. 9(af). In this case the most probable number N * of iterations increases with increasing value of the mass parameter. Indeed, for µ = 0.986173 we have that N * = 8, while for µ = 0.99999 the value of N * is elevated to 10.
When µ = 1 (which means that the central primary body is missing) the five body problem degenerates to the equilateral restricted four-body problem. This fact is also verified through the corresponding Newton-Raphson basins of convergence. In panel (a) of Fig. 10 we provide the CCD for µ = 1. One can easily observe that the overall structure is completely identical to that of Fig. 5 of Zotos (2017a), where the attracting domains of the planar equilateral restricted four-body problem had been investigated.

An overview analysis
The color-coded convergence diagrams on the configuration (x, y) space, presented earlier in Figs. 4 and 7 provide sufficient information regarding the attracting domains, however for only a fixed value of the mass parameter µ. In order to overcome this handicap we can define a new type of distribution of initial conditions which will allow us to scan a continuous spectrum of µ values, rather than few discrete levels. The most interesting configuration is to set x = y, while the value of the mass parameter will vary in the interval (0, 1]. 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 µ is always the ordinate. Fig. 11 shows the basins of attraction on the (x = y, µ) plane, 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 convergence is shown. In panel (a) of Fig. 11 it can be seen very clearly how the convergence properties of the system change when µ ≥ µ * .

Panel (a) of
It would be very interesting to know what happens around the critical value of the mass parameter µ * , where the number of the convergence basins changes from 9 to 15. In the same vein, useful results could be obtained near µ = 1, where the dynamical properties of the system change drastically, as it degenerates to the restricted four-body problem. In panels (a) and (c) of Fig. 12 we present magnifications of panel (a) of Fig. 11, near the two regions with high interest. Now the transition is much more clear, while the changes on the convergence properties of the (x = y, µ) plane are more evident. The corresponding distributions of the required iterations are given in panels (b) and (d) of Fig. 12, respectively.
So far in our work, the degree of fractality of the several types of two-dimensional convergence planes has been discussed only in a qualitative way. We seen that the highly fractal domains are those mainly located in the vicinity of the basin boundaries in which it is almost impossible to predict the final states of the initial conditions. Inside the convergence areas on the other hand, we can easily predict from which attractor (equilibrium point) each initial condition is attracted by. On this basis, it would be an issue of great importance to provide quantitative results about the degree of fractality of the (x = y, µ) plane of panel (a) of Fig. 11. For measuring the degree of fractality we calculated the uncertainty dimension (Ott, 1993) for different values of the mass parameter µ, thus following the computational approach discussed in Aguirre et al. (2001). At this point, it should be emphasized that the uncertainty dimension is completely independent of the particular set of initial conditions used for its calculation.
In Fig. 13 we depict the evolution of the uncertainty dimension D 0 of the (x = y, µ) plane, as a function of the mass parameter µ. It is seen that D 0 ∈ (0, 1) because the calcula-tion of D 0 was performed for only a "1D slice" of initial conditions. It is interesting to note that as the numerical value of µ varies in the interval (0, 1] the degree of fractality fluctuates, while the highest value of D 0 is observed exactly at the critical value of the mass parameter. On the contrary, D 0 is almost zero (display the lowest possible value) when µ = 1, that is when the system degenerates to the circular restricted four-body problem.

Parametric evolution of the basin entropy
In Daza et al. (2016) a new tool for measuring the uncertainty of the basins has been introduced. This new tool is called the "basin entropy" and refers to the topology of the basins, thus describing the notion of fractality and unpredictability in the context of basins of attraction or basins of escape.
Let us briefly recall the numerical algorithm of the basin entropy. We assume that there are N(A) attractors (equilibrium points) in a certain region R of the configuration plane in our dynamical systems. Moreover, R can be subdivided into a grid composed of N square boxes. Each box of the square grid can contain between 1 and N(A) attractors. Therefore we can denote P i, j the probability that inside the box i the resulting attractor is j. Due to the fact that inside the box the initial conditions are completely independent, the Gibbs entropy, of every box i, is given by where m i ∈ [1, N A ] is the number of the attractors inside the box i. The entropy of the entire region R, on the configuration plane, can be computed as the sum of the entropies of the resulting N boxes of the square grid as S = N i=1 S i . On this basis, the entropy relative to the total number of boxes N, which is called basin entropy S b , is given explicitly by the following expression Using the above-mentioned expressions and also adopting the value ε = 0.005, suggested in Daza et al. (2016), we computed the basin entropy S b of the configuration plane, when the mass parameter lies in the interval µ ∈ (0, 1]. Here it should be clarified that in the case where non-converging points are present we count them as an additional type of basin which coexists with the other basins, corresponding to the equilibrium points. In Fig. 14 we present the evolution of the basin entropy as a function of the mass parameter. At this point, it should be noted that for creating this diagram We see that as long as µ > 0 the basin entropy decreases up to about µ = 0.4, while for higher values of the mass parameter the tendency is reversed. The highest value of S b is observed at the critical value µ * , while the basin entropy displays a sudden drop when µ = 0.986173. However with increasing value of µ it continues to gradually grows as the mass parameter tends to one. Interestingly, at exactly µ = 1 (which corresponds to an entire different dynamical system) the basin entropy is suddenly reduced to minimum.
It would be very useful if we could directly relate the results of Fig. 13, regarding the fractal dimension D 0 with those of Fig. 14, about the basin entropy. Unfortunately, the results of the two figures correspond to entirely different types of distributions of initial conditions. In particular, for computing the fractal dimension we used information from the (x = y, µ) plane (that is a one-dimensional set of initial conditions), while for the calculation of the basin entropy we exploited information from the two-dimensional sets of initial conditions on the configuration (x, y) plane. Nevertheless, we cannot ignore the fact that both dynamical quantities (fractal dimension and basin entropy) suggest that the highest degree of fractalization is observed exactly at the critical value of the mass parameter. This implies that for this particular value of µ there is almost complete lack of predictability in a totally deterministic Hamiltonian system.

Discussion and conclusions
We numerically explored the basins of convergence, related to the equilibrium points, in the planar circular restricted five-body problem. More precisely, we demonstrated how the mass parameter µ influences the position as well as the linear stability of the libration points. The multivariate version of the Newton-raphson iterative scheme was used for revealing the corresponding basins of attraction on the configuration (x, y) plane. These attracting domains play a significant role, since they explain how each point of the configuration plane is attracted by the libration points of the system, which act, in a way, as attractors. We managed to monitor how the Newton-Raphson basins of attraction evolve as a function of the mass parameter. Another important aspect of this work was the relation between the basins of convergence and the corresponding number of required iterations and the respective probability distributions.
To our knowledge this is the first time that the Newton-Raphson basins of attraction in the planar circular restricted five-body problem are numerically investigated in a systematic manner. On this basis, the presented results are novel and this is exactly the contribution of our work.
The following list contains the most important conclusions of our numerical analysis.
1. The stability analysis suggests that all most of the equilibrium points of the system are always linearly unstable. Only L 3 , L 4 , and L 5 can be linearly stable and only for extremely small values of the mass parameter.
2. The attracting domains, associated to all the libration points, extend to infinity, in all studied cases. Moreover, the 2π/3 symmetry of the system is present in all the convergence diagrams on the (x, y) plane. 3. Just before the critical value µ * we detected the existence of non-converging initial conditions. Additional computations revealed that these initial conditions are in fact extremely slow converging points, which do converges to one of the equilibrium points of the system but only after many iterations (N 500). 4. Through the classification of the nodes on the several two-dimensional planes we did not encounter any true non-converging initial conditions. 5. The multivariate Newton-Raphson method was found to converge very fast (0 ≤ N < 10) for initial conditions close to the equilibrium point, fast (10 ≤ N < 20) and slow (20 ≤ N < 40) for initial conditions that complement the central regions of the very fast convergence, and very slow (N ≥ 40) for initial conditions of dispersed points lying either in the vicinity of the basin boundaries, or between the dense regions of the libration points. 6. As the value of the mass parameter increases from 0 to µ * the most probable number of required iterations, N * , was found to decrease, while for µ > µ * the tendency is reversed. 7. Both the fractal dimension and the basin entropy were found to display their global maximum at exactly the critical value of the mass parameter µ * . Therefore we have a strong numerical evidence that for this value the system displays its highest degree of fractalization. On the other hand, for µ = 1 we detected the lowest possible value for both of them.
A double precision numerical code, written in standard FORTRAN 77 (Press et al., 1992), was used for the classification of the initial conditions into the different basins of attraction. In addition, for all the graphical illustration of the paper we used the latest version 11.2 of Mathematica (Wolfram, 2003). Using an Intel Quad-Core TM i7 2.4 GHz PC the required CPU time, for the classification of each set of initial conditions, was about 5 minutes.

Future work
In a series of previous papers we numerically investigated the Newton-Raphson basins of convergence in the circular restricted three-body problem (Zotos, 2016a), as well as in the restricted four-body problem Zotos, 2017a). It is in our future plans to devote a new paper in order to demonstrate in detail all the similarities and differences between the several cases, in an attempt to obtain a general overview regarding the convergence properties and the fractal structures, associated to the circular restricted Nbody problem, with N = 3, 4, 5.
In systems where three or more basins of convergence (or escape) coexist one should know whether these basins verify, or not, the so-called Wada property, that is where three or more basins share the same boundary (Poon et al., 1996, e.g.,). An easy way of determining if the Wada property is verified is by plotting the unstable manifold, of a periodic orbit found in the boundary, and monitoring which basins are crossed by it. These calculation, which will allow us to obtain additional relevant information regarding the predictability of the system, will be also conduced in the future paper.
Another interesting aspect would be to use other types of iterative schemes and compare the corresponding similarities and differences related to the basins of attraction of the equilibrium points of the dynamical system. More precisely, we could use iterative methods of higher order, with respect to the classical Newton-Raphson method of second order. In a recent paper (Zotos, 2017b) we deployed a large collection of iterative schemes, of higher order, for revealing and therefore comparing the convergence properties of the Hill problem with oblateness and radiation pressure. So far, all these iterative schemes work only for solving an equation with one variable. Currently, we are trying to expand all these numerical methods so as to be able to use them for solving a system of two nonlinear equations (such as that of Eq. (9)).
We hope that all the above-mentioned ideas would lead to useful, and perhaps unexpected, results in the very active field of attracting domains of equilibrium points in dynamical systems.