On the axisymmetric restricted ve-body problem within the frame of variable mass: the convex case

The present manuscript reveals the existence and stability of the libration points (LPs) in the axisymmetric restricted ﬁve-body problem (AR5BP) in which the mass of the ﬁfth particle varies with respect to time as per Jeans’ law. It has been assumed that the four bodies P i with masses m i , i = 1 , 2 ,..., 4 (with m 3 = m 4 = ˜ m ) form an axisymmetric convex conﬁguration. The systems of differential equations which govern the motion of the test particle (whose mass is variable) moving under the gravitational inﬂuence of the four primaries, have been presented, indeed, these equations are different from those of the test particle with constant mass in the AR5BP. We have determined the in-plane and out-of-plane LPs along with their linear stability. Further-more, it is observed that the existence of these LPs depends not only on the angle parameters α and β but also on the parameters occur due to variation in mass namely γ ( 0 < γ ≤ 1 ) and σ ( 0 ≤ σ ≤ 2 . 2, a proportionality constant occurs in Jeans’ law). Moreover, we have investigated the evolution of the regions of possible motion as function of variable mass parameters where the ﬁfth body can move freely. The topology of the basins of convergence (BoC) linked with the LPs as function of σ are unveiled by deploying the bivari-ate version of the Newton-Raphson (NR) iterative method, and, in order to measure the uncertainty of the basins, the basin entropy is also evaluated. The co-relation between the domains of the basins of convergence (DoBoC) and the required number of iterations to achieve predeﬁned accuracy, and also with the corresponding probability distributions are illustrated.


Introduction
One of the most fundamental problem in the field of celestial mechanics is to determine the motion of the n bodies, with the assumption that the mutual Newtonian gravitational forces act between them. The set of masses m i of P i bodies, where the barycentric position vectors are r i , i = 1, 2, ..., n, the equations of motion are: where G is the Gaussian constant of gravitation and r i j = |r j − r i | is the distance between P i and P j bodies. In general, the system (1) is unsolvable for more than n = 2 bodies, therefore, the researchers and scientists discuss the important subclass, i.e., the central configuration. Moreover, many research papers and results have also been published on the different type of central configurations of four-body problem (see Hampton and Moeckel (2006)).
In the recent paper, the central configuration of four bodies is analyzed byÉrdi and Czirják (2016) where two of the bodies are on the axis of symmetry where the positions of the bodies on the axis of symmetry are shown by the angle co-ordinates with respect to the outer bodies. The authors have shown that there exist two types of configurations: (i) the convex configuration and (ii) the concave configuration, further the concave configuration is divided into two categories, the first concave configuration and second concave configuration.
Mathematically, configuration is known as central if there exist Λ ∈ R,Λ > 0 such thatr i = −Λ r i for all i = 1, 2, ...n. Consequently, in the central case, the Eqs. 1 reduce to n ∑ j=1 i j m j r 3 i j (r j − r i ) = −Λ r i , i = 1, 2, ..., n, further, we assumed that the units of mass, length and time are scaled in such a manner that G = 1. Over the years, the restricted three-body problem (R3BP) remains one of the most well-investigated problems in the field of celestial mechanics and non-linear dynamical system. In the planar problem of three bodies where the primaries move under the Newtonian gravitational law, there exist only two permanent configurations: (i) the collinear (or Eulerian) configuration and (ii) the triangular (or Lagrangian) configuration. The permanent configuration can be defined as the ratio of the distance between the primaries is constant at the time of motion. In the collinear configuration, the primaries lie on a straight line whereas in the triangular configuration, the primaries are lying at the vertices of an equilateral triangle. The various novel results in the study of the R3BP under the different perturbations have motivated the scientists and researchers to extend the study of the dynamical system of restricted problem of three bodies into the restricted problem of four bodies. As one moves from the R3BP model to the restricted four-body problem (R4BP) model, the new challenges and problems increases manifold. The R4BP can be retrieved from the R3BP just by introducing one more primary namely P 3 . Previously, the dynamical system of the restricted three-body problem has been studied largely, and numerous models have been suggested for the research focusing the behavior of real celestial systems.
In Ref. Baltagiannis and Papadakis (2011) have studied the R4BP and discussed the motion of the fourth body with negligible mass, equilibrium points (EPs), and their stability. Though, Jeans (1928) and Meshcherskii (1952) are the pioneers in the field of variable mass, but the work of other authors cannot be underestimated. Srivastava and Ishwar (1983), Singh and Ishwar (1984), and Singh and Ishwar (1985) have done impressive work on the problem of restricted three-body with decreasing mass. Abouelmagd et al. (2015b) have investigated the out-of-plane EPs and the forbidden regions in the restricted three-body problem with variable mass. There are few authors who have applied the concept of variable mass in the R4BP and R5BP like: Mittal et al. (2016), Mittal et al. (2018), Aggarwal et al. (2018), and Suraj et al. (2019a).
In the same vein, one can extend the restricted fourbody problem to the restricted five-body problem. There are various applications of the restricted five-body problem, in celestial mechanics, dynamical astronomy, and galactic dynamics, which have been pre-owned by some researchers like Ollöngren (1988), Papadakis and Kanavos (2007) etc. There are many scientists and researchers who worked on the restricted five-body problem, among them few are as follows: Gao et al. (2017), Suraj et al. (2019b), Suraj et al. (2019c), Suraj et al. (2019d) and Suraj et al. (2019e) along with the references therein.
In this work our main objective is to investigate the relative equilibrium solutions and their stability within the frame of variable mass in the restricted five-body problem for the convex configuration. Further, we have found and drawn the domain of BoC associated with the attractors known as LPs by using the well known NR iterative method for few values of the involved parameters. The basin of attraction (see Croustalloudi and Kalvouridis (2007)) is defined as region in which the set of initial conditions (x, y) converge to a particular attractor.
The ultimate goal of our learn the ropes is to pay attention towards the effect of the parameters on the dynamics of the test particle (fifth particle) in the restricted five-body problem in the convex configuration within a frame of variable mass.
The AR5BP within a frame of variable mass is very important because of its wide applications. This is the main reason, we have studied the restricted five-body problem by taking this configuration. So far, no one is able to determine the solution of the N-body problem with N ≥ 3 in spite of an attractive and burning topic of research. Thus, the researchers put forward their efforts in the restricted N−body problem for N = 3, 4, 5 by taking the different perturbations like oblateness Abouelmagd et al. (2015a), radiation pressure and oblateness Douskos (2010), Zotos (2016), Hill problem with oblateness and radiation Zotos (2017d), Copenhagen problem with oblateness Zotos (2018b), the R3BP with triaxial rigid bodies Elshaboury et al. (2016).
The above mentioned ideas motivate us to introduce the effect of variable mass within a frame of the AR5BP in the convex configuration. Insertion of the variable mass of the infinitesimal body has increased challenges manifold in spite of that there are many results and observations exist. Therefore, this problem is one of the worth studying problems.
This paper is organized as follows: a review of the literatures has been presented in Section 1. The intrinsic properties of the dynamical system along with the equations of motion of the fifth particle are discussed in Section 2. In Section 3, the parametric evolution of the in-plane and outof-plane LPs are drawn along with their numerical results whereas, in Section 4, the zero velocity surfaces as a function of energy constant and parameters, occur due to variation in mass, have been drawn. The numerical simulation of the LPs with the help of basins of attractions and basins entropy have been presented in Section 5, whereas the discussion and conclusions are drawn in Section 6.

