## ORIGINAL RESEARCH article

Front. Astron. Space Sci., 27 April 2022
Sec. Fundamental Astronomy
Volume 9 - 2022 | https://doi.org/10.3389/fspas.2022.877459

# Existence and Stability of Equilibrium Points in the R3BP With Triaxial-Radiating Primaries and an Oblate Massless Body Under the Effect of the Circumbinary Disc A. E. Vincent1 A. E. Perdiou2* E. A. Perdios3
• 1Department of Mathematics, School of Basic Sciences, Nigeria Maritime University, Okerenkoko, Nigeria
• 2Department of Civil Engineering, University of Patras, Patras, Greece
• 3Department of Electrical and Computer Engineering, University of Patras, Patras, Greece

The aim of this article is to study the existence, location, and stability of equilibrium points in a generalized restricted three-body problem (R3BP) that consists of an oblate infinitesimal body when the primaries are radiating sources with triaxiality of the two stars surrounded by a belt (circumbinary disc). The existence, number, location, and stability of the collinear and triangular Lagrangian equilibrium points of the problem depend on the mass parameter and the perturbing forces involved in the equations of motion. We find numerically that four additional collinear equilibrium points Lni, i = 1, 2, 3, 4, exist, in addition to the three Eulerian points Li, i = 1, 2, 3, of the classical case, making up a total of up to seven collinear points. Ln1 and Ln2 result due to the potential from the belt, while Ln3 and Ln4 arise from the effect of triaxiality. The positions of the equilibrium points are affected by the presence of perturbations, since they are deviated from the classical R3BP on the x-axis and out of the x-axis, respectively. The stability of the equilibrium points, for a particular set of the parameters, is analyzed, and it is concluded that all the collinear points are unstable except Ln1, which is always linearly stable. The range of stability of the Lagrangian points L4,5 is determined analytically and found that being stable for 0 < μ < μcrit and unstable for μcritμ ≤ 1/2, where μcrit is the critical mass ratio which depends on the combined effects of the perturbing forces. It is noticed that the critical mass ratio decreases with the increase in the values of the radiation pressure, triaxiality, and oblate infinitesimal body; however, it increases with the increase in the value of mass of the disc. All three of the former and the latter one possess destabilizing and stabilizing behavior, respectively. The net effect is that the size of the region of stability that decreases when the value of these parameters increases. In our model, the binary HD155876 system is used, and it is found that there exists one stable collinear equilibrium point viz. Ln1.

## Introduction

There is no doubt that the restricted three-body problem (R3BP) is certainly the most known and studied problem in celestial mechanics as well as in dynamical astronomy since the age of Euler and Lagrange and still does, due to the fact that the study of this problem in its many variants has had important implications in several scientific fields, including galactic dynamics, chaos theory, molecular physics, and lunar theory, among others (see, for example, Szebehely, 1967; Prosmiti et al., 1996; Vrahatis et al., 2001; Contopoulos, 2002; Marsden and Ross, 2005; Douskos et al., 2007; Sano, 2007; Musielak and Quarles, 2017; Voyatzis and Antoniadou, 2018; Zotos et al., 2020). The circular R3BP models the motion of an infinitesimal mass m3 which moves under the gravitational attraction of two massive bodies m1 and m2, which are called primaries. We refer to these bodies as the bigger and smaller primaries, respectively. The primaries revolve in circular orbits around their common center of mass, and their motion is not affected by the massless body. The classical circular R3BP admits five equilibrium points denoted as Li,  i = 1, 2, ..., 5. Three of such points, Li, i = 1, 2, 3, lie on the line connecting the two massive bodies m1 and m2, and they are called collinear equilibrium points which are unstable, while the other two points form an equilateral triangle with the two massive bodies, called the triangular equilibrium points L4 and L5, and are stable for 0 < μ < μc = 0.03852…, where μ is the mass parameter, and μc is the critical mass parameter (see, for e.g., Valtonen and Karttunen, 2006). These points and their stability are among the most characteristic features of any dynamical system not only due to their immediate applications in real systems but also due to their theoretical importance since the nature of these points characterizes the behavior of nearby orbits (Markellos et al., 1996; Hou and Liu, 2008; Capdevila and Howell, 2018).

During the past, several variants of this classical problem have been proposed in order to make it more applicable to real systems of dynamical astronomy. These modifications include the consideration of one or both primaries being sources of radiation pressure and/or oblateness of the third body and/or triaxial rigid bodies or even some more simplified versions such as the well-known Hill’s problem (see, for e.g., Voyatzis et al., 2012; Singh, 2013; Elshaboury et al., 2016; Amuda et al., 2021; Gao and Wang, 2020; Kalantonis, 2020; Singh et al., 2018; Saeed and Zotos, 2021), among others. Another interesting modification of the R3BP is the Chermnykh’s problem which can be considered a generalization of the Euler’s problem of two fixed gravitational centers and the R3BP when the massless body moves in the orbital plane of a dumbbell which is rotating with constant angular velocity around the center of mass (Chermnykh, 1987; Goździewski and Maciejewski, 1998; Perdios and Ragos, 2004; Papadakis, 2005; Perdiou et al., 2013). In Yeh and Jiang (2006) and Jiang and Yeh (2006), the authors studied the Chermnykh-like problem analytically and numerically, in which an additional gravitational potential from the belt is included in the model. Other works that took into account the gravitational potential from the belt/disc under different assumptions have been presented and studied by Kishor and Kushvah (2013), Abouelmagd et al. (2014), Singh and Leke (2014), Perdios et al. (2015), Singh and Amuda (2019), Kalantonis et al. (2021), Leke and Singh (2021), among others.

In the case of a triaxial rigid primary body, Saeed and Zotos (2021) revealed the way numerically in which the triaxiality parameters affect the existence and linear stability of the libration points and explored the entire range of values of the triaxiality parameters that control the shape of the primary body. It was observed that the number of collinear and noncollinear equilibrium points depends on the triaxiality parameters of the primary. They also showed that the equilibria are always linearly unstable for all possible values of the parameters of the primary body. Recently, Gyegwe et al. (2022) have studied the positions and stability of the triangular equilibrium points in the framework of the R3BP when the primaries are triaxial-radiating rigid bodies, while the third body of infinitesimal mass is an oblate spheroid. They found that the parameters of the problem play a significant role on the regions of stability and in particular that the comprehensive effects of the perturbations have destabilizing tendencies.

In the present work, we aim to extend the work of Gyegwe et al. (2022) by considering not only the radiation pressure and asphericity of the primaries but also the gravitational potential of the belt. Moreover, taking into account these perturbing forces, we investigate the existence, location, and stability of both the collinear and triangular equilibrium points of this system. New equilibrium points, zero-velocity surfaces, and allowed regions of motion make a qualitative difference to the dynamical features of the model. In our study, the triaxiality factors of the primaries are of small quantities, and therefore we choose arbitrary values that are far less than one (see, for e.g., Singh, 2013; Elshaboury et al., 2016, and references therein). Our goal in this article was to study important aspects of the dynamics of the restricted problem with the involved parameters from both numerical and analytical points of view. Several results concerning the equilibrium points and stability are proved analytically, and the study is completed numerically in some cases.

The content of this article is organized as follows: in Equations and Integral of Motion section, we define the mathematical model and the governing equations of motion that we use for our study. In Equilibrium Points and Zero-Velocity Surfaces section, we determine the existence and locations of the equilibrium points as well as the geometric structure of the zero-velocity surfaces and curves. The linear stability of the equilibria is analyzed in Stability of Collinear and Triangular Points section, while in Numerical Application section, numerical and graphical simulation is made by using the physical data of the binary HD155876 system for various values of the parameters under consideration. Finally, Discussion and Conclusion section summarizes the discussion and conclusion of our study.

## Equations and Integral of Motion

