Equilibrium points and basins of convergence in the triangular restricted four-body problem with a radiating body

TL;DR

This paper investigates how radiation effects influence equilibrium points and basins of convergence in a restricted four-body problem with three primaries in an equilateral triangle, revealing instability and reduction of fixed points at high radiation levels.

Contribution

It extends the restricted four-body problem by incorporating radiation pressure, Poynting-Robertson drag, and solar wind drag, analyzing their impact on equilibrium points and basins of convergence.

Findings

Libration points become unstable if radiation factor exceeds 0.01.

Number of equilibrium points decreases as radiation factor increases.

Radiation effects can destabilize the system's equilibrium in realistic celestial configurations.

Abstract

The dynamics of the four-body problem have attracted increasing attention in recent years. In this paper, we extend the basic equilateral four-body problem by introducing the effect of radiation pressure, Poynting-Robertson drag, and solar wind drag. In our setup, three primaries lay at the vertices of an equilateral triangle and move in circular orbits around their common center of mass. Here, one of the primaries is a radiating body and the fourth body (whose mass is negligible) does not affect the motion of the primaries. We show that the existence and the number of equilibrium points of the problem depend on the mass parameters and radiation factor. Consequently, the allowed regions of motion, the regions of the basins of convergence for the equilibrium points, and the basin entropy will also depend on these parameters. The present dynamical model is analyzed for three combinations of mass for the primaries: equal masses, two equal masses, different masses. As the main results, we find that in all cases the libration points are unstable if the radiation factor is larger than 0.01 and hence able to destroy the stability of the libration points in the restricted four-body problem composed by Sun, Jupiter, Trojan asteroid and a test (dust) particle. Also, we conclude that the number of fixed points decreases with the increase of the radiation factor.

Figures

Tables

Table 1. Table 1: Location of the primary bodies according to the four possible cases.

Case	$x_1$	$x_2$	$x_3$	$y_1$	$y_2$	$y_3$
1	$-$	+	+	0	$-$	+
2	$-$	+	+	0	+	$-$
3	+	$-$	$-$	0	$-$	+
4	+	$-$	$-$	0	+	$-$

Table 2. Table 2: Existence of equilibrium points with the variation of the radiation factor β 𝛽 \beta , for m 1 = m 2 = m 3 subscript 𝑚 1 subscript 𝑚 2 subscript 𝑚 3 m_{1}=m_{2}=m_{3} .

Interval	N. Equilibria	Equilibria
$\beta \in [0.000,0.687]$	10	$L_{1,2,3,4,5,6,7,8,9,10}$
$\beta \in [0.688,0.999]$	8	$L_{2,3,4,5,6,7,8,9}$
$\beta =1.000$	4	$L_{3,4,7,8}$

Table 3. Table 3: Existence of equilibrium points with the variation of the radiation factor β 𝛽 \beta , for different values of m 2 = m 3 subscript 𝑚 2 subscript 𝑚 3 m_{2}=m_{3} .

$𝔪$	Interval	Equilibria
$0.05$	$\beta \in [0.000,0.928]$	$L_{1,2,5,6,7,8,9,10}$
	$\beta \in [0.929,0.999]$	$L_{2,5,6,7,8,9}$
	$\beta =1.000$	$L_{7,8}$
$0.10$	$\beta \in [0.000,0.912]$	$L_{1,2,5,6,7,8,9,10}$
	$\beta \in [0.913,0.999]$	$L_{2,5,6,7,8,9}$
	$\beta =1.000$	$L_{7,8}$
$0.15$	$\beta \in [0.000,0.350]$	$L_{2,3,5,6,7,8,9,10}$
	$\beta \in [0.351,0.360]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.361,0.889]$	$L_{1,2,5,6,7,8,9,10}$
	$\beta \in [0.890,0.999]$	$L_{2,5,6,7,8,9}$
	$\beta =1.000$	$L_{7,8}$
$0.20$	$\beta \in [0.000,0.293]$	$L_{2,3,5,6,7,8,9,10}$
	$\beta \in [0.294,0.608]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.609,0.859]$	$L_{1,2,5,6,7,8,9,10}$
	$\beta \in [0.860,0.999]$	$L_{2,5,6,7,8,9}$
	$\beta =1.000$	$L_{7,8}$
$0.25$	$\beta \in [0.000,0.160]$	$L_{2,3,5,6,7,8,9,10}$
	$\beta \in [0.161,0.816]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.817,0.999]$	$L_{2,3,4,5,6,7,8,9}$
	$\beta =1.000$	$L_{3,4,7,8}$
$0.30$	$\beta \in [0.000,0.751]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.752,0.999]$	$L_{2,3,4,5,6,7,8,9}$
	$\beta =1.000$	$L_{3,4,7,8}$
$0.35$	$\beta \in [0.000,0.644]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.645,0.999]$	$L_{2,3,4,5,6,7,8,9}$
	$\beta =1.000$	$L_{3,4,7,8}$