Properties of the dynamical system
We can not use Newton's second law for the equation of motion within a frame of the variable mass system directly as it is true only for the constant mass system (see Plastino and Muzio (1992)). Therefore, the general equation of motion for the variable mass system can be written as where F ext is the net external force on the body, ϑ rel = ν − ϑ is the relative velocity of the escaping or incoming mass with respect to the body and ϑ is the velocity of body in the inertial frame whereas, ν is the velocity of the escaping or incoming mass to the body. In case of rocket, the velocity ϑ is called the effective exhaust velocity. Now, we consider the dynamical model of the AR5BP with variable mass. As a result, we consider that the five bodies P i of masses m i , i = 1, 2, ..., 5, have their respective position vectors as R i from the origin in the inertial frame OUVW .
Let us assume that the fifth body P 5 has variable mass, i.e., the mass which changes with time t, (m 5 = m 5 (t)). If the loss of mass is taken as non-isotropic, we can apply Newton's second law in the modified form as used by Abouelmagd et al. (2015b) to obtain the equation of motion for the fifth body of the infinitesimal mass m 5 in the inertial frame, The convex configuration of the AR5BP where OUVW is fixed and Ouvw is rotating coordinate systems. The E ′ is the distance from point F to P 3 or P 4 , taken as 1 unit.
in the case in which the escaping or incoming mass occurs from n points for the body in the form Here, Eq. 4 represents the equation of motion of the fifth body whose mass varies with respect to time t according to Jeans' law, and ϑ i = ν i −Ṙ 5 where ϑ i and ν i are the relative velocity and velocity respectively of the escaping or incoming mass with respect to the body from the points i(i = 1, 2, 3, ..., n), whereasṘ 5 is the velocity of fifth particle in the inertial frame. The last term of the above equation vanishes in two cases: if the value of the sum is equal to zero or if ν i =Ṙ 5 , which means the loss of mass is isotropic in both the cases. Now, we further assume that the escaping or incoming mass has zero momentum then ν i = 0 ∀ i. If we replace R 5 by R and m 5 by m then Eq. 4 can be rewritten as The Eq. 5 represents the vector equation of motion of the fifth body whose mass varies with respect to time t.
In the present paper, the case of convex configuration is considered in the AR5BP with in a frame of variable mass. According toÉrdi and Czirják (2016), the conditions to be simultaneously satisfied are as follows: Moreover, the quantity ∆ is the angular velocity of the synodic frame and given by the formula and for all other values of the involved parameters see Ref. Gao et al. (2017), Suraj et al. (2019b).