We adopt the usual barycentric, rotating, and dimensionless coordinate system Oxyz, in which its origin is located at the center of mass of the two primary bodies with masses m1 = 1 − μ and m2 = μ, μ = m2/(m1+m2) ∈ (0, 1/2], have fixed positions at the Ox-axis. Then, the coordinates of the bigger and smaller primaries are (−μ, 0) and (1 − μ, 0), respectively, while the coordinates of the third body of negligible mass are (x, y, z). Following Gyegwe et al. (2022) and Singh and Leke (2014), the equations of motion of the massless body under the influence of its shape, forces of radiation pressure, and triaxiality coupled with gravitational potential from a belt around the stars are:

$x¨−2ny˙=∂U∂x=Ux; y¨+2nx˙=∂U∂y=Uy(1)$

where

$Ux =n2x−(1−μ)(x+μ)q1r13−μ(x+μ−1)q2r23−3(1−μ)(2σ1−σ2)(x+μ)q12r15−3μ(2σ1′−σ2′)(x+μ−1)q22r25+15(1−μ)(σ1−σ2)(x+μ)y2q12r17+15μ(σ1′−σ2′)(x+μ−1)y2q22r27−3(1−μ)(x+μ)A32r15−3μ(x+μ−1)A32r25−Mbx(r2+T2)32,Uy =n2y−(1−μ)yq1r13−μyq2r23−3(1−μ)(σ1−σ2)yq1r15−3μ(σ1′−σ2′)yq2r25−3(1−μ)(2σ1−σ2)yq12r15−3μ(2σ1′−σ2′)yq22r25+15(1−μ)(σ1−σ2)y3q12r17+15μ(σ1′−σ2′)y3q22r27−3(1−μ)yA32r15−3μyA32r25−Mby(r2+T2)32,(2)$

are the partial derivatives of the potential function,

$U =n22(x2+y2)+(1−μ)q1r1+μq2r2+(1−μ)(2σ1−σ2)q12r13+μ(2σ1′−σ2′)q22r23−3(1−μ)(σ1−σ2)q1y22r15−3μ(σ1′−σ2′)q2y22r25+(1−μ)A32r13+μA32r23+Mb(r2+T2)12,(3)$

with

$n=1+32(2σ1−σ2)+32(2σ1′−σ2′)+2Mbrc(rc2+T2)32,0<σi,σi′,Mb,T<<1,i=1,2,r1=[(x+μ)2+y2]12,r2=[(x+μ−1)2+y2]12,r=(x2+y2)12,rc=(1−μ+μ2)12.(4)$

Here, we denote by ri, i = 1, 2, the distances of the third body from the bigger and smaller primaries, respectively, by 0 < q i= 1 − δi ≤ 1, i = 1,   2, the radiation factor of each primary, where δ1 and δ2 are the ratios of the radiation force Fri to the gravitational force Fgi which results from the gravitation due to the primary bodies m1 and m2, respectively, by A3, the oblateness coefficient of the test body defined by the formula A3 = (AE2AP2)/5R2 << 1, where AE and AP are the equatorial and polar radii of the said third body, respectively, and R is the dimensionless distance between the primaries. The triaxiality of the bigger and smaller primaries is defined with the help of the parameters $σi,σi′<<1,i=1,2$, where

$σ1=h2−f25R2, σ2=b2−f25R2, σ1′=h′2−f′25R2, σ2′=b′2−f′25R2,$

with h, b, f as lengths of the semiaxis of the bigger primary and h, b, f as those of the smaller primary. Also, 0 < Mb << 1 is the total mass of circular cluster of material points, and r is the radial distance of the test body while rc is that in the classical case. T = a + b is the density profile of the accumulated materials chosen to the value T = 0.01 (Jiang and Yeh, 2006), while a and b are the flatness and core parameters, respectively, while n is the mean motion of the primaries. We assume that the admissible region for A3 in this study is [0, 0.1], whereas we restrict the triaxiality coefficients $σi,σi′,i=1,2$, to the interval [0, 0.12], while the radiation factors q1 and q2 vary in the interval qi ∈ (0, 1], i = 1, 2, and the region of Mb is [0, 0.15]. Such a region should be enough to show the influence of the parameter on the equilibria in which the upper boundary can be found in that given in Yeh and Jiang (2006). The Jacobian-type integral of this problem has the following expression:

$C=2U−(x˙2+y˙2),(5)$

where C is the Jacobian constant. Since the square of the particle velocity must be positive or at least zero, relation (Eq. 5) may provide the following equation:

$C=2U(x,y),(6)$

where its plots are known as zero-velocity curves (ZVC), and they separate the (x, y) plane in regions in which the particle motion is allowed and in regions in which this is not allowed to take place. If we consider a third dimension that counts the Jacobian constant C, then we obtain the so-called zero-velocity surface (ZVS). Here, we note that the aforementioned model reduces to Gyegwe et al. (2022) when Mb = 0, while for Mb = 0 and A3 = 0 some of the results can be found in Singh (2013).

## Equilibrium Points and Zero-Velocity Surfaces

To find the positions of the equilibrium points of Eq. 1, we have to set $x˙=y˙=x¨=y¨=0$ in the equations of motion. Precisely, we have the following results:

### Collinear Equilibrium Points (y = 0)

The collinear (or Eulerian) equilibrium points lie on the Ox-axis and in this case y = 0 which can be determined by solving the nonlinear equation numerically for x:

$g(x,0)=n2x−(1−μ)(x+μ)q1r13−μ(x+μ−1)q2r23−3(1−μ)(2σ1−σ2)(x+μ)q12r15−3μ(2σ1′−σ2′)(x+μ−1)q22r25−3(1−μ)(x+μ)A32r15−3μ(x+μ−1)A32r25−Mbx(r2+T2)32=0,(7)$

with r1 = |x + μ| and r2 = |x + μ − 1| for certain values of the parameters of the problem A3, Mb, qi, σi, $σi′,i=1,2,$ and μ. It has been known that there may be five collinear equilibrium points for the case of mass of disc (Yeh and Jiang, 2006) or triaxiality parameters of the massive primary body (Saeed and Zotos, 2021). In this article, such investigations will be extended to the case of both primaries with radiation, triaxiality, mass of disc, and oblate infinitesimal body.

In Figure 1, we present the exact locations of the collinear equilibrium points of the problem, for two different values of Mb, which have been obtained by solving numerically Eq. 7. In particular, Figure 1A illustrates the position of the three classical collinear points Li, i = 1, 2, 3, for the case where Mb = 0.00005. For Mb = 0.05, we observe in Figure 1B that the problem has two more collinear points located between the primaries to be referred to as Lni, i = 1, 2 (the first either for x > 0 or x < 0 depending on the relative values of q1, q2, while the second one is for x < 0 regardless of the values of the parameters q1, q2). In both cases, we have fixed the remaining parameters to μ = 0.455, σ1 = 4 × 10−5, σ2 = 3 × 10−5, $σ1′= 2×10-5$, $σ2′= 10-5$, A3 = 0.02, T = 0.01, q2 = 0.99, and q1 = 0.8. Figure 1C shows the case of seven collinear points (including the two collinear points Ln3 and Ln4 on either side of the bigger primary body m1 between Ln2 and L3), similar to panel (B) but for σ1 = 0.01, σ2 = 0.03, σ3 = 0.05, σ4 = 0.07, and A3 = 0.001. From these figures (Figures 1A,B), we can observe that for some values of the system parameters, the additional equilibria do not in general exist even in the presence of the circular cluster of material points, indicating that such equilibria only exist depending on the mass parameter and mass of the enclosed material points. Additionally, the mass parameter does not affect significantly the existence of the new collinear equilibria Ln3 and Ln4 but due to the triaxiality of the primaries. As a summary of Figure 1C, the corresponding locations of the seven libration points using the special values of the parameters of the problem mentioned previously can be identified as follows:

$L1(0.139966,0), L2(1.19480,0), L3(−1.08861,0),Ln1(−4.05331×10−6,0), Ln2(−0.133979,0), Ln3(−0.339265,0), Ln4(−0.569789,0).$

FIGURE 1