$0.40$	$\beta \in [0.000,0.430]$	$L_{1,2,3,4,5,6,7,8,9,10}$
	$\beta \in [0.431,0.999]$	$L_{2,3,4,5,6,7,8,9}$
	$\beta =1.000$	$L_{3,4,7,8}$
$0.45$	$\beta \in [0.000,0.999]$	$L_{1,2,3,4,5,6,7,8}$
	$\beta =1.000$	$L_{3,4,7,8}$

Table 4. Table 4: Existence of equilibrium points with the variation of the radiation factor β 𝛽 \beta , for the Sun-Jupiter-Trojan Asteroid system.

Interval	N. Equilibria	Equilibria
$\beta \in [0.0000,0.0029]$	8	$L_{1,2,3,4,5,6,7,8}$
$\beta \in [0.0030,0.9887]$	6	$L_{1,2,3,5,7,8}$
$\beta \in [0.9888,0.9940]$	4	$L_{1,2,5,8,}$
$\beta \in [0.9941,1.0000]$	2	$L_{2,5}$

Equations

$$m\ddot{\mathbb{R}}$$

$$\beta =\frac{F_r}{F_g}=\frac{SA{Q}_{\text{pr}}{R}_1^2}{cGM{m}_1},$$

$$\ddot{x}$$

$$\ddot{y}$$

$$F_x$$

$$F_y$$

$$N=(x-x_1)\dot{x}+(y-y_1)\dot{y},$$

$$r_i={\left[(x-x_i{)}^2+(y-y_i{)}^2\right]}^{1/2};i=1,2,3.$$

$$y_1$$

$$m_1{x}_1+m_2{x}_2+m_3{x}_3$$

$$m_1{y}_1+m_2{y}_2+m_3{y}_3$$

$$(x_2-x_1{)}^2+(y_2-y_1{)}^2$$

$$(x_3-x_2{)}^2+(y_3-y_2{)}^2$$

$$(x_1-x_3{)}^2+(y_1-y_3{)}^2$$

$$x_1$$

$$x_2$$

$$x_3$$

$$y_1$$

$$y_2$$

$$y_3$$

$$K_1=\pm \sqrt{m_2^2+m_2{m}_3+m_3^2};K_2=m_1+m_2+m_3.$$

$$U=\frac{1}{2}\left(x^2+y^2\right)+\frac{m_1}{r_1}(1-\beta )+\frac{m_2}{r_2}+\frac{m_3}{r_3},$$

$$\ddot{x}-2\dot{y}$$

$$\ddot{y}+2\dot{x}$$

$$\dot{\mathcal{C}}=2(1+sw)\left(F_x\dot{x}+F_y\dot{y}\right).$$

$$\mathcal{C}=2U-\left({\dot{x}}^2+{\dot{y}}^2\right).$$

$$x-\frac{m_1(x-x_1)}{r_1^3}(1-\beta )+(1+sw)\frac{\beta m_1(y-y_1)}{c{r}_1^2}$$

$$-\sum _{i=2}^3\frac{m_i(x-x_i)}{r_i^3}=0,$$

$$y-\frac{m_1(y-y_1)}{r_1^3}(1-\beta )-(1+sw)\frac{\beta m_1(x-x_1)}{c{r}_1^2}$$

$$-\sum _{i=2}^3\frac{m_i(y-y_i)}{r_i^3}=0,$$

$$C(t)$$

$$a_1$$

$$a_2$$

$$a_3$$

$$C=2U-2(1+sw)\frac{\beta m_1}{c}\arctan \left(\frac{y-y_1}{x-x_1}\right).$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{\xi }$$

$$\dot{\eta }$$

$$\mathbb{A}=\left(\begin{array}[]{cccc}0&0&1&0\\ 0&0&0&1\\ A_{11}&A_{12}&A_{13}&A_{14}\\ A_{21}&A_{22}&A_{23}&A_{24}\end{array}\right)$$

$$A_{11}$$

Full text

1 Introduction

The $n$-body problem is one of the oldest and most researched problems in astrophysics, it deals with the motion of bodies that are subject to their mutual gravitational attractions (or to other forces as well). The study of the $n$-body problem is essential to the understanding of the motions of celestial bodies and is of particular interest in spacecraft navigation [1]. Since there is no general analytical solution to the $n$-body problem for $n\geq 3$, several simplifications have been introduced with the most prominent being the restricted three and four body problems [2]. In both cases, the mass of one of the bodies (test particle) is negligible in comparison to the others (primaries) such that it does not affect their motion. In the restricted three-body problem the primaries move in a circular or elliptic trajectory (solutions of a two-body problem), while in the restricted four-body problem the primaries move in a straight-line equilibrium configuration or in an equilateral triangle configuration (solutions to the three-body problem) [3].

The restricted four-body problem has many applications in celestial mechanics, dynamical astronomy, and galactic dynamics, mainly motivated by the fact that there are many astrophysical systems that can be roughly approximated to that model, for example, the Sun-Jupiter-Saturn-satellite system, the Sun-Jupiter-Trojan-spacecraft system, the Sun-Jupiter-Earth-satellite system, the Saturn-Tethys-Telesto-spacecraft system [4], or even any of the four-star systems of our galaxy [5]. Aiming to get a better understanding of the problem, several modifications to the effective potential in the classical restricted four-body problem have been proposed, see e.g. [6, 7, 8, 9], where it was found that the introduction of the additional terms significantly modifies the existence, position, and stability of the equilibrium points, and therefore the overall dynamics of the system.