The equations of motion
The dynamical system, which is to be studied, is the AR5BP in convex configuration under the effect of variable mass. In this problem, we have considered four primaries P i which move in the circular orbits around their common center of mass. Further, we have considered that the fifth body has negligible mass m which varies according to Jeans' law.
In the planar motion of the test particle, we consider OUVW as inertial frame of reference, whereas, Ouvw is taken as synodic frame. Let us assume (U i ,V i ,W i ) and (U,V,W ) be the coordinates of the primaries m i and test particle of mass m, respectively in the inertial (sidereal) coordinates system whereas, in the rotating (synodic) coordinates system, the coordinates of the primaries and infinitesimal body are (u i , v i , w i ), (i = 1, 2, 3, 4) and (u, v, w), respectively.
We have considered that two masses m 1 and m 2 are on an axis of symmetry whereas, other two primaries of equal masses m 3 = m 4 =m are situated symmetrically with respect to the axis of symmetry. Further, we have assumed that the origin of both coordinates system is the center of mass of the system such that Ou−axis goes through the first two primaries m 1,2 with positive direction from m 2 to m 1 and Ov−axis is taken as line perpendicular to the Ou−axis and goes through the center of mass of the system.
We consider that the distances from primaries m 1,2 to origin O are a and b respectively, whereas c is the distance from origin to middle point F of the line joining the masses m 3 to m 4 . Therefore, the equations of motion of the fifth body, whose mass varies in the inertial coordinates system, where If the rotating frame is rotating about w−axis with the angular velocity ω, using the relations between the inertial and rotating coordinates, from the Eqs. 8a − 8c, the equations of motion, in the rotating coordinates system, of an infinitesimal body with variable mass, can be written as: where and (v,ǔ,w) is the dimensional coordinate of m 5 . Now, we establish the equations of motion in the rotating coordinate system using dimensionless variables. We obtain the simplest form of the differential equations of motion. Taking E ′ as the unit of length, shown in Fig. 1, the non-dimensional coordinates are The Eqs. 9a − 9c becomë where According to Jeans' law dm dt = −σ m s where σ is the constant coefficient and s ∈ [0.4, 4.4]. Now, introduce the space-time transformations (u = γ −q x, v = γ −q y, w = γ −q z, dt dτ = γ −k ) such that γ = m m 0 , m 0 is the mass of fifth body at the initial time (t = 0). According to Singh and Ishwar (1984) the applicable values for (s, k, q) are s = 1, k = 0 and q = 1 2 . Thus, the components of velocity and acceleration arė where On substituting Eqs. 13a − 13 f into 12a − 12c, we geẗ where The Eqns. 14a − 14c are the equations of motion of the restricted five-body problem when variation of mass of the fifth body is non-isotropic. Let us consider that the variation of mass of fifth body originates from one point only, i.e., n = 1, the Eqns. 14a − 14c become: where The Eqs. 15a − 15c are equations of motion, with variable mass, which are similar to the case of classical form of the R5BP. In order to save the space we are not writing the expression of the involved parameters in Eq. 16 (for the expression of the involved parameters see Ref. Érdi and Czirják (2016), Gao et al. (2017), Suraj et al. (2019b)). Multiplying the Eqs. 15a − 15c byẋ,ẏ andż, respectively, and adding as well as integrating with respect to t, we obtain a Jacobian quasi-integral or an invariant integral relation in the form: whereẋ,ẏ, andż are components of the velocity , whereas the Jacobian constant is represented by C which is conserved.

The existence of the LPs and their stability
The necessary and sufficient conditions for the existence of LPs arė x =ẏ =ż =ẍ =ÿ =z = 0.
The coordinates (x, y, z) of the EPs can be evaluated by solving the system of following equations, numerically where In the same vein, the second-order partial derivatives of the effective potential function, which will be used later on to illustrate the stability along with BoC associated with the LPs, are: It can be observed that the Eq. 19a-19b reduces to the classical case of the axissymetric restricted five-body problem when σ = 0 (i.e., γ = 1) (see Suraj et al. (2019a)). Moreover, the LPs obtained in this case are ceases to be classical ones, since they do not meatẋ 0 =ẏ 0 = 0 when σ 0 and these are not the special solutions of Eq. 15a-15b. But these LPs converted to the classical one naturally if the variable mass problem degenerates to constant mass. Therefore, it is still meaningful and necessary to evaluate and analyze these unconventional LPs.

In-plane LPs
The LPs which lie in xy-plane (z = 0) are called in-plane LPs. The in-plane LPs can be classified as the LPs lie on the x−axis (known as collinear LPs), and lie off the x−axis (known as triangular LPs). In the present paper, we determine the total number of LPs in the coplanar AR5BP with variable mass when the parameters α, β , and γ are fixed and σ is varying. In an attempt to unveil the effect of the parameters, we have considered three sets of values of the parameters α, β , γ, and σ .

The collinear LPs
The collinear LPs are lying on the x-axis. These LPs are obtained from Eqns. 19a,19b,and 19c by taking y = z = 0. The collinear LPs are the roots of the equation There exist three collinear LPs for the values of parameters α = 43.5 • , β = 41 • , γ = 0.1, and σ ∈ [0.1, 2.2] which are denoted by L 1 , L 2 , and L 3 .

The non-collinear LPs
The non-collinear LPs are lying on the xy-plane but not on x-axis. The non-collinear LPs can be obtained from the Eqs. 19a,19b,and 19c by taking z = 0. The Eq. 19c is always true for z = 0, therefore these LPs can be obtained from Eqs. 19a and 19b by taking z = 0. Thus, non-collinear LPs are the roots of the equations: There exist ten, six, or two non-collinear LPs which are labeled by L i (i = 4, ..., 13) when α = 43.5 • , β = 41 • , γ = 0.1, in the following intervals of σ : In Fig. 2, we depict how the intersections of the equations Ω x = 0, (the darker green line) and Ω y = 0 (the magenta line) define the positions of the LPs on the configuration (x, y) plane. The numerically evaluated LPs are represented by black dots while blue dots represent the center of four primaries. It is observed that the total number of the LPs as well as their positions are changed for the different combination of the parameters α, β , γ, and σ . We have observed that there exist either five, nine, or thirteen LPs in total for various combinations of the parameters used. The parametric evolution of positions of the LPs are shown in Fig. 3 when α = 43.5 • , β = 41 • , γ = 0.4, and for varying values of parameter σ ∈ (0, 2.2]. The color codes for the LPs are as follows: L 1 (purple); L 2 (crimson red); L 3 (green); L 4 (sacramento green); L 5 (pink); L 6 (gray); L 7 (yellow); L 8 (ruby red); L 9 (aqua); L 10 (magenta); L 11 (mustard yellow); L 12 (dark yellow); L 13 (blue).
The red dot represents the critical value of σ = 1.82501 when the LPs decrease from thirteen to nine. The main observations can be summarized as follows: † The collinear libration point L 1 is moving from right to left and L 2 and L 3 are moving from left to right as the value of σ increases. Indeed, the libration point L 1 is moving towards the primary m 1 and libration point L 3 is moving towards the primary m 2 , and L 2 is moving to-wards the primary m 1 and away from the primary m 2 (Fig. 3). † The non-collinear LPs L 4,5 ,L 6,7 ,L 8,9 ,L 10,11 ,and L 12,13 are symmetric about x-axis. † The LPs L 4,5,6,7,8,9 move towards the origin whereas the LPs L 10,11,12,13 move away from it for the increasing values of σ .