FIGURE 1. (A) Positions of the three collinear equilibrium points (red dots) Li i = 1, 2, 3, for Mb = 0.00005. (B) Five collinear points Li, i = 1, 2, 3, Lni, i = 1, 2, for Mb = 0.05. In both frames, we have fixed μ = 0.455, σ1 = 4 × 10−5, σ2 = 3 × 10−5, σ3 = 2 × 10−5, σ4 = 10−5, A3 = 0.02, T = 0.01, q2 = 0.99, and q1 = 0.8. (C) similar to panel (B), but for the case of seven collinear points in which σ1 = 0.01, σ2 = 0.03, σ3 = 0.05, σ4 = 0.07, and A3 = 0.001. The positions of the primary bodies, mi, i=1,2, are denoted by black dots.

We note that the triaxiality of the primaries has a high degree of complexity in a dynamical system, and any subtle changes in different parameters may substantially change the third body’s dynamic behavior (see, for e.g., Saeed and Zotos, 2021). The stability of the system equilibria is a significant problem for the understanding of the dynamical behaviors around these points. Based on the results given by Saeed and Zotos (2021), it can be predicted that the equilibria Ln3 and Ln4 are unstable. As a preliminary study, the discussions about these new equilibria resulting from triaxiality of the primaries will be ended by only examining the existence and stability of these points. Further studies can be made if needed, especially when these parameters are “pushed” at their limits.

### Triangular Equilibrium Points (y ≠ 0)

For zero velocity and acceleration components in the equations of motion (Eq. 1) with y ≠ 0, the positions of the equilibria are the solutions of the nonlinear algebraic system Ux = Uy = 0, which yields:

$n2x−(1−μ)(x+μ)q1r13−μq2(x+μ−1)r23−3(1−μ)(2σ1−σ2)(x+μ)q12r15−3μq2(2σ1′−σ2′)(x+μ−1)2r25+15(1−μ)(σ1−σ2)y2q1(x+μ)2r17+15μ(σ1′−σ2′)y2q2(x+μ−1)2r27−3(1−μ)(x+μ)A32r15−3μ(x+μ−1)A32r25−Mbx(r2+T2)32=0,(8)$
$n2y−(1−μ)q1yr13−μq2yr23−3(1−μ)(2σ1−σ2)q1y2r15−3(σ1−σ2)(1−μ)q1yr15−3μ(2σ1′−σ2′)q2y2r25 −3(σ1′−σ2′)μq2yr25+15(1−μ)(σ1−σ2)q1y32r17 +15μ(σ1′−σ2′)q2y32r27−3(1−μ)A3y2r15 −3μyA32r25−Mby(r2+T2)32=0.(9)$

The real solutions (x, y) of Eqs 8, 9 for those parameter values for which they exist provide the coordinates of the equilibrium points. Notice that

$Ux(x,−y;μ,q1,q2,σ1,σ2,σ1′,σ2′,A3,Mb,T)=Ux(x,y;μ,q1,q2;σ1,σ2,σ1′,σ2′,A3,Mb,T)$

and

$Uy(x,−y;μ,q1,q2,σ1,σ2,σ1′,σ2′,A3,Mb,T)=−Uy(x,y;μ,q1,q2,σ1,σ2,σ1′,σ2′,A3,Mb,T),$

and thus for each solution (x, y) of the equations with y > 0, (x, − y) is also a solution. So, the triangular equilibrium points come in pairs. Using the perturbation technique as that in Gyegwe et al. (2022), Eqs 8, 9 yield:

$x=12−μ−(1−q1)3+(1−q2)3−12(1μ−14)σ1+12(1μ+34)σ2+12(μ1−μ+34)σ1′−12(μ1−μ+74)σ2′;y=±32[1−2(1−q1)9−2(1−q2)9−4Mb(2rc−1)9(rc2+T2)32+13((−234+1μ)σ1+(194−1μ)σ2+(μ1−μ−194)σ1′+(154−μ1−μ)σ2′+2A3)],(10)$

$0 Hence, Eq. 10 are the required coordinates of the triangular equilibrium points denoted by L4,5 (x, ± y). It follows directly that the positions of the triangular points depend on the mass parameter μ and the perturbing forces involved in the equations of motion.

Moreover, we have verified these positions through numerical computations in which we plot in common diagram by making zero the equations of motion Eq. 1 in order to obtain the equilibria which appear as the intersection points of the curves on the xy–plane. In Figure 2, we show as a matter of interest, the location of the collinear and triangular equilibrium points as well as the fixed location of the primaries resulting from the solutions of Eq. 1 for the parameter values μ = 0.4583, Mb = 0.05, σ1 = 4 × 10−5, σ2 = 3 × 10−5, σ3 = 2 × 10−5, σ4 = 10−5, A3 = 0.02, T = 0.01, q2 = 0.99, and q1 = 0.8. Black dots represent primaries, red dots stand for the equilibrium points, and blue curves and orange curves represent the curves Ux = 0 and Uy = 0, respectively. Moreover, our numerical computations (not shown here for latter studies) suggest the influence of triaxiality parameters of the primaries by revealing additional cases, regarding the total number of noncollinear equilibria, especially when these parameters are “pushed” at their limits. Detailed discussions can be made in terms of the variation of such parameters in further studies if needed.

FIGURE 2

FIGURE 2. Plots of Ux = 0 (blue curves) and Uy = 0 (orange curves) depicting the locations of the seven equilibrium points (red dots) on the x y–plane for μ = 0.4583, Mb = 0.05, σ1 = 4 × 10−5, σ2 = 3 × 10−5, σ3 = 2 × 10−5, σ4 = 10−5, A3 = 0.02, T = 0.01, q2 = 0.99, and q1 = 0.8.

The Jacobian integral given by Eq. 6 for μ = 0.455 and $σi=σi′=0,i=1,2,Mb=A3=0,q1=q2=1$ (gravitational case) provides the ZVCs and the corresponding surface shown in Figure 3A,B, respectively, in which the equilibria appear either as intercrossing points between two or more branches of the equipotential curves for a given C or as reducing points after a process of shrinking of closed zero-velocity curves when C varies. The infinitesimal particle is permitted to move inside the “chimney”-like infinitely high peaks at m1 and m2 where C tends to infinity and below it (Figure 3B).

FIGURE 3

FIGURE 3. (Color figure online) (A) Zero-velocity curves for μ = 0.455, $σi=σi′=0$, i = 1, 2, Mb = A3 = 0, q1 = q2 = 1, and Jacobi constant values: $C=CL1$ (red line), $C=CL2$ (green line), C = CL3 (black line), C = 3 (yellow line), and C = 2.8 (blue line). The two primary bodies are denoted by mi i = 1, 2, while the five equilibrium points Li, i = 1,...,5, are indicated by small circles. (B) Surface and the corresponding zero-velocity curves of the problem in the gravitational case.

## Stability of Collinear and Triangular Points

We investigate the linear stability of the infinitesimal motions when small displacements ξ and η are given to the coordinates of an equilibrium point (x0, y0) such that:

$ξ=x−x0, η=y−y0.(11)$

Substituting these values in system Eq. 1, we obtain the linearized system:

$ξ¨−2nη˙=Uxx(0)ξ+Uxy(0)η,η¨+2nξ˙=Uxy(0)ξ+Uyy(0)η.(12)$

Here, only linear terms in ξ and η have been taken. The second partial derivatives of U are denoted by subscripts, while the superscript “0” indicates that the derivatives are to be evaluated at the equilibrium point (x0, y0). Now, the characteristic equation of the linearized system Eq. 12 is:

$λ4+(4n2−Uxx(0)−Uyy(0))λ2+Uxx(0)Uyy(0)−[Uxy(0)]2=0.(13)$

For the collinear equilibrium points $Uxy(0)=Uyx(0)=0$, and stability occurs when all four roots of the characteristic equation for λ are pure imaginary or complex roots all of which have negative real parts. This happens if the following three necessary as well as sufficient conditions are satisfied simultaneously:

$[4n2−Uxx(0)−Uyy(0)]2−4Uxx(0)Uyy(0)>0, 4n2−Uxx(0)−Uyy(0)>0, Uxx(0)Uyy(0)>0.(14)$