Since it is a well-known fact that in the $n$-body with $n\geq 3$ it is not possible to find explicit formulae for the position of the equilibrium points, numerical methods become the natural and the most efficient way of finding the locations of the fixed points in many different dynamical systems. One of the most popular methods to find roots of multivariate functions is the Newton-Raphson algorithm, this method allows us to get a global picture of the set of initial conditions that lead to a particular fixed point. The final picture produced after evolving (via the root-finding algorithm) the whole set of initial conditions in a certain region $\Omega$ of the configuration space is called the basin of convergence (see e.g. [10, 11]).

During the past few years, the Newton-Raphson basins of convergence have been extensively investigated in many different versions of the restricted four-body problem. As characteristic examples of these works, may be mentioned the one by Baltagiannis & Papadakis [12], who determined the position and stability of the equilibrium points in the equilateral triangle configuration of the four-body problem using different combinations of primaries, and the study carried out by Zotos [13], who extended the study of the basins of convergence for cases not considered by Bataglianis. Moreover, Suraj et. al. [14] studied the existence and locations of libration points along with the Newton-Raphson basins of convergence for the same model, but using triaxial rigid bodies as primaries, and later studied the effect of small perturbations in the inertial forces on the Newton-Raphson basins of convergence [15]. Further, the linear equilibrium configuration of the restricted four-body problem was considered by Zotos in order to determine how the mass parameter and angular velocity influence the geometry of the Newton-Raphson basins of convergence [16], while the photo-gravitational restricted four-body problem with variable mass was studied in [17].

In the current paper, we aim to extend the work of Baltagiannis [12] and Zotos [13] by performing a full analysis of the location, stability, and basins of convergence of the equilibrium points, associated to the inclusion of a radiating body in the equilateral triangle configuration of the four-body problem. The inclusion of the radiation and drag forces in the restricted four-body problem, allow us to model in a more realistic way the dynamics of a test particle in presence of an astrophysical system with, e.g., an active star. Unlike the gravitational force, radiation and drag forces are generally non-conservative, causing a loss of orbital energy such that the particles will spiral toward the source. This new effect can significantly modify the dynamics of the model and deserves a complete study in the context of the restricted four-body problem.

The present paper has the following structure: the equations of motion for a test particle are derived in section 2. In section 3, we analyze the position, linear stability, zero velocity surfaces and basins of convergence of the equilibrium points, as a function of the radiation parameter. Here, we study three combinations of mass for the primary bodies: equal masses, two equal masses, different masses. In Section 4 the complexity (unpredictability) of the basins is analyzed through the basin entropy, showing its dependence with the radiation parameter. Finally, the main conclusions of this work are drawn in section 5.

2 Equations of motion

Consider the equilateral restricted four-body problem (henceforth ERFBP). Let $m_{1},m_{2}$, and $m_{3}$ be the masses of the primaries and let $m$ be the mass of the test particle, which satisfies $m\ll m_{i}$ with $i=1,2,3$. The primaries revolve in the same plane with uniform angular velocity, and regardless of the mass distribution of the primaries, they always will lie at the vertices of an equilateral triangle. If the primary $m_{1}$ is a radiating body, the equations of motion for the test particle $m$ in an inertial frame of reference $\mathbb{R}=(X,Y)$, can be written as [18]:

[TABLE]

where $S$ denotes the solar energy flux density, $A$ the geometric cross-section of the test particle, $Q_{\text{pr}}$ the radiation pressure coefficient, $c$ the speed of light, and $sw$ the ratio of solar wind to Poynting-Robertson drag.

In Eq. (1), the first term on the right-hand side represents the influence of the gravitational forces due to the three massive bodies, while the second term includes the effect of radiation pressure and drag forces. The ratio of force due to radiation pressure compared to gravity, allow us to introduce a dimensionless radiation factor (see e.g. [8])

[TABLE]

For simplicity, in all that follows we shall use canonical units, such that the sum of the masses, as well as the distance between the primaries, the angular velocity, and the gravitational constant, are set to 1. Additionally, as we consider the non-relativistic limit of the model, the speed of light will be chosen to the value $c=1\times 10^{4}$ [19, 20], unless otherwise is specified. Taking into account the previous definitions, the equations of motion in a synodic frame of reference $\mathbb{r}=(x,y)$ read as

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

The coordinates of the primaries can be fully determined by placing the center of mass of the ERFBP at the origin of the synodic frame of reference, with $m_{1}$ lying on the $x$-axis, and taking each side of the equilateral triangle equal to unity (see Eq. (13) in Ref. [21]), i.e.,

[TABLE]

This set of equations lead to the following solutions for the coordinates of the primaries

[TABLE]

with

[TABLE]