The out-of-plane LPs
The LPs which lie in xz− plane (y = 0), but not on x−axis, are called out-of-plane LPs. The out-of-plane LPs can be obtained from Eqs. 19a,19b,and 19c when x 0, z 0, and y = 0. There exist two out-of-plane LPs which are symmetrical with xz-plane, denoted by L + z and L − z and shown in   direction of the movement of LPs when the value of σ increases. It can be noticed that the out-of-plane LPs L ± z move towards the x−axis as the value of σ increases.

The stability of the LPs
The libration point is said to be linearly unstable if its motion is a rapid departure from its vicinity. The libration point is supposed to be a stable position if the infinitesimal mass oscillates around it. To study the stability of the libration point, we give displacement to this libration point as: (23) where ξ , η, and ζ are the coordinates relative to (x 0 , y 0 , z 0 ). x 0 , y 0 , and z 0 are time (t) dependent as they are functions of γ. Therefore, we shall follow the procedure discussed in Zhang et al. (2011) to evaluate the stability of the LPs. Thus, expanding the equations of motion 15a, 15b, and 15c up to the first-order terms, with respect to ξ , η, and ζ , we geẗ where A i j represent the second order partial derivatives of Ω , with respect to x, y, and z, evaluated at the respective LPs and their values are given by with ρ 2 0i =x 2 0i +ỹ 2 0i +z 2 0i , The system 24 is a linear, nonhomogeneous system of equations and known as system of variational equation corresponding to the LPs (x 0 , y 0 , z 0 ). This system reduces to a linear homogeneous system of equations for σ = 0 (which is corresponding to constant mass). When σ > 0, the coordinates (x, y, z) of the fifth body vary with time t, and their distances to the LPs (x 0 , y 0 , z 0 ) changes with time in the system Oxyz . Therefore, linear stability of these LPs can't be determined in ordinary method.
As the mass of the test particle varies with time, to find the linear stability, the system 24 is transformed into a set of first-order equations by usinġ ξ = ξ 1 ,η = η 1 ,ζ = ζ 1 , Using inverse transformation of Meshcherskii (1949), and taking then, the variational equations can be written aṡ The origin is not the solution of Eq. 32 if σ > 0. Indeed, origin is the solution corresponding to σ = 0 has been converted into a nontrivial solution (Lu (1990)). To determine the linear stability of this nontrivial solution corresponds to the LPs depend on the boundedness of the solution of linear, homogeneous system of equations associated to Eq. 32. Then linear stability of the LPs can be determined by finding the characteristic roots of the coefficient matrix in the linear segment of Eq. 32. One can simply obtain the characteristic equation linked to the LPs. Therefore, in order to study these LPs numerically we fix the value of the parameter γ and we have determined the linear stability of the LPs, by evaluating the characteristic roots of the coefficient matrix A, numerically. The characteristic equations of matrix A in Eq. 33a is where ϕ 6 = 1, The libration point is claimed to be stable if all the six roots of characteristic equation, i.e., Eq.34, at the respective libration point are pure imaginary. The values of A 11 , A 22 , A 33 , and A 12 are given by the Eqs. 25-30. The characteristic roots of the Eq. 34 have been calculated at various LPs in the range 0 < γ < 1, 0 < σ ≤ 2.2, and fixed values of the angle parameters (coordinates) α = 43.5 • , β = 41 • . It has been observed that there always exist at least one positive real characteristic root at each libration point. Hence, it can be concluded that all the LPs are unstable for all values of σ > 0.