In all other cases, the equilibrium points are unstable. In the present problem, the collinear points L1,2,3,n2,n3,n4 with respect to mass parameter μ and relative to involved perturbing forces are unstable since Eq. 13 has two real eigenvalues λ1,2 = ±a and two imaginary eigenvalues λ3,4 = ±ib, while the point Ln1 is stable since the conditions (Eq. 14) are satisfied simultaneously, i.e., Eq. 13 has four purely imaginary eigenvalues of the form λ1,2 = ±ia, λ3,4 = ± ib. Studying the stability of the collinear equilibria for the specific values μ = 0.455, σ1 = 0.01, σ2 = 0.03, σ3 = 0.05, σ4 = 0.07, A3 = 0.001, T = 0.01, Mb = 0.05, q2 = 0.99, and q1 = 0.8, we found that all the collinear equilibria are unstable except point Ln1, which is stable.

For the triangular equilibrium points, the characteristic equation is also given by Eq. 13 but now $Uxy(0)=Uyx(0)≠0,$ so the corresponding inequalities for purely imaginary roots are given by:

$(4n2−Uxx(0)−Uyy(0))2−4[Uxx(0)Uyy(0)−(Uxy(0))2]>0, 4n2−Uxx(0)−Uyy(0)>0, Uxx(0)Uyy(0)−[Uxy(0)]2>0.(15)$

Now, at the triangular points which are given by Eq. 10, we have the following analytical formulas depending only on the parameters of the problem:

$Uxy(0)=Uyx(0)=332{12−μ−[19+μ9](1−q1)+[29−μ9](1−q2)+[4724−89μ24−13μ]σ1+[−38+37μ24+13μ]σ2+[2512−85μ24+μ6(1−μ)+μ26(1−μ)]σ1′+[−32+11μ8−μ6(1−μ)−μ26(1−μ)]σ2′+A3(13−2μ3)+Mb2[11(2rc−1)9(rc2+T2)32+1(rc2+T2)52]−μMb[11(2rc−1)9(rc2+T2)32+1(rc2+T2)52]},Uxx(0)=34−[12−3μ2](1−q1)+[1−3μ2](1−q2)+[5716+45μ16−32μ]σ1+[−316−93μ16+32μ]σ2 +[398−69μ16−3μ22(1−μ)]σ1′+[−92+117μ16+3μ22(1−μ)]σ2′+5Mb(2rc−1)4(rc2+T2)32+3Mb(14−μ+μ2)(rc2+T2)52Uyy(0)=94+3A3+[12−3μ2](1−q1)−[1−3μ2](1−q2)+[8716−45μ16+32μ]σ1+[−2116+45μ16−32μ]σ2 +[338+135μ16−33μ8(1−μ)+45μ28(1−μ)]σ1′+[−135μ16+33μ8(1−μ)−45μ28(1−μ)]σ2′ +7Mb(2rc−1)4(rc2+T2)32+9Mb4(rc2+T2)52(16)$

Substituting these values in Eq. 13, the characteristic equation after normalizing becomes:

$λ4+b1λ2+b2=0,(17)$

where we have abbreviated,

$b1=4n2−Uxx(0)−Uyy(0), b2=Uxx(0)Uyy(0)−Uyx(0)Uxy(0),(18)$

with

$b1=1−3A3+3σ1+(−92+3μ)σ2+3σ1′−3(μ+12)σ2′+Mb(2rc+3)(rc2+T2)32−3Mbrc2(rc2+T2)52,$
$b2=μ(1−μ)[274+32(1−q1)+32(1−q2)+9A3]−[458−891μ16+801μ216]σ1 +[458−423μ16+333μ216]σ2+[711μ16−801μ216]σ1′−[243μ16−333μ216]σ2′ +33μ(1−μ)(2rc−1)Mb2(rc2+T2)32+27μ(1−μ)Mb4(rc2+T2)52.$

Its roots are:

$λ2=−b1±Δ2,(19)$

where the discriminant Δ which defines the nature of roots is expressed as:

$Δ=b12−4b2=μ2(27+36A3+6(1−q1)+6(1−q2)+8014(σ1+σ1′)−3334(σ2+σ2′)+66Mb(2rc−1)(rc2+T2)32+27Mb(rc2+T2)52)−μ(27+36A3+6(1−q1)+6(1−q2)+8914σ1−4474σ2+7114σ1′−2194σ2′+66Mb(2rc−1)(rc2+T2)32+27Mb(rc2+T2)52)+1−6A3+572σ1−632σ2+6σ1′−3σ2′+2Mb(2rc+3)(rc2+T2)32−6Mbrc2(rc2+T2)52$

Now

$(Δ)μ=0=1−6A3+572σ1−632σ2+6σ1′−3σ2'+2Mb(2rc+3)(rc2+T2)32−6Mbrc2(rc2+T2)52>0,$

and

$(Δ)μ=12=−234−15A3−32(1−q1)−32(1−q2)−56516σ1+5716σ2−52516σ1'+5716σ2'−33Mb(2rc−1)2(rc2+T2)32+2Mb(2rc+3)(rc2+T2)32−27Mb4(rc2+T2)52−6Mbrc2(rc2+T2)52<0.$

Evidently, the values of Δ when $μ=0$ and $μ=12$ have opposite signs. Hence, there is a value of $μ$ in the interval $0≤μ≤12$ at which Δ = 0 is called the critical mass parameter denoted by μcrit. Hence, we solve Δ = 0 for $μ$ and simplify to finally get:

$μcrit=μo+μb+μr+μt+μd,(20)$

with

$μo=12(1−2327);μb=−22A3969;μr=−22769[(1−q1)+(1−q2)],μt=(512+591869)σ1−(1936+851869)σ2−(512−591869)σ1′+(1936−851869)σ2′,μd=(76−8rc)Mb2769(rc2+T2)32−(1+6rc2)Mb369(rc2+T2)52.(21)$

Equation 20 represents the perturbed critical mass ratio with five parts (i.e., in terms of all the perturbing parameters). The first part μ0 represents the classical case without perturbation, the second part μb is related to the oblate infinitesimal body, the third part μr is related to radiation pressure of the primaries, and the fourth part μt is connected to the triaxiality of the primaries, while the fifth part μd is related to the circular cluster of material points. Therefore, it is clearly observed that stability regions depend on the perturbing forces of our model. We remark that for Mb = 0, the results agree fully with previously published results (see Gyegwe et al., 2022 and references therein), while in case where $σi =σi′ =0, i=1, 2,$ our results coincide with those by Amuda et al. (2021).

There are three cases (1–3) to determine the linear stability of the triangular equilibrium points. 1) When 0 < μ <μcrit, Δ > 0, the roots Eq. 19 are distinct pure imaginary numbers, and the solution is bounded. Hence, the triangular points are stable in this region. 2) When $μcrit<μ≤12,$ Δ < 0, then the real parts of two of the roots (Eq. 19) are positive and the triangular points are unstable. 3) When μ = μcrit, Δ = 0, then the roots in Eq. 19 are double which leads to instability of the triangular points. Therefore, we established that the triangular points are stable only for 0 < μ < μcrit and unstable for μcritμ ≤ = 1/2. Hence, the stability region is:

$0<μ<μcrit=μo−22A3969−22769[(1−q1)+(1−q2)]+(512+591869)σ1−(1936+851869)σ2−(512−591869)σ1′+(1936−851869)σ2′+(76−8rc)Mb2769(rc2+T2)32−(1+6rc2)Mb369(rc2+T2)52,$

where $μo=12(1−2327)≈0.03852089$ is Routh’s critical value.

This critical mass parameter value gives the measure of the size of the stability region. Computations of Eq. 20 clearly reveal that the critical mass ratio decreases with the increase in the values of the radiation pressure, triaxiality, and oblate infinitesimal body, as shown in Gyegwe et al. (2022); however, it increases with the increase in the value of mass of the disc. It is noticed that all three of the former and the latter one possess destabilizing and stabilizing behavior, respectively. Moreover, under the combined effect of the parameters, the size of the region of stability decreases with a simultaneous increase in the values of the parameters. Remarkably, the present results of the triangular points are in conformity with Yeh and Jiang (2006) for qi = 1,  i = 1, 2, A3 = 0, and $σi=σi′=0,i=1,2,$ in the present work and for Mb = 0, also agree with those of Gyegwe et al. (2022) and references therein; the difference in configuration of the primary bodies is responsible for the difference in sign.