Therefore, there are four possible cases for the location of the three massive bodies. Here, we will consider the fourth case (see Table 1).

It can be easily noted that by defining the effective potential

[TABLE]

the equations of motion (3) and (4) can be written in the compact form

[TABLE]

Due to the existence of dissipative terms, the Jacobi constant is no longer conserved and varies with time according to the following relation

[TABLE]

In absence of dissipative forces ($\beta=0$), Eqs. (7) and (8) reduce to the classical ERFBP (see e.g. [21, 12, 13]), and hence the Jacobi constant (9) is conserved

[TABLE]

3 Libration points and zero velocity surfaces

In what follows, we will determine how the radiation factor $\beta$ affects the position, stability, and basins of convergence of the libration points, when using one of the three possible combinations of mass for the primaries, $m_{1}=m_{2}=m_{3},m_{1}\neq m_{2}=m_{3},$ and $m_{1}\neq m_{2}\neq m_{3}$. To do so, let us derive some general expressions in terms of the masses of the primary bodies.

First, the location of the libration points can be determined by solving the system of equations $\dot{x}=\dot{y}=\ddot{x}=\ddot{y}=0$, which leads to the following algebraic system of equations

[TABLE]

whose solution shall depend on the parameters $\beta$ and $sw$, as well as the values of mass, where the last ones will modify also the position of the primaries $(x_{i},y_{i})$ according to (5).

In spite of the fact that the Jacobi constant is not a conserved quantity when radiation terms are included ($\beta\neq 0$), an analytic expression for the zero velocity surfaces (ZVS) can be derived as follows. From Eq. (9), we get

[TABLE]

with

[TABLE]

and setting $\dot{x}=\dot{y}=0$, time dependent terms vanish and consequently Eq. (13) reduces to

[TABLE]

Therefore, Eq. (17) defines the zero velocity surfaces of the problem, or in other words, it determines the forbidden regions of motion for a test particle. Note that for $\beta=0$, the expression for ZVS coincides with the usual expression for the classical ERFBP.

Finally, concerning the stability of the fixed points, it is a well-known fact that detailed information about the motion of a system can be obtained by linearizing the equations of motion about a fixed point, $(x^{*},y^{*})$. For this, we first introduce the following changes of variable: $\dot{x}\rightarrow\xi$, and $\dot{y}\rightarrow\eta$, such that the equations of motion (7-8) can be written as the first-order system

[TABLE]

Then, the system (18-21) can be linearized by means of the coefficient matrix

[TABLE]

with

[TABLE]

and therefore, the characteristic polynomial will be given by the quartic equation

[TABLE]

The equilibrium points are stable if all the roots of the characteristic polynomial evaluated at $(x^{*},y^{*},0,0)$ are pure imaginary roots or complex roots with negative real parts; otherwise, they are unstable.

3.1 Case 1: $m_{1}=m_{2}=m_{3}$

Let us start considering the case in which the three primary bodies have the same value of mass, i.e, $m_{1}=m_{2}=m_{3}=1/3$. According to Eqs. (5), the coordinates of the primaries $(x_{1},y_{1}),(x_{2},y_{2}),$ and $(x_{3},y_{3})$, are given respectively by

[TABLE]

In Fig. 1, we have plotted the location and evolution of the equilibrium points for $\beta$ increasing from [math] to $1$. Blue dots are the solutions for $\beta=0$, while red dots denote the coordinates of the three massive bodies111Along the paper, we shall use the customary value for the ratio of solar wind to Poynting-Robertson drag, $sw=0.35$.. It can be seen that as $\beta$ increases, $L_{1}$ and $L_{10}$ reach both the same point along the $x$-axis and completely disappear for $\beta\approx 0.687$. Also, it is observed that the libration points $L_{2}$, $L_{5}$, $L_{6}$, and $L_{9}$, move gradually toward the position of the radiating body until they disappear for $\beta\approx 0.999$. On the other hand, $L_{3}$, $L_{4}$, $L_{7}$, and $L_{8}$, exhibit a small displacement from their initial positions (see Table 2 for a detailed description).

On the subject of the stability of the fixed points, it is found that in accordance with Ref. [12], for $\beta=0$, all the equilibrium points are unstable. In the interval $\beta\in(0,1]$, the stability of the fixed points remains unaltered, since the form of the roots for all the equilibria does not change, i.e., for $L_{1}$, $L_{3}$, $L_{5}$ and $L_{6}$, the characteristic equation (23) gives place to complex eigenvalues of the form $\lambda_{1,2,3,4}=\pm a\pm ib$, while for the libration points $L_{2}$, $L_{4}$, $L_{7}$, $L_{8}$, $L_{9}$ and $L_{10}$, all roots takes the form $\lambda_{1,2}=\pm\,ib$ and $\lambda_{3,4}=\pm\,a$. So, we may conclude that if the primary bodies are equal mass, the stability of the libration points does not change with the radiation factor $\beta$.