Regions of motion and zero velocity surfaces
Multiplying the Eqns. 15a, 15b, and 15c byẋ,ẏ, andż, respectively, adding and integrating with respect to the time, we get where V 2 =ẋ 2 +ẏ 2 +ż 2 . The integral part of the Eq. 35 can be re-written in the form: Now from the Eqns. 35 and 36, we have The Eq. 37 can be written as: We shall adopt the methodology of Lukyanov (2009), the Eq. 38 satisfies the inequalities For the possible motion we have V 2 ≥ 0. The Eq. 38 can be re-written as: where On substituting the Eq. 39 into the Eq. 40, we have It can easily be verified that the following relations are valid On substituting the Eqns. 43, 44 into the Eqns. 42a, 42b, respectively, the conditions of motion can be written as Therefore, the conditions of possible regions of motion in Eqs. 45a and 45b in more compact form can be written as: where In Eq. 46c, C is the energy constant which depends on the parameter γ and the initial conditions. Thus, we may infer that the regions of possible motions are indeed explained only by one inequality Ω (γ) ≥ C(γ).
If we compare these regions of possible motions with those of the regions of possible motions in the restricted problem of five bodies which evolves with constant mass, we have observed that the analytical definition for their explanation nearly agree. The only dissimilarity is the existence of variable mass parameter. Therefore, regions of possible motions change with time, whereas these regions do not vary for constant masses, since C(γ) =const. Therefore, these regions change in such a way that at each time they coincide with those for constant masses corresponding to the changed value of γ. In Fig. 5, the boundary of the surface is defined by ∂ Ω ∂ γ = 0. It separates the regions where the energy of the fifth particle is increasing or decreasing. The energy of the fifth particle is decreasing around the circular regions near the primaries m 1 , m 2 , m 3 , and m 4 , i.e., the regions are shown by mustard yellow color which corresponds to ∂ Ω ∂ γ < 0 whereas, this energy is increasing in the region outside the primaries, i.e., the regions shown in the light green color which corresponds to ∂ Ω ∂ γ > 0. In Fig. 6, we have drawn the ZVSs as a function of energy constant C which in turn depend on the various combination of angle parameters and the parameters occur due to variation in mass of the test particle. The Fig. 6 has been drawn for the fixed parameters α = 43.5 • , β = 41 • , γ = 0.8, σ = 1.9, and various values of the Jacobi constant C. The panel-(a) has been drawn for the value of the Jacobi integral C = 9.77592, and shows that there exist small circular regions in white color around each of the primaries where the fifth body can freely move whereas it can not move from one primary to any other primaries. In panel-(b), we have decreased the value of Jaco-bian constant C to 8.09926 and observed that the various opening in the vicinity of the LPs L 1,2,8,9 occur, and therefore, the test particle can orbit freely from one primary to other and vice-versa but not able to orbit to the LPs L 2,4,5,6,7 . As we decrease the energy C, the regions of possible motion increase, and consequently, the forbidden regions shrink significantly. There exist only a rectangular area of the forbidden motion consisting the libration point L 2 (see panels-d).
In Fig. 7, we have illustrated the ZVCs as a function of the parameter σ for α = 43.5 • , β = 41 • , γ = 0.1, C = 3.31691. The shaded regions show the forbidden regions where the motion of the test particle is not possible. We have observed that as the value of the parameter σ increases, the regions of possible motion also increase.
Additionally, We have fixed the value of the parameters α, β , γ in the interval for which there exist thirteen LPs and keeping the Jacobian constant C as constant. It is observed that when σ = 0.1 (see Fig. 7a), four oval shaped regions containing the primaries m i exist where the motion of the test particle is possible, however, the test particle is free to  Fig. 5 The sections of F: F = ∂ Ω ∂ γ = 0 (black boundaries), ∂ Ω ∂ γ > 0 (light green region) , and ∂ Ω ∂ γ < 0 (mustered yellow region ). move on entire configuration plane outside the outer circle. Consequently, the fifth body can not communicate from one primary to either of the primaries and the test particle can not move from inner regions to the outer region and vice-versa, i.e., there are no communication channels between them.
In Fig. 7b, when σ = 0.4, the limiting situations at the LPs L 10 and L 11 occur and consequently the test particle can orbit to communicate between the primaries m 1 , m 3 and m 4 inside the interior region but it can not communicate to m 2 as well as to exterior region and vice-versa. Further increase in the value of the parameter σ to 0.8, the forbidden region constitutes two branches containing the LPs L 1,4,5 and the LPs L 2,3,6,7,12,13 , respectively. Moreover, two more channels around the LPs L 8,9 appear so that the test particle can communicate from interior region to the exterior region and viceversa but it can not move around m 2 ( see Fig. 7c). As we further, increase the value of σ , in the panel-d, two more channels in the vicinity of the LPs L 1,3 occur and the forbidden region constitutes three branches containing the LPs L 4 , L 5 and L 2,6,7 , respectively. Now, the test particle can orbit to communicate between the primaries m i , , i = 1, 2, 3, 4 as well as from interior region to exterior region. Moreover, for σ = 1.0 (see panel-e), two more channels around the LPs L 12,13 exist and the forbidden region constitutes five branches in the vicinity of each of the LPs L 2,4,5,6,7 . In the panel-f, it is observed that as the parameter σ increases to 1.9, the forbidden regions shrink significantly and constitute only one branch containing the libration point L 2 . The regions of possible motion are illustrated for the specific values of the parameters α, β , γ for which thirteen LPs exist and observed that as the value of parameter σ increases, the regions of possible motion also increase but even for higher values of σ , the forbidden regions do not disappear completely.
In Fig. 8, we have illustrated the out-of-plane ZVCs for different values of the parameter σ and for fixed values of the energy constant C = 6.3 and parameters α = 43.5 • , β = 41 • , γ = 0.5. For σ = 0.8 (gray), we have observed that the infinitesimal mass can move only in the oval shaped regions around m 1 and m 2 , therefore, the infinitesimal mass can never approach from the primary m 1 to m 2 and viceversa. For σ = 0.905 (purple), we have observed that there is a small opening at the libration point L 1 , so that the infinitesimal mass can move from m 1 to entire (x, z) plane as well as an oval shaped region around m 2 except the lobes shaped forbidden region. For σ = 0.91, 1.2 shown in green and orange color respectively, it has been observed that the ZVCs constitute two branches containing the LPs L + z and L − z , so that the infinitesimal mass can move more freely through the corridors around the primaries m 1 and m 2 and vice-versa. It is observed that as the value of the parameter σ increases, the regions of possible motion also increase but the forbidden regions do not disappear completely in the out-plane case also.