## Numerical Application

In this section, we apply the data of the binary HD155876 system given in Table 1 to the results obtained in our previous sections (Singh et al., 2018). The influence of the involved perturbations on the zero-velocity surface and the positions of the collinear and triangular points will be studied in this section. The numbering of the equilibrium points and the relative positions of the primaries are illustrated according to the diagram shown in Figure 2 for the same mass ratio (μ = 0.4583) of the primaries. Eq. 10 has been used to compute the coordinates of the triangular points.

TABLE 1

TABLE 1. Numerical data for the binary HD155876 system (Singh et al., 2018).

### Effects of the Perturbations on the Zero -Velocity Surfaces

We plot in Figure 4 the zero-velocity surface Eq. 6 for the HD155876 binary system for the purely gravitational case and five different cases of the involved parameters: Panel (A) is for the gravitational case, panel (B) shows the oblateness effect of the third body, panel (C) shows the effect of enclosed material points in which the immediate formation of a new “chimney” around the bigger primary is obvious, panel (D) shows the effect of triaxiality of the bigger primary, panel (E) shows the effect of triaxiality of the smaller primary, and panel (F) shows the effect of the radiation coefficients of the binary system, while panel (G) shows the combine effects of triaxiality and radiation pressure of both primaries, oblateness of the third body, and enclosed material points of the primaries. From these frames, we can observe two main results. The first is a change of the Jacobian constants, which correspond to all equilibrium points with or without involved parameters due to which results of zero-velocity surfaces show the chaotic behavior. Note the different values of the Jacobian constant C of the panels (B)–(G) w.r.t panel (A). The second is a shrinking or enlargement of the “chimneys” around each primary which limit the area of the permitted third body’s motion in the neighborhood of the two primaries for each of the panels (B)–(G) w.r.t panel (A). We note two collinear equilibria which emerge on side of primary body $m1,$ in addition to the three classical points [see panels (C) and (G)]. It is also remarkable that for a given value of the Jacobian constant C and as the involved parameters increase, the allowed regions of the third body motion around the two bodies expand [panel (G)]. This means that the third body can move freely around the primaries for bigger and bigger values of the Jacobian constant C as the parameters increase.

FIGURE 4

FIGURE 4. (Color figure online) Zero-velocity surfaces of the HD155876 (μ = 0.4583) binary system for (A) A3 = 0, σ1 = σ2 = 0, $σ1′=σ2′=0$, Mb = 0, and q1 = q2 = 1 (classical case); (B) A3 = 0.09, σ1 = σ2 = 0, $σ1′=σ2′=0$, Mb = 0, and q1 = q2 = 1 case; (C) Mb = 0.05, T = 0.01, A3 = 0, σ1 = σ2 = 0, $σ1′=σ2′=0$, and q1 = q2 = 1 case; (D) σ1 = 0.02, σ2 = 0.01, $σ1′=σ2′=0$, A3 = 0, Mb = 0, and q1 = q2 = 1 case; (E) $σ1′=0.02$, $σ2′=0.01$, σ1 = σ2 = 0, A3 = 0, Mb = 0, and q1 = q2 = 1 case; (F) q1 = 0.5, q2 = 0.8, A3 = 0, σ1 = σ2 = 0, $σ1′=σ2′=0$, and Mb = 0 case; (G) q1 = 0.45, q2 = 0.75, A3 = 0.1, σ1 = 0.07, σ2 = 0.05, $σ1′=0.03$, $σ2′=0.02$, Mb = 0.055, and T = 0.01 case. The presence of the new collinear points is obvious in panels (C) and (G). Motion is permitted inside the chimneys and below it. Note the different scales of Jacobian constant C in the first and the next frames.

It is obvious from these figures that the perturbing forces under consideration have a significant effect in which the motion of the test particle is allowed or forbidden in the vicinity of the HD155876 binary system, comparing the first frame with the six last ones. Remarkably, the cluster of material points of the primaries provides a new observation of its topological structure resulting in the neighborhood of $m1$ primary a pair of new collinear equilibrium points. Here, we note that since the two primaries have unequal masses and unequal assumed model parameters, the two “chimneys” are asymmetric with respect to both (C, y) and (C, x) planes.

### Influence of [q1, q2] on the Equilibria

For the arbitrary chosen values of the model parameters $A3,σi,σi′,i=1,2,$ and Mb, the locations of the collinear and triangular equilibrium points are, respectively, presented in Tables 2 and 3 as a function of q2 w.r.t. different values of q1. Evidently, for fixed values of q1 and increasing values of q2, the coordinates x and y of the equilibria decrease or increase. Here, the increase or decrease is a measure of the deviation from the classical estimations. Specifically, with the increase of q2, the collinear equilibrium positions L2 and L3 constantly shift their distance and move closer to the position where primary bodies $m1$ and $m2$ are located, correspondingly. L1 and Ln1 go closer to $m2$ and $m1$, correspondingly, while at the same time collinear point Ln2 moves closer to the origin. At the same time, the two triangular points L4 and L5 decrease with the line connecting the two primaries.

TABLE 2

TABLE 2. Positions (x0, 0) of the collinear points of the binary HD155876 system for the fixed values Mb = 0.05, T = 0.01, σ1 = 4 × 10−5, σ2 = 3 × 10−5, $σ1′=2×10−5$, $σ2′=10−5$, A3 = 0.02, and for two fixed values of q1 when q2 radiates.

TABLE 3

TABLE 3. Positions (x0, ± y0) of the triangular points of the binary HD155876 system for three fixed values of q1 when q2 radiates. The disc, oblateness, and triaxiality parameters are kept as constants as in the collinear case (Table 2).

Such a variational trend can also be seen from Figures 5 and 6 in which the equilibria positions have been marked with curves corresponding to the cases with q1 = 1, q1 = 0.5, and q1 = 0.1 for the increasing values of radiation factor q2 of the smaller primary. Here, we note from the tables and figures that the positions of Ln1, L4, and L5 are decreasing and increasing functions of q2. A similar phenomenon is observed under the combined effects of q1 and q2 in which the positions of Ln1, L4, and L5 are monotonic functions of q1 and q2.

FIGURE 5

FIGURE 5. (A–E) Positions of the collinear equilibrium points Li, i = 1, 2, 3, Lnj, j = 1, 2, of the HD11565 system, as a function of q2, which varies in the interval (0, 1] for fixed values of q1 equal to 1 (1), 0.5 (2), and 0.1 (3). The disc, oblateness, and triaxiality parameters are kept as constants as given in Table 2. The numbering of equilibrium points is according to the diagram shown in Figure 2.

FIGURE 6

FIGURE 6. Positions of the two triangular equilibrium points L4 and L5 on the xy–plane as a function of q2, which varies in the interval (0, 1] for fixed values of q1 equal to 1 (1), 0.5 (2), and 0.1 (3). The disc, oblateness, and triaxiality parameters are kept as constants as given in Table 3.

Similarly, the locations of the collinear and triangular points are shown in Tables 4 and 5, respectively, as a function of q1 with respect to different values of q2. We observe that as the radiation coefficient q2 increases for varying q1, collinear equilibrium positions L2 and L3 decrease and move closer to the primaries $m2$ and $m1$, correspondingly, and points Ln1 and Ln2 both tend to $m2$ and $m1$, respectively, while point L1 moves toward the origin. At the same time, the two triangular equilibria L4 and L5 go toward the line connecting the primaries as q1 tends to zero for fixed values of q2.

TABLE 4

TABLE 4. Positions (x0, 0) of the collinear points of the binary HD155876 system for fixed values Mb=0.05, T = 0.01, σ1 = 4 × 10−5, σ2 = 3 × 10−5, $σ1′=2×10−5$, $σ2′=10−5$, A3 = 0.02, and for two fixed values of q2 when q1 radiates.

TABLE 5

TABLE 5. Positions (x0, ± y0) of the triangular points of the binary HD155876 system for three fixed values of q2 when q1 radiates. The disc, oblateness, and triaxiality parameters are kept as constants as in the collinear case (Table 4).