Moreover, in Fig. 2, we show the parametric evolution of the ZVS with $\beta$, for different values of the Jacobi constant. Here, it can be noted that the regions of allowed motion are substantially modified for larger values of $\beta$. For example, when $C=2.92$ (blue contour) and for $\beta=0.125$, the forbidden region is made up of three small islands located just around the libration points $L_{3}$, $L_{5}$, and $L_{6}$. Then, the test particle will be free to move almost without any restriction around the primary bodies. However, as $\beta$ increases, these small islands get bigger, firstly, forming a small connection near $m_{1}$ until completely surround it. For $\beta\rightarrow 1$, the forbidden region in left-hand-side joins the one surrounding $m_{1}$, until the allowed regions of motion around $m_{2}$ and $m_{3}$, become almost isolated. This effect can be explained if we consider that, $\beta\approx 1$, implies that the force due to radiation pressure is comparable to the force of gravity, being able to push the test particle away from $m_{1}$.

We finish this subsection by discussing the so-called basins of convergence and their evolution with the radiation factor $\beta$. As pointed out in section 1, by basin of convergence, we refer to the set of points that after successive iterations converge to a specific fixed point. The numerical procedure for obtaining such basins is the multivariate Newton-Raphson method, which can be defined by the map

[TABLE]

where ${\bf{x}}=(x,y)$, $f({\bf{x_{n}}})$ represents the system of equations (11-12), and $J^{-1}$ is the inverse Jacobian matrix. The initial guesses are defined within the region enclosed by $x\in\left[-8,5;\,8,5\right]$ and $y\in\left[-8,5;\,8,5\right]$, the step size is equal to $10^{-2}$, and tolerance of the order $10^{-12}$. Each initial condition leading to a given fixed point is plotted using the following color code: $L_{1}\rightarrow$ dark green, $L_{2}\rightarrow$ dark grey, $L_{3}\rightarrow$ black, $L_{4}\rightarrow$ orange, $L_{5}\rightarrow$ purple, $L_{6}\rightarrow$ red wine, $L_{7}\rightarrow$ light blue, $L_{8}\rightarrow$ blue, $L_{9}\rightarrow$ green and $L_{10}\rightarrow$ yellow; while the libration points are denoted by cyan dots.

In Fig. 3, we plot the basins of convergence of the equilibrium points, for nine different values of the radiation factor $\beta$. Here, it can be easily observed that some of the equilibrium points disappear as $\beta$ increases, or in other words, the number of colors for $\beta=1$ is reduced compared to the case $\beta=0$. It should be noted that as soon as $\beta>0$ the $2\pi/3$ symmetry is broken and replaced by a reflection symmetry with respect to the $x$-axis. Additionally, the complexity of the basins increases for higher values of the radiation parameter, despite the fact that the number of libration points is much smaller than in the case $\beta=0$. The last statement can be clearly envisioned in the last panel of Fig. 3, where the right-hand side region becomes very noisy, suggesting an increasing unpredictability of the basins.

3.2 Case 2: $m_{1}\neq m_{2}=m_{3}$

As a second case, we have assumed that $m_{2}=m_{3}={\mathfrak{m}}$. Accordingly, the mass parameter of the radiating body will be given by $m_{1}=1-2{\mathfrak{m}}$, and the coordinates of the primaries $(x_{1},y_{1}),(x_{2},y_{2}),$ and $(x_{3},y_{3})$, read as

[TABLE]

The last result implies that not only the position of the primaries depends on the parameter ${\mathfrak{m}}$, but the existence and location of libration points will also depend on this parameter and on $\beta$, according to Eqs. (11, 12). In Table 3, we present how the total number of equilibrium points varies on different intervals of the radiation factor and for different values of ${\mathfrak{m}}$. It can be observed that for larger values of $\beta$ the total number of fixed points is reduced, where the final number of surviving points increases for large values of ${\mathfrak{m}}$.

In Figure 4, we present a graphic representation of the evolution of equilibrium points with the variation of the radiation factor $\beta$, for different values of $m_{2}=m_{3}={\mathfrak{m}}$. This figure can be analyzed as follows: (i) in the case $m_{1}\gg{\mathfrak{m}}$, it can be noted that the collinear points $L_{3}$ and $L_{4}$ do not exist. For larger values of $\beta$, the equilibria $L_{1}$ and $L_{10}$ approximate each other until they join near the $x$-axis, while the points $L_{2}$, $L_{5}$, $L_{6}$, and $L_{9}$ move toward the radiating body until they disappear for $\beta=1$. (ii) As expected, when $m_{1}\approx{\mathfrak{m}}$ the results are exactly like the ones described in the previous subsection. (iii) For $m_{1}\ll{\mathfrak{m}}\approx 0.5$, the non-collinear points $L_{9}$ and $L_{10}$ do not exist. As $\beta$ increases, $L_{3}$, $L_{4}$, $L_{7}$, and $L_{8}$ exhibit very small displacements with respect to their initial positions, while the libration points $L_{1}$, $L_{2}$, $L_{5}$, and $L_{6}$ move abruptly toward the radiating body, until they disappear for $\beta=1$.