The Newton-Raphson basins of convergence
The NR iterative method has become as one of the most accurate tool to solve the system of non-linear equations. The associated multivariate iterative scheme can be read as: where f (x n ) represents the system of equations, while the corresponding inverse Jacobian matrix is represented by J −1 .
In the present section, we shall only deal with the BoCs on the configuration (x, y) plane. In this dynamical system, the system of equations in the configuration plane (x, y) can be read as: Moreover, for each coordinate (x, y), the multivariate version of the iterative scheme in Eq.47 may be decomposed into two parts as: Here, the values of x and y are x n and y n at the n-th step of the iterative process of the NR method whereas the corresponding partial derivatives of the effective potential function are represented by the subscripts of Ω (x, y). Recently, the above mentioned iterative scheme has been used to unveil the BoC in various types of dynamical systems such as R5BP (e.g., Zotos (2018a)), the R4BP (?, Zotos (2017a), Zotos (2017b)), collinear restricted four-body problem, in the Copenhagen case with a repulsive quasi-homogeneous Manev-type potential (e.g., ?), in the Sitnikov problem (e.g., Zotos (2018c,d)), and in the R3BP (e.g., Zotos (2017c)). The philosophy that works behind the NR iterative scheme can be summarized as follows: the numerical code is activated with an initial condition (x 0 , y 0 ), on the configuration plane, whereas the iterative scheme continues until a libration point is achieved, with the predefined desired accuracy. We state that the iterative scheme converges for a fixed ini-tial condition on the configuration (x, y) plane, only when, it leads to one of the attractors, no matter what the state of stability of the LPs is. The collection of all such initial conditions, which converge to the same LPs (i.e., attractors), constitute the domains of convergence, always referred as the NR basins of attraction.
The NR BoC associated with the LPs plays a crucial role as the most intrinsic properties of the dynamical system are reflected by the attractive domains. The first and second order partial derivatives of effective potential in the iterative formulae (49a, 49b) lead to the fact that they are associated with the dynamics of the fifth particle's orbit together with the stability properties. It can be observed that the first-order partial derivatives of the effective potential are directly associated with the equations of motion of the fifth body, on the other hand, the second-order partial derivatives of effective potential are used to calculate the variational equations. Moreover, in the monodromy matrix of the periodic orbits, these variational equations appear and they are also associated with the stability properties of the fifth particle. These facts provide the sufficient reasons to study the basins of attraction associated with the LPs.
In the following subsections we will illustrate how the parameter σ affects the topology of the NRBoC in the axisymmetric R54BP, by taking three different cases regarding the total number of the LPs which act as attractors. For the classification of the nodes on the configuration (x, y) plane we will use the color-coded diagrams (CCDs), where each pixel is linked with different color, according the final state of the associated initial condition.

Case I: thirteen LPs exist
We start with the case where thirteen LPs exist. We have drawn the BoC associated with the LPs for the fixed values of the parameters α = 43.5 • , β = 41 • , γ = 0.1, and three different values of the parameter σ = 0.1, 1.8, 1.824. In Fig.  9, the BoC, when thirteen LPs exist, have been drawn. Fig.  9a is drawn for α = 43.5 • , β = 41 • , γ = 0.1, and σ = 0.1, where three collinear and ten non-collinear LPs exist.
It can be seen that the well formed BoC spread over the configuration (x, y) plane, in which the domain of the BoC associated with the central libration point L 2 has infinite extent whereas the domains of convergence associated with all the remaining LPs are finite. In addition, the neighborhood of the basin boundaries (i.e., the various local areas) are composed by extremely fractal mixture of initial conditions. It can be noted that the word "fractal" refers to the particular local region which displays a fractal-like geometry. The final state for an initial condition (x 0 , y 0 ) which lie in the vicinity of basins boundaries is highly sensitive. In particular, a small change in the initial conditions leads to entirely different attractors or final state. Thus, it is very difficult to prognosticate from which of the LPs (attractors) each initial condition will be attracted by.
In Fig. 9(a, d, g), it can be observed that the domain of convergence associated with the libration point L 1,3,8,9 (purple, green, ruby red, aqua) looks like exotic bugs with many legs and antennas. The DoBoC associated with the libration point L 4,5,6,7 looks like a mud lamp in vicinity of these LPs whereas in the vicinity of LPs L 10,11,12,13 looks like eyeball shaped region. We observed that, as the parameter σ increases the overall finite DoBoC shrinks and the basins boundaries become highly chaotic which leads to the fact that the predictability of the initial conditions falling inside these regions is next to impossible. In panels-a, d, g, it is also observed that the BoC linked with the LPs L 10,11,12,13 decrease relatively fast from the other DoBoC. Moreover, these BoC are symmetrical with respect to x−axis.
In Fig. 9(b, e, h), the distribution of the associated number N of iterations requires to reach the preset accuracy are depicted, using tones of different colors. It can be noticed that initial conditions lying inside the basins of attraction has fast convergence for N < 10, whereas the value of N is high for the initial conditions which lean in the neighborhood of the basins boundaries. For Fig.9(e, h), we have observed that there exist considerable amount of initial conditions which converge relatively very slow as for these panels the values of σ are very close to critical value.
Further, Fig. 9(c, f, i) has been drawn corresponding to the probability distribution of the iterations required. We define the probability P as follows: if N 0 initial conditions (x 0 , y 0 ) converge, after N iterations, to one of the attractors, i.e., the LPs, then P = N 0 /N t , where N t denotes the total number of nodes in every CCD. In all the panels of Fig.  9(c, f, i), the tails of histograms extended to cover 95% of the associated distributions of the iterations. The rubin red lines depict the most probable number of iterations N * . We further find that the value of N * is constant (N * = 8) for σ = 0.1, 1.8, 1.824 respectively in these panels.