Such a variational trend can also be seen from Figures 7 and 8 in which the equilibria positions have been marked with curves corresponding to the cases with q2 = 1, q2 = 0.5, and q2 = 0.1 for values of the radiation factor of the bigger primary q1 ∈ (0, 1]. We also observe from the tables and figures that the positions of Ln1, L4, and L5 are decreasing functions of q1. A similar phenomenon is observe under the combined effects of q1 and q2 in which the positions of Ln1, L4, and L5 are monotonic functions of q1 and q2. It is also remarkable that under the combined effects of the radiation factors q1 and q2, the position (x, 0) of the equilibrium points L1 and Ln2 is an increasing function and that of L2 and L3 is a decreasing function, while the x coordinate of Ln1, L4, and L5 is a monotonic function of q1 and q2 (Tables 4, 5 and Figures 7, 8).

FIGURE 7

FIGURE 7. (A–E) Positions of the collinear equilibria Li, i = 1, 2, 3, Lnj, j = 1, 2, of the HD11565 system, as a function of q1, which varies in the interval (0, 1] for fixed values of q2 equal to 1 (1), 0.5 (2), and 0.1 (3), respectively. The disc, oblateness, and triaxiality parameters are kept as constants as given in Table 4. The numbering of equilibrium points is according to the diagram shown in Figure 2.

FIGURE 8

FIGURE 8. Positions of the two triangular equilibrium points L4 and L5 on the xy–plane as a function of q1, which varies in the interval (0, 1] for fixed values of q2 equal to 1 (1), 0.5 (2), and 0.1 (3). The disc, oblateness, and triaxiality parameters are kept as constants as given in Table 5.

### Influence of [A3, Mb] on the Equilibria

In order to analyze the effects of oblate infinitesimal mass A3 and mass of the disc Mb, we have obtained the coordinates of the five collinear and the two triangular equilibria for varying oblateness and varying mass of the disc at fixed values of q1, q2, σ1, σ2, $σ1′$, and $σ2′$ (Tables 6, 7). From the results presented in Table 6, it is clear that with an increase in the value of A3 for fixed values of Mb, points L2 and L3 tend to infinity (L2 → ∞, L3 → − ∞), and L1 and Ln2 tend closest to the origin, while the inner point Ln1 tends to the position of $m1,$ correspondingly. On the other hand, with an increase in the value of A3 for fixed values of Mb, the two triangular points L4 and L5 both tend to infinity (L4 → ∞, L5 → − ∞) (Table 7).

TABLE 6

TABLE 6. Positions (x0, 0) of the collinear points of the binary HD155876 system for fixed values q1 = 0.979950, q2 = 0.983912, T = 0.01, σ1 = 4 × 10−5, σ2 = 3 × 10−5, $σ1′=2×10−5$, $σ2′=10−5$, and for two fixed values of Mb when the value of A3 varies.

TABLE 7

TABLE 7. Positions (x0, ±y0) of the two triangular points L4 and L5 of the binary HD155876 system for three fixed values of Mb when A3 varies. The radiation pressure and triaxiality parameters are kept as constants as in the collinear case (Table 6).

Such variational trends for graphical solutions of the collinear and triangular equilibria are shown in Figures 9 and 10, respectively. From Tables 6 and 7, it is clearly seen that for fixed values of A3 and varying values of Mb, points L1 and L2 move closest to the position of $m2,$ L3 and L1 tend to the position of $m1,$ and L1 moves closer to the origin, while L4 and L5 tend to the line joining the primaries. Due to the limiting space of the article, we do not present here this case as the presentation is obvious with regard to the previous case. Moreover, in the absence and presence of the remaining parameters, changes in the value of Mb or A3 do not influence the x-coordinates of the triangular points analytically (Table 7 and Figure 10). However, under the joint action of both parameters, the situation is analogous for variation of Mb and fixed values of A3. This indicates that the presence of the disc plays a certain degree of “correction” to the two primaries’ dynamic behavior, so that they can maintain regular motion as much as possible.

FIGURE 9

FIGURE 9. (A–E) Positions of the collinear equilibrium points Li, i = 1, 2, 3, Lnj, j = 1, 2, of the HD155876 system, as a function of A3, which varies in the interval (0, 0.1] for fixed values of Mb equal to 0.06 (1), 0.1 (2), and 0.15 (3), respectively. The radiation pressure and triaxiality parameters are kept as constants as in the previous case (Table 6). The numbering of equilibrium points is according to the diagram shown in Figure 2.

FIGURE 10

FIGURE 10. Positions of the two triangular equilibrium points L4 and L5 on the xy–plane as a function of A3, which varies in the interval (0, 0.1) for fixed values of Mb equal to 0.06 (1), 0.1 (2), and 0.15 (3). The radiation pressure and triaxiality parameters are kept as constants as in Table 7.

### Influence of Triaxiality of the Bigger Primary (σ1, σ2) on the Equilibria

Tables 8 and 9 evaluate triaxiality effects of the bigger primary in the binary system on the collinear and triangular equilibrium points, respectively, while in Figures 11 and 12, we summarize these effects on their positions, respectively. Evidently, an increase in the value of σ1 for fixed values of σ2, points L1 and L2 tend toward the smaller primary m2, Ln1 and Ln2 approach each other, while L3 tends to minus infinity, whereas the two triangular points L4 and L5 both tend to the Ox–axis.

TABLE 8

TABLE 8. Positions (x0, 0) of the collinear points of the binary HD155876 system for fixed values q1 = 0.979950, q2 = 0.983912, A3 = 0.02, Mb = 0.06, T = 0.01, $σ1′=2×10−5$, $σ2′=10−5$, and for two fixed values of $σ2$ when $σ1$ varies.

TABLE 9

TABLE 9. Positions (x0, ± y0) of the two triangular points L4 and L5 of the binary HD155876 system for three fixed values of σ2 when σ1 varies. The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the smaller primary are kept as constants as in the collinear case (Table 8).

FIGURE 11

FIGURE 11. (A–E) Positions of the collinear equilibrium points Li, i = 1, 2, 3, Lnj, j = 1, 2, of the HD155876 system, as a function of σ1, which varies in the interval (0.02, 0.12) for fixed values of σ2 equal to 0.01 (1), 0.05 (2), and 0.09 (3), respectively. The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the smaller primary are kept as constants as in the previous case (Table 8). The numbering of equilibrium points is according to the diagram shown in Figure 2.

FIGURE 12

FIGURE 12. Positions of the two triangular equilibrium points L4 and L5 on the xy–plane as a function of σ1, which varies in the interval (0.02, 0.12) for fixed values of σ2 equal to 0.01 (1), 0.05 (2) and 0.09 (3). The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the smaller primary are kept as constants as in the previous case (Table 9).

### Influence of Triaxiality of the Smaller Primary ($σ1′,σ2′$) on the Equilibria

Tables 10 and 11evaluate triaxiality effects of the smaller primary in the binary system on the collinear and triangular equilibrium points, respectively, while in Figures 13 and 14, we summarize these effects on their positions, respectively. Evidently, an increase in the value of $σ1′$ for fixed values of $σ2′$, points L1 and Ln1 both tend toward the origin, L3 and Ln2 both approach the bigger primary m1, while L2 tends to plus infinity, whereas the two triangular points L4 and L5 both move toward the Ox–axis.

TABLE 10

TABLE 10. Positions (x0, 0) of the collinear points of the binary HD155876 system for fixed values q1 = 0.979950, q2 = 0.983912, A3 = 0.02, Mb = 0.06, T = 0.01, σ1 = 4 × 10−5, σ2 = 3 × 10−5, and for two fixed values of $σ2′$ when $σ1′$ varies.

TABLE 11

TABLE 11. Positions (x0, ±y0) of the two triangular points L4 and L5 of the binary HD155876 system for three fixed values of $σ2′$ when $σ1′$ varies. The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the bigger primary are kept as constants as in the collinear case (Table 10).

FIGURE 13