Concerning the stability of the fixed points, we start considering the case $\beta=0$. Here, we found that if $\mathfrak{m}\in(0,0.0027]$ the libration points $L_{3}$, $L_{5}$, and $L_{6}$ are stable, for $\mathfrak{m}\in[0.0027,0.0188]$, only $L_{5}$, and $L_{6}$ remain stable, but for $\mathfrak{m}>0.0188$ all the equilibria are unstable. All those results agree with the ones given in Ref. [12]. Once we start varying $\beta$ in steps of $\Delta\beta=1\times 10^{-3}$, we get: (i) for $\beta=0.001$ and $\mathfrak{m}\in(0,0.0024]$ the libration points $L_{3}$, $L_{5}$ and $L_{6}$ are stable; for $\mathfrak{m}\in[0.0025,0.0185]$ only $L_{5}$ and $L_{6}$ are stable, while for $\mathfrak{m}>0.0185$ all the equilibria are unstable. (ii) For $\beta=0.002$ and $\mathfrak{m}\in(0,0.0011]$ the fixed points $L_{3}$, $L_{5}$ and $L_{6}$ are stable; for $\mathfrak{m}\in[0.0012,0.0174]$ only $L_{5}$ and $L_{6}$ are stable, but for $\mathfrak{m}>0.0174$ all the fixed points are unstable. (iii) For $\beta=0.003$ and $\mathfrak{m}\in(0,0.0147]$ only $L_{5}$ and $L_{6}$ are stable, while for $\mathfrak{m}>0.0147$ all the equilibria are unstable. (iv) For $\beta=0.004$ and $\mathfrak{m}\in(0,0.0093]$ $L_{5}$ and $L_{6}$ are stable, but for $\mathfrak{m}>0.0093$ all the fixed points are unstable. (v) Finally, for $\beta>0.004$, all fixed points are unstable regardless the value of $\mathfrak{m}$.

From the previous results, we can infer that the stability of the libration points remains unaltered, only for values of radiation factor of the order $10^{-3}$, however, for larger values of $\beta$, all the equilibria become unstable. At this point, it should be noted that the radiation factor of the solar system is of the order $10^{-1}$, which implies that, in a realistic system, such factor should be able to destroy the stability of the fixed points.

On the other hand, in Figure 5, we present the ZVS by using different values of the Jacobi constant. The tendency in this figure is very similar to the one observed in Fig. 2. For example, in the case $C=2.64$ and for small values of $\beta$, say $\beta<0.25$, the test particle is free to move in the whole configuration space, however, as $\beta$ increases the libration points $L_{5}$ and $L_{6}$ become isolated points. Such islands start growing in size until the ZVS surrounds completely $m_{1}$, in particular, for $\beta=0.875$ the motion of the test particle is limited to a small region around $m_{1}$, or to a larger zone that allows collisions with $m_{2}$ and $m_{3}$, but imposes a barrier in the vicinity of $m_{1}$. This behavior is consistent with the idea that, for $\beta\rightarrow 1$, the drag forces are strong enough to push away the test particle from $m_{1}$.

In Figure 6 we use the same color code of Fig. 3. Here, it can be noted that in the case $m_{1}=0.5$, $m_{2}=m_{3}=0.25$, the basins of attraction exhibit a very complex pattern symmetric in shape, but asymmetric in color. Due to the fact that the number of fixed points grows for intermediate values of $\beta$ ($\beta\in[0.25,0.75]$) and decreases for larger values of the radiation factor ($\beta>0.875$), the highly fractal structures are always present. Related to the extent of the basins of convergence, it is seen that for $\beta\in[0,0.375]$ the area corresponding to libration point $L_{9}$ and $L_{10}$ predominates, while for $\beta\in[0.5,0.75]$ the area of the set of points converging to the fixed point $L_{1}$ is bigger than the one for $L_{9}$ and $L_{10}$. Special attention deserves the case $\beta=0.625$, where the basins seem stretched, here there are three prevailing colors, the ones corresponding to $L_{1},L_{9}$ and $L_{10}$, i.e., the fixed points that are closer to the origin of coordinates. Such a structure allows us to prognosticate a reduction of the fractality of the basin.

Finally, for $\beta>0.875$, the area corresponding to the basin for $L_{9}$ increases, while for $\beta=1$ it seems to be dominated by $L_{3}$, however, in the inset of the lower right panel, it is observed that the black region is, in fact, a very noisy region composed of initial conditions tending to all the libration points. From the previous description, we may infer that in the case of two equal masses, the determinations of the final state via the root finding algorithm of Newton-Raphson is highly sensitive, mainly, for lower and higher values of $\beta$.

3.3 Case 3: $m_{1}\neq m_{2}\neq m_{3}$

As a last scenario, we consider the equilateral triangle configuration Sun-Jupiter-Trojan Asteroid, in addition with the P-R and solar wind drag forces. In the normalized units, the masses of the primaries are given by $m_{1}=m_{S}=0.999046321943$, $m_{2}=m_{J}=0.000953678050$ and $m_{3}=m_{A}=6.99996\times 10^{-12}$, where $m_{3}$ corresponds to the mass of 624 Hektor, an actual asteroid of the Trojan group222Unlike the previous cases, all calculations of this subsection were performed in quadruple precision and the speed of light was set as in Baltagiannis [22].