Case II: nine LPs exist
In this case, we consider those values of the parameters for which nine LPs exist. The topology of the BoC for the fixed values of the parameters α = 43.5 • , β = 41 • , γ = 0.1, and various values of the parameter σ = 1.825, 1.829, 2.02, are drawn in Fig. 10. When we slightly increase the value of the parameter σ = 1.825, the number of LPs changes and the LPs L 6,7,12,13 disappear and the BoCs linked with these LPs also disappear and large area of these BoCs are covered by the BoC linked to the central libration point. Moreover, these areas become highly chaotic and composed with the various initial conditions which converge to other attractors. If the parameter σ increases, the regions near the basin boundaries of the domain of the BoC become highly chaotic which lead to the fact that it is almost impossible to predict that which initial conditions falling inside these regions will converge to which attractors. The corresponding number (N) of needed iterations for the coveted accuracy is illustrated in Fig.10(b, e, h), using tones of different colors whereas the probability distribution of iterations is depicted in Fig. 10(c, f, i). It is clear that to obtain the coveted accuracy, the iterative scheme need not requires more than 20 iterations of the vast majority of the initial conditions. A strange phenomenon occur in these panels, the chaotic sea which occur due to the disappearance of LPs L 6,7,12,13 , in this case, the number of required iteration increases significantly. Indeed, the initial conditions falling inside these chaotic sea have extremely slow rate of convergence.
Furthermore, in this case also the average value of the most probable number N * = 8 (see panels-(c, f, i) of Fig. 10). Moreover, it is also evident that there exist a large amount of initial condition for which the number of required iteration N > 40. Fig. 10 The NR basins of attraction on the configuration (x, y) plane for the increasing values of the parameter σ , when nine LPs exist for α = 43.5 • , β = 41 • , and γ = 0.1: frame-(a) σ = 1.825, frame-(d) σ = 1.829, and frame-(g) σ = 2.02. The color codes for LPs L 1,2,...,9 in Frames-(a, d, g) are same as in Fig. 9(a, d, g). Frames-(b, e, h) are the distribution of the corresponding number N of required iterations for obtaining the NR BoC, shown in Fig. 10(a, d, g). Frames-(c, f, i) represent the corresponding probability distribution of required iterations to obtain the NR BoC. The vertical dashed red line indicates, in each case, the most probable number (N * ) of iterations ( Fig. 10(c, f, i)). The black color dots show the positions of the LPs. (Color figure online).

Case III: five LPs exist
We end our analysis of BoC with the case where five LPs exist in the axisymmeric R5BP in which three are collinear and two are non-collinear. In this subsection, we have drawn the BoC for the fixed values of the parameters α = 43.5 • , β = 41 • , γ = 0.1, and two different values of parameter σ = 2.022, 2.2. In Fig. 11(a, c), the DoBoC linked to the central libration point is infinite whereas for other libration it is finite. The BoC linked with the LPs for which the DoBoC are finite looked like the exotic bugs with various legs and antennas. However, these regions occur in the vicinity of the LPs and the remaining configuration plane is occupied by those initial conditions which converge to central libration point. Moreover, a circular bounded area of the BoC is highly chaotic in nature and composed of the mixture of initial conditions which converge to different attractors. Consequently, the final states of the ICs falling inside these chaotic sea are extremely impossible to predict. In Fig. 11(b, e), the corresponding number (N) of required iterations are depicted, using tones of different colors. One can notice that the number of the iteration required to the achieve the desired accuracy is extremely high inside the circular finite DoBoC which, further, decrease as the values of σ increase.
Moreover, the associated probability distribution of iterations is presented in Fig. 11(c, f). It is observed that the most probable number N * is slightly decreases as the values of parameter σ increase from 2.022 to 2.2.

Basin entropy
The analysis of the CCDs in the above section reveals the fact that the vicinity of the basins boundaries are highly fractal. Consequently, the unpredictability of the final state of ICs falling inside these regions is extremely sensitive, indeed, a slight change in the ICs leads to the completely different attractors. However, the term "fractal" simply means that the specific areas show the fractal like geometry. In Daza et al. (2016) a new tool was introduced to measure the uncertainty of the basins. This new tool named as "basin entropy" refers to the geometry of the basins, consequently, describes the unpredictability in the context of BoC.
Let us recall briefly the numerical algorithm to compute the basin entropy. If in a regtangular region R = [−x, x] × [−y, y] of the convergence plane, there exist N A attractors (the LPs which act as attractors), we partitioned the region R into a grid composed of N square boxes where each box can contain 1 and N A attractors. Consequently, we compute the p i j , the probability that inside the box i resulting attractor is j. Moreover, the ICs inside these boxes are completely independent. The Gibbs entropy, of every box i, is given by where n i ∈ [1, N A ] is the number of attractors inside the box i. The entropy of the entire region R, on the configuration (x, y) plane, is evaluated by the sum of the entropies of the N square boxes as On this basis, the so called basin entropy S b , which is the entropy of the total number of the boxes N, can be defined by the following expression This quantitative measure of the basin entropy can provide the answer of the question that: how the one basin is more unpredictable than the other? In previous section, we have observed that there exist fractal structures for the various values of σ in the DoBoC. Consequently, we discuss the basin entropy to quantify the degree of fractality for various values of σ . However, in some cases, we may only aim in the uncertainty of the basins' boundaries, in particular, we frequently wish to find whether a boundary is fractal or not. For this purpose, we can restrict the calculations of the basin entropy to the boxes falling in the boundaries of basins, that is, we can determine the basin entropy only for those boxes N b which contain more than one color, then the basins boundaries entropy is given as: where S is calculated from Eq. 51. The number S bb is the boundary basin entropy, nature of this quantity is different from the basin entropy S. Therefore, considering a sufficient number of boxes in the boundaries, we can affirm that if the boundary basin entropy is larger than ln 2, then the boundary is fractal, which can be expressed as S bb > ln 2. This is a sufficient but not necessary condition which means that there may be fractal boundaries even for S bb < log2.
Using the algorithm stated above, we have evaluated both entropies, i.e., the basin entropy S b and boundary basins entropy S bb for different values of the parameter σ and shown in Table-1. We can notice that the values of S b and S bb are very high when values of σ < 1.85, therefore in this case the basins boundaries are highly fractal. However, the value of S bb is less than ln 2 for the value of σ = 2.022 and 2.2. It is very crucial to note that in the scenario where the nonconverging initial conditions, or the points those tend to asymptotically to infinity, continue, we calculate them as supplementary basins which co-exist with the other basins linked to the LPs.