FIGURE 13. (A–E) Positions of the collinear equilibrium points Li, i = 1, 2, 3, Lnj, j = 1, 2, of the HD155876 system, as a function of $σ1′$, which varies in the interval (0.01, 0.06) for fixed values of $σ2′$ equal to 0.02 (1), 0.04 (2), and 0.06 (3), respectively. The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the bigger primary are kept as constants as in the previous case (Table 10). The numbering of equilibrium points is according to the diagram shown in Figure 2.

FIGURE 14

FIGURE 14. Positions of the two triangular equilibrium points L4 and L5 on the xy–plane as a function of $σ1′$, which varies in the interval (0.01, 0.06) for fixed values of $σ2′$ equal to 0.02 (1), 0.04 (2), and 0.06 (3). The radiation pressure, oblateness, mass of the disc, and triaxiality parameters of the bigger primary are kept as constants as in the previous case (Table 11).

Now, we shall analyze the linear stability of the collinear and triangular equilibrium points by looking at the eigenvalues of the characteristic Eq. 13. In this article, for a particular example, we compute the stability of the collinear and triangular points under the joint effect of oblate infinitesimal body and mass of the disc of the binary system HD155876 when the remaining parameters q1, q2, σ1, σ2, $σ1′$, and $σ2′$ are fixed. The results are summarized in Table 12 in which the presence of positive real roots or positive real part in complex roots affirms the instability of L1,2,3,n2 and L4(5) of the problem when applied to the HD155876 system. However, for point Ln1, we get all purely imaginary roots implying that due to all imaginary roots at this point are stable.

TABLE 12

TABLE 12. Locations and eigenvalues λ1,2, λ3,4 of Eq. 13 in the HD155876 binary system with varying oblateness and mass of the disc for q1 = 0.979950, q2 = 0.983912, σ1 = 4 × 10−5, σ2 = 3 × 10−5, $σ1′=2×10−5$, and $σ2′=10−5$.

## Discussion and Conclusion

Following the equations of motion of an oblate body with negligible mass moving under the gravitational attraction of radiating primaries coupled with triaxiality of the two stars together with a circular cluster of materials, we determined the positions of the collinear and triangular equilibrium points. The resulting equations of motion differ from those described by Gyegwe et al. (2022) due to the potential from the belt. It was observed that the positions of these points are affected by the mass parameter and perturbing forces involved in the equations of motion. Moreover, it was found that in addition to the collinear equilibrium points of the classical R3BP, there emerge four additional collinear points considering very small ranges of asphericity of the primaries and potential from the belt. It was also noted that if the potential from the belt is neglected, the obtained semi-analytical formulas which describe the coordinates of the triangular equilibrium points coincide with those presented by Gyegwe et al. (2022).

Despite the presence of the perturbing forces, all the collinear points remain unstable except Ln1. The triangular points are stable for 0 <μ < μcrit and unstable for μcritμ ≤ 1/2, where μcrit is the critical mass ratio which depends on the combined effect of the perturbing forces. The resultant Eq. 20 for the critical mass clearly describes the effects of the perturbing forces on it. Again, if the gravitational potential from the belt is neglected, μcrit confirms the result of Gyegwe et al. (2022). Equation 20 indicates that the critical mass ratio decreases with the increase in the values of the radiation pressure, triaxiality, and oblate infinitesimal body, while it increases with the increase in the value of mass of the disc (namely the range of stability). Hence, all three of the former and the latter one possess destabilizing and stabilizing behavior, respectively.

Finally, a numerical exploration, using the binary HD155876 system, was performed to locate the collinear and triangular points of the system as well as the allowed regions of motion of the test body. These points were shown numerically and graphically, thus highlighting the effects of the involved parameters. It was found that the position of the collinear and triangular points is affected in presence of perturbations because it is deviated from the classical R3BP on the x–axis and out of the x–axis, respectively. Likewise, allowed regions of motion of the infinitesimal body are significantly affected with respect to the gravitational case. The collinear points L1,2,3,n2,n3,n4 remain unstable except Ln1, while triangular points L4,5 are stable under certain conditions.

In conclusion, we would like to report that the mass parameter μ has a minor impact on the evolution of the equilibria, with respect to the influence of the triaxiality parameters contrary to the mass of the disc. The presented results may help to analyze more generalized problems of few bodies under the influence of different kinds of perturbations such as Stokes drag, Poynting–Robertson drag, and Yarkovsky or Albedo effects, which are in relevance to this study, in a general theoretical sense of studying nongravitational nature phenomena regarding nonclassical effects (see, for e.g., Ershkov, 2012; Abouelmagd and Sharaf, 2013; Idrisi and Ullah, 2020; Ershkov et al., 2021).

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

## Author Contributions

AV contributed to conception and investigation of the study as well as wrote the first draft of the manuscript. AP and EP contributed to formal analysis, validation, and supervision of the study. All authors contributed to manuscript revision and read and approved the submitted version.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Abouelmagd, E. I., Awad, M. E., Elzayat, E. M. A., and Abbas, I. A. (2014). Reduction the Secular Solution to Periodic Solution in the Generalized Restricted Three-Body Problem. Astrophys. Space Sci. 350, 495–505. doi:10.1007/s10509-013-1756-z

Abouelmagd, E. I., and Sharaf, M. A. (2013). The Motion Around the Libration Points in the Restricted Three-Body Problem with the Effect of Radiation and Oblateness. Astrophys. Space Sci. 344, 321–332. doi:10.1007/s10509-012-1335-8

Amuda, T. O., Singh, J., and Oni, L. (2021). Motion Around Equilibrium Points of an Oblate Body in the PR3BP with Disc. Indian J. Phys. 95, 1305–1315. doi:10.1007/s12648-020-01799-z

Capdevila, L. R., and Howell, K. C. (2018). A Transfer Network Linking Earth, Moon, and the Triangular Libration point Regions in the Earth-Moon System. Adv. Space Res. 62, 1826–1852. doi:10.1016/j.asr.2018.06.045

Chermnykh, S. V. (1987). Stability of Libration Points in a Gravitational Field. Leningradskii Universitet Vestnik Matematika Mekhanika Astronomiia 2, 73–77.

Contopoulos, G. (2002). Order and Chaos in Dynamical Astronomy. Berlin: Springer-Verlag.

Douskos, C., Kalantonis, V., and Markellos, P. (2007). Effects of Resonances on the Stability of Retrograde Satellites. Astrophys. Space Sci. 310, 245–249. doi:10.1007/s10509-007-9508-6

Elshaboury, S. M., Abouelmagd, E. I., Kalantonis, V. S., and Perdios, E. A. (2016). The Planar Restricted Three-Body Problem when Both Primaries Are Triaxial Rigid Bodies: Equilibrium Points and Periodic Orbits. Astrophys. Space Sci. 361, 315. doi:10.1007/s10509-016-2894-x

Ershkov, S., Leshchenko, D., and Abouelmagd, E. I. (2021). About Influence of Differential Rotation in Convection Zone of Gaseous or Fluid Giant Planet (Uranus) onto the Parameters of Orbits of Satellites. Eur. Phys. J. Plus 136, 387. doi:10.1140/epjp/s13360-021-01355-6

Ershkov, S. V. (2012). The Yarkovsky Effect in Generalized Photogravitational 3-body Problem. Planet. Space Sci. 73, 221–223. doi:10.1016/j.pss.2012.09.002

Gao, F., and Wang, R. (2020). Bifurcation Analysis and Periodic Solutions of the HD 191408 System with Triaxial and Radiative Perturbations. Universe 6, 35. doi:10.3390/universe6020035

Goździewski, K., and Maciejewski, A. J. (1998). Nonlinear Stability of the Lagrangian Libration Points in the Chermnykh Problem. Celest. Mech. Dyn. Astron. 70, 41–58.

Gyegwe, J. M., Vincent, A. E., and Perdiou, A. E. (2022). “On the Stability of the Triangular Equilibrium Points in the Photogravitational R3BP with an Oblate Infinitesimal and Triaxial Primaries for the Binary Lalande 21258 System,” in Approximation and Computation in Science and Engineering. Editors N. J. Daras, and Th. M. Rassias (Springer Optimization and Its Applications). to be published.