According with Eqs. (5), the coordinates of the primary bodies are now

[TABLE]

In agreement with the seminal paper by Baltagiannis [22], in absence of radiation pressure and drag forces ($\beta=0$), the Sun-Jupiter-Trojan Asteroid system admits eight non-collinear equilibrium points, three of them ($L_{6}$, $L_{7}$, and $L_{8}$) linearly stable. The adjective non-collinear is due to the fact that $L_{1}$, $L_{2}$, and $L_{3}$, do not lie exactly on the $x$-axis. As $\beta$ increases, the equilibria $L_{4}$ and $L_{6}$, reach both the same point and completely disappear for $0.0029<\beta<0.0030$. The same applies for $L_{3}$ and $L_{7}$ for $0.9887<\beta<0.9888$, and for $L_{1}$ and $L_{8}$ for $0.9940<\beta<0.9941$ (we refer the reader to Table 4 for details). On the other hand, $L_{2}$ and $L_{5}$ hardly move, approaching each one to its nearest primary (See Fig. 7).

On the stability subject, for $\beta\in(0,0.0029]$ the equilibria $L_{6}$, $L_{7}$, and $L_{8}$ are linearly stable, however, as $\beta$ grows the stability of the fixed points also changes, i.e., for $\beta\in$ $[0.0030,0.0104]$, only $L_{7}$ and $L_{8}$ remain stable, while for $\beta>0.0104$, all the fixed points become unstable. Here, it can be noted that $\beta\approx 10^{-2}$ corresponds to the largest value of the radiation factor (of the three considered cases) in which we still having stable fixed points. So, taking into account that $\beta\approx 10^{-1}$ in the Solar system, we may conclude that in the restricted four-body problem, Sun-Jupiter-Trojan Asteroid-Spacecraft (test particle), all the libration points must be unstable.

In Fig. 8 we present the parametric evolution of the ZVS with $\beta$, for different values of the Jacobi constant. In absence of the P-R and drag forces, the forbidden regions of the Sun-Jupiter-Trojan Asteroid system are made of horseshoe-shaped contours, very similar to the ones observed in the Sun-Jupiter system in the circular restricted three-body problem (see e.g. [23]). The circular horseshoe is open on $m_{2}$, such that this opening will be greater at lower values of $C$. Keeping a fixed value of $C$, as $\beta$ increases, the mentioned horseshoe closes on it selves forming a donut-shaped contour in which the allowed motion region becomes smaller and smaller around the radiating body. For larger values of $\beta$ ($\beta\rightarrow 1$), the circular region around $m_{1}$ disappears, and the forbidden region of motion becomes a circle containing the Lagrangian points $L_{2,5}$ and the three masses. Here, we would like to highlight that the location of $m_{3}$ almost coincides with the position of $L_{5}$, for this reason, the red dot of $m_{3}$ is not seen in Fig. 8.

On the other hand, in Fig. 9, we plot the basins of convergence of the equilibrium points for the Sun-Jupiter-Trojan Asteroid system. The color code used in this case is the same as in figure 3. For $\beta=0$, the frame exhibits the eight colors of the corresponding libration points, albeit the most noticeable basins are the ones corresponding to $L_{1},L_{3},L_{5},L_{6},L_{7}$, and $L_{8}$. In general, it is observed that along the negative $x-$axis the basin associated to $L_{1}$ predominates, while along the positive $x-$axis the one associated to $L_{3}$ does. In the remaining frames, the number of equilibria decreases to 6 ($L_{1}$, $L_{2}$, $L_{3}$, $L_{5}$, $L_{7}$, and $L_{8}$), then 4 ($L_{1}$, $L_{2}$, $L_{5}$, and $L_{8}$) and finally 2 ($L_{2}$ and $L_{5}$). However, due to the location of the equilibria $L_{5}$ (superimposed with $m_{3}$), no initial condition of our mesh converges to this point and therefore, the purple basin neither appears.

In accordance with cases 3.1 and 3.2, the complexity of the basins increases for larger values of $\beta$, this is so because, for $\beta\rightarrow 1$, $L_{3}$ (whose initial position starts alongside the x-axis) moves downward, while $L_{1}$ moves upward, $L_{7}$ moves in counterclockwise and $L_{8}$ moves toward the radiating body $m_{1}$. Such displacements shall break the symmetry of the basins of convergence and hence increasing the intricacy of them (see e.g. the lower-mid panel in Fig. 9). Lastly, when $\beta\rightarrow 1$ we observe a practically white configuration space, this color is assigned to non-converging initial conditions, however, the inset of this panel shows a region composed by scattered dark gray dots. This results can be understood by considering that for $\beta=1$ there exist only 2 fixed points, one of them corresponds to $L_{5}$ which, due to its superposition with the primary $m_{3}$, do not allow for convergence of any initial condition of our mesh after $10^{5}$ iterations.

4 Basin entropy