Discussion and conclusions
The objective of this manuscript is to illustrate the existence, locations and stability of the LPs in the AR5BP under the effect of variable mass in convex-configuration. We have drawn the regions of possible motion where the test particle is free to move. We have illustrated that how the parameter σ influences the locations of the LPs and their linear stability in both the cases i.e., in-plane and out-of-plane cases. The evolution of the topology of the BoC on the configuration (x, y) plane is explored when the values of the parameter σ vary. Moreover, the number of required iterations N to obtain the desired accuracy, is also monitored in order to classification of the nodes. Moreover, the probability distribution is also illustrated. These BoC linked to the LPs play an important role in any dynamical system as they explain how each initial condition in the configuration plane attracted by the attractors (LPs) of the system. Furthermore, both the basin entropy S b and the boundaries of basins entropy S bb are also illustrated as function of the parameter σ varies and depicted in Table-1. This problem so far under the effect of variation of mass of the test particle in the AR5BP has been studied numerically in a systematic and thorough manner. Therefore, the presented outcomes are novel and this is exactly the contribution of our work.
The main outcomes of our results are as follows: † It has been observed that the existence and number of LPs are strongly depend on the variation of the parameters used. There exist either five, nine or thirteen LPs for various values of considered parameters. † In Figs. 6 and 7, we have plotted the ZVSs in the xyplane and have observed that the regions of possible motion for the test particle increase as the energy constant C decreases (see Fig. 6). Additionally, the regions of possible motion also increase as the parameter σ increases (see Fig. 7). † In the xz-plane, if we increase the parameter σ , the regions of possible motion of the test particle increase (see Fig. 8). † None of the LPs are stable for all permissible values of the used parameters. So we conclude that variation in mass of the test particle destroy the stability of the LPs. † All the BoC associated with the libration point L 2 extend to infinity whereas for all other LPs these basins have finite extent. Further, we have observed that the x−axis is the axis of symmetry for all the BoC associated with the LPs. † None of the initial conditions are there which correspond to the non-converging points. † We have further determined that the huge amount of ICs have very slow rate of convergence when the value of σ is very close to its critical value. Indeed, these ICs are not non-converging points, on the contrary they converge to some attractors sooner or later. † we have determined that the multivariate version of NRiterative scheme is relatively very fast(N < 10), for those initial conditions which lie in the neighborhood of the LPs. We also found that there are very slow convergence for those initial conditions which are lying in the vicinity of the basins boundaries. † It may also be noticed that the value of the basin entropy S b is high when the value of σ is close to zero. However, the value of basins boundaries S bb is greater than ln 2 for the lower values of σ which shows that the basins boundaries are highly fractal in nature(see Table  1). Moreover, when value of σ > 2, the value of basins boundaries S bb is less than the values of ln 2. We believe that this happen due to the fact that for these values of σ there exist only five LPs.
All the numerical calculations and the graphical results are done in this paper by using the latest version 12 of Mathematica . It is very interesting to study the topology of the BoC linked to the LPs in the axisymmetric R5BP in the concave configuration when the mass of the test particle is not constant. Therefore, in future course of work we shall study the problem for the concave case. Further, in future work, one can also find the similarities and differences of the BoC linked to the LPs of the dynamical system by using other types of iterative methods in the axisymmetric R5BP in both the cases, i.e., the convex case and concave case. One should investigate the Wada properties of the basins to study more information about the unpredictability of the system. We say that a basin has the Wada property if any of the initial conditions lie on the boundary of one basin, also lie simultaneously on the boundary of three or more basins. Currently, Zotos (2017d) has used the different iterative methods of one-variable, of higher order (2 or more), to reveal the BoC and compare their convergence properties in the Hills problem with perturbations. Thus, it is very interesting, if one tries these numerical methods to solve the system of two non-linear equations. We suppose that all the above-illustrated ideas would reach to useful, and unexpected, results in the study of the LPs, regions of possible motion and in the active field of DoBoC linked with the LPs in the dynamical system where the mass of the test particle is not constant.

Compliance with Ethical Standards
-Funding: The authors state that they have not received any fund.
-Conflict of interest: The authors declare that they have no conflict of interest.

Data Availability Statements
All data generated or analysed during this study are included in this published article (and its supplementary information files).