Hou, X. Y., and Liu, L. (2008). Vertical Bifurcation Families from the Long and Short Period Families Around the Equilateral Equilibrium Points. Celest Mech. Dyn. Astr 101, 309–320. doi:10.1007/s10569-008-9147-4

Idrisi, M. J., and Ullah, M. S. (2020). A Study of Albedo Effects on Libration Points in the Elliptic Restricted Three-Body Problem. J. Astronaut. Sci. 67, 863–879. doi:10.1007/s40295-019-00202-2

Jiang, I.-G., and Yeh, L.-C. (2006). On the Chermnykh-like Problems: I. The Mass Parameter μ = 0.5. Astrophys. Space Sci. 305, 341–348. doi:10.1007/s10509-006-9065-4

Kalantonis, V. S. (2020). Numerical Investigation for Periodic Orbits in the Hill Three-Body Problem. Universe 6, 72. doi:10.3390/universe6060072

Kalantonis, V. S., Vincent, A. E., Gyegwe, J. M., and Perdios, E. A. (2021). “Periodic Solutions Around the Out-Of-Plane Equilibrium Points in the Restricted Three-Body Problem with Radiation and Angular Velocity Variation,” in Nonlinear Analysis and Global Optimization. Editors Th. M. Rassias, and P. M. Pardalos (Cham: Springer), 251–275. Springer Optimization and Its Applications. doi:10.1007/978-3-030-61732-5_11

Kishor, R., and Kushvah, B. S. (2013). Linear Stability and Resonances in the Generalized Photogravitational Chermnykh-like Problem with a Disc. Mon. Not. R. Astron. Soc. 436, 1741–1749. doi:10.1093/mnras/stt1692

Leke, O., and Singh, J. (2021). Exploring Effect of Perturbing Forces on Periodic Orbits in the Restricted Problem of Three Oblate Spheroids with Cluster of Material Points. Int. Astron. Astrophysics Res. J. 2 (4), 48–73.

Markellos, V. V., Papadakis, K. E., and Perdios, E. A. (1996). Non-linear Stability Zones Around Triangular Equilibria in the Plane Circular Restricted Three-Body Problem with Oblateness. Astrophys. Space Sci. 245, 157–164. doi:10.1007/bf00637811

Marsden, J., and Ross, S. (2005). New Methods in Celestial Mechanics and mission Design. Bull. Amer. Math. Soc. 43, 43–73. doi:10.1090/s0273-0979-05-01085-2

Musielak, Z., and Quarles, B. (2017). “Three Body Dynamics and its Applications to Exoplanets,” in Springer Briefs in Astronomy (Cham, Switzerland: Springer). doi:10.1007/978-3-319-58226-9

Papadakis, K. E. (2005). Numerical Exploration of Chermnykh's Problem. Astrophys. Space Sci. 299, 67–81. doi:10.1007/s10509-005-3070-x

Perdios, E. A., Kalantonis, V. S., Perdiou, A. E., and Nikaki, A. A. (2015). Equilibrium Points and Related Periodic Motions in the Restricted Three-Body Problem with Angular Velocity and Radiation Effects. Adv. Astron 2015, 473–483. doi:10.1155/2015/473483

Perdios, E. A., and Ragos, O. (2004). Asymptotic and Periodic Motion Around Collinear Equilibria in Chermnykh's Problem. A&A 414, 361–371. doi:10.1051/0004-6361:20031619

Perdiou, A. E., Nikaki, A. A., and Perdios, E. A. (2013). Periodic Motions in the Spatial Chermnykh Restricted Three-Body Problem. Astrophys. Space Sci. 345, 57–66. doi:10.1007/s10509-013-1368-7

Prosmiti, R., Farantos, S. C., Guantes, R., Borondo, F., and Benito, R. M. (1996). A Periodic Orbit Analysis of the Vibrationally Highly Excited LiNC/LiCN: A Comparison with Quantum Mechanics. J. Chem. Phys. 104, 2921–2931. doi:10.1063/1.471113

Saeed, T., and Zotos, E. E. (2021). On the Equilibria of the Restricted Three-Body Problem with a Triaxial Rigid Body - I. Oblate Primary. Results Phys. 23, 103990. doi:10.1016/j.rinp.2021.103990

Sano, M. M. (2007). Dynamics Starting from Zero Velocities in the Classical Coulomb Three-Body Problem. Phys. Rev. E Stat. Nonlin Soft Matter Phys. 75, 026203. doi:10.1103/PhysRevE.75.026203

Singh, J., and Amuda, T. O. (2019). Stability Analysis of Triangular Equilibrium Points in Restricted Three-Body Problem under Effects of Circumbinary Disc, Radiation and Drag Forces. J. Astrophys. Astron. 40, 5. doi:10.1007/s12036-019-9573-6

Singh, J., and Leke, O. (2014). Analytic and Numerical Treatment of Motion of Dust Grain Particle Around Triangular Equilibrium Points with post-AGB Binary star and Disc. Adv. Space Res. 54, 1659–1677. doi:10.1016/j.asr.2014.06.031

Singh, J., Perdiou, A. E., Gyegwe, J. M., and Perdios, E. A. (2018). Periodic Solutions Around the Collinear Equilibrium Points in the Perturbed Restricted Three-Body Problem with Triaxial and Radiating Primaries for Binary HD 191408, Kruger 60 and HD 155876 Systems. Appl. Maths. Comput. 325, 358–374. doi:10.1016/j.amc.2017.11.052

Singh, J. (2013). The Equilibrium Points in the Perturbed R3BP with Triaxial and Luminous Primaries. Astrophys. Space Sci. 346, 41–50. doi:10.1007/s10509-013-1420-7

Szebehely, V. (1967). Theory of Orbits. The Restricted Problem of Three Bodies. New York: Academic Press.

Valtonen, M., and Karttunen, H. (2006). The Three-Body Problem. Cambridge: Cambridge University Press.

Voyatzis, G., and Antoniadou, K. I. (2018). On Quasi-Satellite Periodic Motion in Asteroid and Planetary Dynamics. Celest. Mech. Dyn. Astr. 130, 59. doi:10.1007/s10569-018-9856-2

Voyatzis, G., Gkolias, I., and Varvoglis, H. (2012). The Dynamics of the Elliptic Hill Problem: Periodic Orbits and Stability Regions. Celest. Mech. Dyn. Astr. 113, 125–139. doi:10.1007/s10569-011-9394-7

Vrahatis, M. N., Perdiou, A. E., Kalantonis, V. S., Perdios, E. A., Papadakis, K., Prosmiti, R., et al. (2001). Application of the Characteristic Bisection Method for Locating and Computing Periodic Orbits in Molecular Systems. Comput. Phys. Commun. 138, 53–68. doi:10.1016/s0010-4655(01)00190-4

Yeh, L.-C., and Jiang, I.-G. (2006). On the Chermnykh-like Problems: II. The Equilibrium Points. Astrophys. Space Sci. 306, 189–200. doi:10.1007/s10509-006-9170-4

Zotos, E. E., Chen, W., Abouelmagd, E. I., and Han, H. (2020). Basins of Convergence of Equilibrium Points in the Restricted Three-Body Problem with Modified Gravitational Potential. Chaos, Solitons & Fractals 134, 109704. doi:10.1016/j.chaos.2020.109704

Keywords: R3BP, equilibrium points, stability, triaxiality, circumbinary disc

Citation: Vincent AE, Perdiou AE and Perdios EA (2022) Existence and Stability of Equilibrium Points in the R3BP With Triaxial-Radiating Primaries and an Oblate Massless Body Under the Effect of the Circumbinary Disc. Front. Astron. Space Sci. 9:877459. doi: 10.3389/fspas.2022.877459

Received: 16 February 2022; Accepted: 24 February 2022;
Published: 27 April 2022.

Edited by:

Elbaz I. Abouelmagd, National Research Institute of Astronomy and Geophysics, Egypt

Reviewed by:

Sergey Ershkov, Plekhanov Russian University of Economics, Russia
Rajiv Aggarwal, University of Delhi, India
Amit Mittal, University of Delhi, India

Copyright © 2022 Vincent, Perdiou and Perdios. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: A. E. Perdiou, aperdiou@upatras.gr