So far, the analysis of the basins of convergence has been purely qualitative, in order to make quantitative our study, in this section, we shall use a quantity recently introduced in [24] which allow us to measure the uncertainty (complexity) of a given basin. The new dynamical quantity is termed basin entropy and provides a useful method to explore the differences observed in our system when the radiation and mass parameters are varied.

For the sake of completeness, we will briefly describe the idea behind the calculation of the basins’ entropy. Assuming that the phase space contains $N_{A}$ different final states (or colors in our case), we divide the space of states into a grid of $N$ square cells, such that each one of these cells contains at least one of the $N_{A}$ states. Defining $p_{i,j}$ as the probability to detect a state $j$ in the $i-$th cell and by applying the Gibbs entropy definition to that set, the entropy for the $i-$th cell can be expressed as

[TABLE]

Therefore, the basin entropy for the total number of cells $N$ in the basin is calculated as an average, i.e.,

[TABLE]

Strictly speaking, the average of this quantity must depend on the number of considered cells, such that for a larger value of $N$ the result for the basin entropy $S_{b}$ should be more precise. However, a bigger $N$ can be reached only with a smaller size of the square cells, which also have a minimum size to contain at least one of the $N_{A}$ states. To solve this issue, we follow the procedure outlined in [25], in which the square cells are randomly picked in the space of states through a Monte Carlo procedure, allowing us to increase the number of cells $N$ as necessary. In our particular problem, we find that the final value for the basin entropy keeps constant for a number of cells larger than $3\times 10^{5}$, hence, in the three cases, we used $N=3.5\times 10^{5}$ cells.

In Fig. 10, we present the parametric evolution of the basin entropy as a function of the radiation parameter $\beta$, for the considered cases. Our results indicate that, in case 1, the basin entropy increases almost monotonically with the radiation parameter $\beta$. A very similar behavior is observed for case 3 in the interval $\beta\in[0.025,0.975]$, yet for $\beta\in[0,0.025]$ and $\beta\in[0.975,1]$, the basin entropy decreases abruptly. Interestingly, the opposite tendency is observed for case 2, where the basin entropy decreases with the radiation parameter $\beta$ in the interval $\beta\in[0,0.825]$, but increases in the interval $\beta\in[0.825,1]$ from $S_{b}\approx 0.1$ up to a value of $S_{b}\approx 1.14$. In general, our findings suggest that in cases 1 and 2, the unpredictability associated to the NR basins of convergence for the triangular restricted four-body problem with an extreme radiating body ($\beta=1$) is larger in comparison with the non-radiating case. Nevertheless, in the case of Sun-Jupiter-Trojan Asteroid-Spacecraft system (case 3), the exact opposite situation occurs. These results can be explained if we consider that in case 3, the final number of libration points is one-fourth of the initial points, i.e., the number of final states $N_{A}$ is considerably diminished, and hence, the value of Eq. (26) is also significantly modified.

5 Discussion and conclusions

In the present paper, we numerically investigated the location, stability, and basins of convergence of the equilibrium points, in the equilateral triangle configuration of the four-body problem with a radiating body. Specifically, we showed how the radiation parameter influences the dynamics of the system, in three different combinations of mass for the primaries: equal masses, two equal masses, and three different masses.

Following the tendency in this field, we used a multivariate Newton-Raphson method in order to calculate the corresponding basins of convergence. Such basins give a global view of the set of initial conditions that, after an iterative process, show a tendency to a certain fixed point. In each of the considered cases, we examine the influence of the radiation parameter on the Newton-Raphson basins of convergence. Aiming to quantify the uncertainty (complexity) of the basins, we calculated the basin entropy, monitoring its variation with the radiation parameter.

The main conclusions of this work can be summarized as follows:

In the three considered cases, the total number of libration points decreases as $\beta\rightarrow 1$. 2. 2.

The stability analysis suggests that, in the case of equal masses and for $\beta\in[0,1]$, the equilibrium points of the system are always linearly unstable. In the rest of the cases, the libration points are unstable if $\beta>0.01$. 3. 3.

Taking into account the radiation factor for the solar system, we conclude that the radiation pressure and drag forces should be able to destroy the stability of the fixed points in the restricted four-body problem composed by Sun, Jupiter, Trojan asteroid and a test (dust) particle. 4. 4.

Only in case 3 (Sun-Jupiter-Trojan Asteroid - Spacecraft system), we detected the existence of a considerable amount of non-converging initial conditions after 50000 iterations. This is due to the fact that there exists only one effective fixed point that can be reached after the iterative process. 5. 5.

The lowest value of the basin entropy was found in the case of two equal masses near $\beta=0.825$, while on the other hand the highest value of $S_{b}$ was measured in the same case for $\beta=1$. 6. 6.

It is found that in the cases of two or three equal masses, the unpredictability associated with the basin of convergence is larger for the extreme radiating case in comparison with the non-radiating one.

Acknowledgments

FLD acknowledge financial support from Universidad de los Llanos, Grant No. CDP 2478. FLD and GAG gratefully acknowledge the financial support provided by COLCIENCIAS (Colombia), Grants Nos. 8840 and 8863.