Abstract

The paper investigates motion of a dust grain around collinear equilibrium points (CEPs) of the circular restricted three-body problem (R3BP), with variable masses and disc. The masses of the binary vary with time as defined by the combined Mestschersky law (CML) and their motion is governed by the Gylden-Mestschersky problem (GMP), with the assumption that the system is capsulated by a disc of dust whose masses are also assumed to vary with time under the CML. The equations of motion are deduced, and the CEPs are investigated. It is seen that there are four or five CEPs depending jointly on the mass variation parameter and mass of the disc. These points are characterized by the mass parameter, mass of disc and the mass variation constant. The time dependent CEPs are different from those of the autonomized locations by a function of time. Additionally, the stability of the CEPs of the autonomized and non-autonomous systems is investigated, and it was seen that these points are unstable. Illustrative numerical examples are given to further explore the dynamics of the particles around the CEPs of the binary MAXI J1659-152. The positions, zero velocity curves (ZVCs), Poincare surface of section (PSS) and stability are analyzed numerically, and the impact of the variable disc is shown. We conclude that the accreting disc in the presence of mass variations of the stars does not automatically result in the existence of two additional CEPs and at some points, the locations of the CEPs remain unchanged and are not influenced by mass gain of the binary and disc. Further, it seen for the ZVCs, that the accreting disc reduces region where motion in allowed and the PSS reveals that, orbits around the CEPs for very small mass of the disc is unstable. This model can be used to study the long-term motion of satellites and planets in binary systems with gradual accretion of disc of dust in the environment.

Keywords:R3BP; Equilibrium Points; Disc; Variable masses; ZVC; PSS; MAXI J1659-152

Introduction

One of the basic model systems of celestial mechanics with many applications is the restricted three-body problem (R3BP). This problem makes up one of the most studied areas in dynamical system theory and is the most searched and interesting problem for astrophysicists. The study of the R3BP in its many classifications has had important implications in several scientific fields, including galactic dynamics, celestial mechanics, molecular physics and chaos theory, among others. Also, solar and stellar systems dynamics and motion of artificial satellites form some of the areas where the R3BP is applied. The classical R3BP has five equilibrium points [1]. Three points lie on the line joining the primaries and two form triangles with the primaries. Now, the classical R3BP considered the masses of the bodies to be constant, however, in the real sense “nothing” is indeed constant. Masses of celestial bodies change with time, especially in binary system where mass changes are pronounced. The model of motion of astronomical objects with variable mass has lot of interesting applications in stellar, galactic, and planetary dynamics. As an example, we can mention the motion of a satellite travelling around a radiating star enclosed by a cloud or disc dust and varying its mass due to particles around it. Another fascinating example of mass change is the real physical scenario of those transiting exoplanets such as HD 209458b whose atmospheres are escaping because of the severe levels of energetic radiation, coming from their very close parent stars, hitting them.

Various studies under different characterizations have investigated R3BP with mass variations of both primaries. Some of these studies are confined on when motion of the primaries is defined by the Gylden-Mestschersky problem (GMP) [2,3] while their mass variations take place within the dictates of the combined Mestschersky Law (CML) [4,5] investigated the R3BP of variable mass in which the primary bodies move within the framework of the GMP and assumed that the isotropic mass variation of the masses of the primaries occurs in accordance with CML. He established the existence of three collinear and two triangular EPs analogous to the classical type. The existence of a pair of out-of-plane or coplanar EPs was found by [6] with same formulation of Gelf’gat (1973). Other researchers that have adopted the same formulations under different characterizations include [7-13]. Recently, [14] and, [15] investigated the Out-ofplane and CEPs when there is a disc to interact with and when the bigger primary is a triaxial body, respectively.

Further, the classical R3BP did not include the presence of circular clusters of materials, which are often referred to as disc. Circumbinary discs (CBD hereafter) are discs of matter that orbits externally to a central binary system that is typically composed of two stars or black holes. They may form in a variety of astrophysical systems, including when stellar binaries capture material in dense star forming regions or in galactic centres where, following a galaxy merger, two super massive black holes (SMBH) accrete gas from the host galaxy. An accretion disc is a swirling disc of gas and dust orbiting a star or black hole. CBD accretion plays an important role in the evolution of many types of astrophysical systems, ranging from young binary stars, main sequence, and post main sequence binaries to super massive binary black holes. Circumbinary accretion has long been suggested to exist around binary massive black holes (MBHs) following galaxy mergers [16]. An additional significant clue for the existence and importance of CBD is given by observed circumbinary planets. These are planetary systems where the planets do not orbit one star, but instead a binary star system. An example is Kepler-16b which orbits its host system consisting of a K-type main-sequence star and an M-type red dwarf, in what is known as a P-type orbit [17]. Numerous discs have been observed around binaries like L14481RS3B (Tobin et al 2016). Also, around young stellar binaries, like GG Tau, DQ Tau and UZ Tau E (Phuong et al 2020) and around binaries like IRAS16293-2423A (Maureira et al 2020).

[18] discussed the orbital evolution of binaries with CD and concluded that for common CD parameters, binaries with non-extreme mass ratios are expected to shrink over time. [19] maintained that for a CD, the results are parameter dependent, with thinner discs causing the binary to shrink, while thicker discs cause the binary to expand. The importance of the inclusion of disc in the study of the R3BP have been investigated by several researchers, amongst whom are where it was shown that the presence of a disc resulted in the existence of two additional CEPs. Inspired by the study of Leke and Singh (2023) where the Out-ofplane EPs of extra-solar planets in the central binaries PSR B1620- 26 and Kepler-16 with cluster of material points and variable masses, was investigated, and the work of Leke and Orum (2024), in which the collinear EPs of the R3BP with variable masses and a triaxial primary was explored. In this paper, our attention is drawn to the study of motion of dust grain particle around CEPs of the R3BP with variable masses and disc. The structure of the paper is as follows. The dynamical equations are deduced in Section 2, while the locations of the CEPs are investigated in Section 3. Section 4 deals with the stability of the CEPs while the numerical simulations of the problem and the discussion are presented in Section 5. The conclusion is given in Section 6.

Equations Of Motion

To deduce the equations of motion of the R3BP, we consider a rotating frame of reference Oxyz , where O is taken as the origin and (x, y, z) is the coordinate of. Next, we let r₁ and r₂ be the distances between the third body and the primaries m₁ and m₂ , respectively. The primaries are placed along x − the axis and positioned at (x₁, 0, 0) and, (x₂, 0, 0) respectively. We assume that there is a disc of dust around the configuration, having variable masses. On this basis, the equations of motion have the form (Leke and Singh 2023):

M , ℜ , a and b are the varying mass, radius, radial scalelength, and vertical scale-height, of the disc, respectively, while is the gravitational constant and is the product of sum of the masses of the primaries and the gravitational constant and c r is the time dependent radial distance of the infinitesimal mass κ (0 <κ < ∞) while is a constant of a particular integral of the GMP and represents the mass variations of the primaries and the disk. The masses of the primaries vary in accordance with the CML while motion of the primaries is described by the GMP. Equations (3) and (4) are different from those of Bekov (1988), Luk’yanov (1989), Singh and Leke (2010, 2012, 2013), Leke and Orum (2024), and, Leke et al. (2025). When there is no disk, the equations of the dynamical systems in Bekov (1988), Luk’yanov (1989) can be recovered from equations (3) and (4). Next, the autonomized system with constant coefficients is given by (Leke and Singh 2023):

υ is the mass parameter, which is defined as the ratio of the mass of the smaller primary to the total mass of the primaries while M_d is the mass of the disk ρ_c and is the radial distance of the infinitesimal body in the autonomized R3BP with constant coefficients.

Locations Of Collinear Equilibrium Points

To obtain the CEPs, we need to solve equations (10) when ζ = 0 . In this premise, the system of equations (4) reduce to

When κ = 1, equations (13-15) are reduced analogously to those of the R3BP of Kushvah (2008), Singh and Taura (2013) [23], Singh and Leke (2014). Equations (13-15) are the autonomized equations of motion of the infinitesimal mass in the gravitational field of the primaries with masses κ (1−υξ) an=d − Kυ 1 when there is a disc of mass κM_d enclosing the system, while the two primaries are placed along the axis ξ −of the synodic frame of reference at the positions ξ₁ = −υ and ξ₂ = 1−υ such that the origin always coincides with the barycenter of the system. Therefore, as increases, the masses of the primaries and the disc also increase and vice versa.

Equations (13) admits the Jacobi integral

where C is the Jacobi or energy constant.

Collinear Equilibrium Points

The CEPs of the autonomized system are the solutions when and of equations (13) when the components of velocity and acceleration are zero. That is, we need to solve the equation

We denote equation (17) by f (ξ ), so that

The abscissas of the collinear points are the roots of (18). Now, observe that

Equation (20) is positive if

Therefore, the real roots of the equation (18) will lie in the open’s interval (−υ, −υ −1), (1−υ,0)and (2 −υ, 1−υ ). These roots correspond to the collinear points and we shall denote them appropriately when they are found.

Now, rewriting equation (18), we have

Now, observe that equation (22) can have the following forms depending on positions of the primaries. When ξ < −υ , −υ <ξ < 1−υ and ξ > 1−υ , then equation (22), give respectively

Hence, in order to investigate the position of collinear LPs, we divide the orbital plane Oxy into three parts with respect to the primaries: ξ < −υ , −υ <ξ < 1−υ and ξ > 1−υ .

The first case ξ < −υ , we have ξ +υ −1 < 0 and so In the second case, when −υ <ξ < 1−υ , we have ξ +υ > 0 and ξ +υ −1 < 0 . The third case ξ > 1−υ , implies that ξ +υ −1 > 0 and if this happens, then ξ +υ > 0 also follows. We shall consider these one after the other.

Case 1

In this region, the CEP lies in the open interval ξ ∈(− ∞, −υ ). In this case, if we take the derivative of equation (24) with respect to ξ , we see that P′(ξ ) > 0 and since , this means P(ξ ) is a strictly increasing function in this interval. Similarly, from equation of (23), we see that . It follows that Q(ξ ) is strictly decreasing function in the interval (− ∞, −υ ). Hence ; so there exist a unique point say ξ₁ in this open interval at which and we denote this point as 1 L .

Case 2

In this case, the CEP lies in the interval ξ ∈(υ, 1−υ ). Here, two cases arise: the first one is when ξ ∈(0,υ ) and when it lies in the interval ξ ∈(0, 1−υ ). We discuss them separately.
a. When ξ ∈(0,υ ),
In this case, we have from (18):

Now, for the chosen parametric values and, it is seen that f (0) < 0 and since , then there exists a point say ξ₂(0,υ ) ∈ at which f(ξ₂)=0. We denote this point by L₂. b. When ξ ∈(υ, 1−υ ); In this case, we suppose that T < (1−υ ) 2 and so, we consider the two cases: when

The existence of the points L₂₁ and L₂₂ depend only on the mass parameter υ, the mass of the disk d M and the parameter T which determine the density profile of the disk.

Case 3

In this region, we have ξ > 1−υ this means that the CEP will lie in the interval ξ ∈(1−υ,∞ ). In this case, if we take the derivative of equation (26) with respect to , we see evidently, that for 0 <υ ≤ 1/2 , we have P′(ξ > 0) and since . This means P(ξ ) that is strictly increasing in the interval. Also, because Q′(ξ ) > 0 , (1−υ )Q(ξ ) < 0 and , then it follows Q(ξ ) that is also a strictly increasing function for ξ ∈(1−υ,∞ ). Hence, . Hence, there exists a point say ξ₃ ∈(1−υ,∞) at which ( 0)

We call this point L₃ .

Thus, in addition to the usual three CEPs of the autonomized R3BP with variable masses, due to the presence of the disk and the mass parameter, there exist two additional CEPs L₂₁ and L₂₂ . Hence, we have the points L₁ and L₂ located to the left of the center of mass while L₂₁ , L₂₂ and L₃ are positioned at the right.

In order to make a comprehensive numerical analysis, we represent the CEPs respectively as:

Then, we substitute these in equations (24), (25) and (26), respectively, to get the respective polynomials given in equations (29- 32).

We shall explore the polynomials numerically to determine the roots ε_i(i =1,2,3) and ε_2i(j=1,2).

These roots give the locations of the CEPs L₁ , L₂ , and L_2j =(j=1,2), respectively away from the positions of the bigger and smaller primary. We shall determine the roots of the polynomials numerically.

The CEPs L_i(x_i, 0:i =1,2,3) and of the system of equations of the non-autonomous system (1) are obtained using the Mestschersky transformation in the form (Luk’yanov 1990):

where are the five CEPs of autonomized system (13). All the CEPs in this case are function of time, and so the locations of the CEPs of equations of motion with variable coefficients always change with time.

Stability osf Collinear Equilibrium Points

The stability of linear systems of ordinary differential equations with constant coefficients is determined completely by the Eigen values of the coefficient matrix. Due to the perturbations induced the location of the test particle would be displaced from the CEPs. If the ensuing motion of the infinitesimal is a rapid disappearance from the vicinity of the point, we can call such a location an “unstable one”, if conversely, the body only oscillates about the CEPs, it is said to be a “stable position” (in the sense of Lyapunov).

The variational equations of motion are

where u and v are small displacements in the position of the CEP (ξ₀, 0,0 ), respectively and the second order partial derivatives of are denoted by subscripts and the superscript indicates that the derivatives are to be computed at the CEPs.

Now, the corresponding characteristic equation of the variational equations (34) is

The second order partial derivatives will be evaluated at the IEPs At the collinear points we have , so that the partial derivatives at the CEPs are:

Now, the nature of the roots of the polynomial (37) will depend on the discriminant together with the variations in a and b which solely depends on the mass variation parameter κ , while the values are affected by the mass parameter and disc parameter. If one of the roots of the characteristic equation (37) is positive, then the point will be unstable. If however, the four roots are different pure imaginary roots, then the point will be a stable collinear point. We shall numerically examine this outcome.

In the case of the stability of the CEP L₃ of the non-autonomous system, it is seen that since this EP is time dependent, we use the definition of a Lyapunov stable solution (Krasnov et al. 1983) [24,25], to get

Equation (38) proves the instability of the solution according to the Lyapunov’s theorem. Numerical Exploration

The CEPs of R3BP when the binary is surrounded by a disc and undergo mass variation with time in accordance with the CML (1952) and motion takes place within the frame of the GMP, has been investigated. In this section, we devote our efforts to the numerical simulation of the obtained analytical results. All numerical explorations have been carried out with the help of the software package, Mathematica (Wolfram 2020). For our numerical explorations, we consider motion of a dust grain around the CEPs of a rapidly rotating binary star MAXI J1659-152, in the constellation Ophiuchus, which consists of a black hole having mass three times the mass of the Sun while its companion is a red dwarf star having mass 20% the Sun’s mass. In this case, since the mass variation constant lies in the interval 0 <κ < ∞ , consequently, for the binary star MAXI J1659-152, we will have . Now, for the autonomized system the mass variation parameter depicts the variations in the masses of the binary and the disc. Hence, the mass of disc varies as the mass variation constant κ varies.

Locations of Collinear Equilibrium Points

The equations (29-32) give the respective roots corresponding to the five CEPs. We shall numerically compute these locations for the binary in Table 1.

Table 1 is the numerical computations of the positions of the CEPs of a dust grain in the gravitational environment of MAXI J1659- 152 system. We classify mass loss or gain on the basis of the mass variation parameter. When is tending to zero, we assume that mass is lost. When it is moving away from zero, we assume that the stars gain mass. It is seen that when mass loss lies between 1×10 −9 ≤κ < 0.1, only four CEPs exist. However, as the binary and disc gain mass, in the interval 0.1 ≤κ < ∞ , five CEPs exist. Also, we see that as the binary and the disc gain mass, the position of the CEPs L₁ and L₂ move closer to the more massive star, while the point L₃ moves nearer to the less massive star. It is noted that 0.9 ≤κ < ∞ when the positions of the CEP L₃ , are no longer affected by the mass variations of the stars and the disc. In the case of the CEP L₁, L₂ and L₂₂ the positions remain unchanged as the masses and the disc gains mass ten times the initial mass. Below in Figure 1, we draw the locations of these points when the mass of the CBD is 0.00001.

Figure 1, gives the locations of all the five CEPs for the MAXI J1659-152 when the mass variation parameter is 0.1 and when the disc mass is 0.00001. It is seen that the location of the CEP L₁ is positioned on the left of the more massive star while the points L₂₁₂ , L and L₂₂ are located between the stars and the point L₃ lies to right of the less massive star.

Zero Velocity Curves and Poincare Surface of Sections

The zero velocity curves (ZVC) play important role in the qualitative investigation of motion of the infinitesimal mass and are the solution of equation (16) when the velocity components are zero. In other words, we have

where C is the Jacobi constant which determines the topology of the zero velocity.

The Poincaré surface of section (PSS) is the locus of plot of all points in the x, x − plane at the location where the orbit intersects the plane. The PSS provides insights into the stability and nature of orbits in the restricted problems. Stable orbits appear as clusters or continuous curves, while chaotic regions appear as scattered points. In order to plot the PSS of the dust grain, we require location (x, y) and the velocity coordinates (x, y) at different instant to give a 4-D dimensional phase space (x, y, x, y). With the help of the Jacobi integral (16), we express the variable in terms of the other three variables and then substitute it in the equations of motion for each occurrence of y , we reduce the phase space to a 3-D (x, y, x) subspace of the original 4-D. Now, we follow many orbits numerically and determine the values of the remaining variables (x, x) each time the orbit passes through the surface of section plane with positive y . The restriction y > 0 is needed so that structures on the plot don’t overlap. Next, in Table 2, we compute the numerical estimates of the energy constants for the ZVC and Poincare Surface of Section of CEPs in MAXI J1659-152 System for 0 <κ ≤100 and 0 < ≤ 0.01 M_d .

Table 2, provides the numerical estimates of the energy constants given in equation (39) for CEPs in MAXI J1659-152 System under effects of the mass variation and disc parameters. C₁, C₂, and C₃ corresponds to the energy constants of the CEPs L₁, L₂, and L₃, respectively. We use these values to plot the ZVC in Figure 2-5 and the PSS in Figure 6-9, for the CEPs L₁, L₂, and L₃, respectively

The ZVC around the CEPs L_1,2,3 and L₂₂ have been drawn for different values of the mass of the disc and energy constants, in Figure 2-5. The yellow, orange, green and brown points represents, the positions of the CEPs L₁, L₂, L₂₂ and L₃, respectively, while the black and red spots designates the bigger and smaller star, respectively. The blue region depicts region where motion is allowed while the red region indicates region where motion is forbidden. It is seen that as the mass of the disc increases, the region where motion of the infinitesimal is allowed reduces. In the case of the PSS, we used the method “EventLocator” and integrating orbits for very long period of time, to plot the PSS in Figures 6 to 9. These plots have been drawn for different mass of disc and the energy constant C given in Table 2. The red dots represent the intersections of the test particle’s trajectory with the PSS at specific conditions. It is seen that there are no visible islands noticed when the mass of the disc is very small as seen in Figure 6a, 7a, 8a and Figure 9a, as regions sparsely populated near the primaries suggest areas where the gravitational influence is strong enough to either capture or eject the particle, leading to unstable regions. The clustering of red dots in certain areas as seen in Figure 6, 7, 8 and Figure 9 suggests the presence of stable or quasi-periodic orbits. These regions correspond to regular motion where the particle revisits similar states in phase space. Hence, the region of the orbital stability is controlled by the energy constant which is defined by the disc and mass variation parameters.

Stability of Collinear Equilibrium Points

The characteristic equation corresponding to the CEP L₃ has been obtained in equation (37). In this subsection of the paper, we shall numerically determine the four roots of equation (37) corresponding to each CEPs in Table 3-6, when the dust grain particle is gravitating under the influence of the binary and the CBD. From the numerical estimates of the four roots, it is seen that, there is at least one positive root which causes instability at the location of the CEPs. This outcome is different from that of Leke and Orum (2024) where the point can be stable. Same numerical procedures show that there exists at least a positive root in the cases of the other collinear equilibrium points.

Conclusion

• The paper investigated dynamics of a dust grain around CEPs of the circular R3BP, with variable masses and disc. The masses of the binary and the disc vary with time as defined by the combined Mestschersky law (CML) and the motion of the binary is determined by the Gylden-Mestschersky problem (GMP), with the assumption that the system is enclosed by a variable mass accretion disc of dust. The equations of motion of the nonautonomous and the autonomized systems were deduced and the CEPs were investigated. It was observed that there exist four or five CEPs depending jointly on the mass variation parameter and mass of the disc. These points were defined by the mass parameter , mass of disc and the mass variation constant . Further, the stability of the CEPs of the autonomized and non-autonomous systems was investigated and it was seen that these points are unstable due to a positive root and time, respectively. Numerical illustrations were rendered for a passively gravitating dust particle in the gravitational environment of the binary MAXI J1659-152. The positions, ZVC, PSS and stability are analyzed numerically and the impacts of the perturbing forces are shown. We conclude that the presence of the accreting disc in the presence of mass variations of the stars does not automatically result in the existence of two additional CEPs and at some points, the locations of the CEPs remain unchanged and are not influenced by the mass gain of the binary and disc. This model can be used to model the long-term motion of satellites and planets in binary systems with gradual accretion of disk of dust in the environment.

References

1. Szebehely VG (1967) Theory of Orbits. Academic Press, New York.
2. Gylden H (1884) Die Bahnbewegungen in EinemSysteme von zweiKörpern in demFalle, dass die MassenVerNderun- Gen Unterworfen Sind, AstronomischeNachrichten 109(1): 1-6.
3. Mestschersky IV (1902) Ueber die Integration der Bewegungs- gleichungenimProbleme zweierKörper von vernderli- cher Masse, AstronomischeNachrichten 159(15): 229-242.
4. Mestschersky IV (1952) Works on the mechanics of bodies of variable mass, GITTL, Moscow p: 205.
5. Gelf’gat BE (1973) Current Problems of Celestial Mechanics and Astrodynamics, Nauka, Moscow.
6. Bekov AA (1988) Libration points of the restricted problem of Three Bodies with variable Mass. Soviet Astronomy Journal 33: 92-95.
7. Luk’yanov LG (1989) Particular solutions in the restricted problem of three bodies with variable masses. Astronomical Journal of Academy of Sciences of USSR 66: 180-187.
8. Luk’yanov LG (1990) Stability of libration points in the restricted three-body problem of variable masses. Soviet Astronomical Journal 67: 167-172.
9. Singh J, Leke O (2010) Stability of the photogravitational restricted three-body problem with variable masses. Astrophysics and Space Science 326: 305- 314.
10. Singh J, Leke O (2012) Equilibrium points and stability in the restricted three- body problem with oblateness and variable masses. Astrophysics and Space Science 340(1): 27-41.
11. Singh J, Leke O (2013) Effects of oblateness, perturbation, radiation and varying masses on the stability of equilibrium points in the restricted three-body problem. Astrophysics and Space Science 344(1): 51-61.
12. Gao F, Wang Y (2020) Approximate analytical periodic solutions to the restricted three-body problem with perturbation, oblateness, radiation and varying Mass. Universe 6(8): 110.
13. Taura JJ, Leke O (2022) Derivation of the dynamical equations of motion of the R3BP with variable masses and disk. FUDMA Journal of Sciences 6: 125-s133.
14. Leke O, Singh J (2023) “Out-of-plane equilibrium points of extra-solar planets in the central binaries PSR B1620-26 and Kepler-16 with cluster of material points and variable masses. New Astronomy Elsevier 99: 101958.
15. Leke O, Orum SA (2024) Motion and zero velocity curves of a dust grain around collinear libration points for the binary IRAS 11472-0800 and G29-38 with a triaxial star and variable masses. New Astronomy 108: 102177.
16. Begelman MC, Blandford RD, Rees MJ (1980) Massive black hole binaries in active galactic nuclei. Nature 287: 307.
17. Doyle Laurance R, Carter Joshua A, Fabrycky Daniel C, Slawson Robert W, Howell Steve B, et al. (2011) Kepler-16: A transiting circumbinary planet. Science 333(6049): 1602.  
18. Heath RM, Nixon CJ (2020) On the orbital evolution of binaries with circumbinary discs A&A 641: A64.
19. Tiede C, Zrake J, MacFadyen A, Haiman Z (2020) Gas-driven inspiral of binaries in thin accretion disc. APJ 900(1): 43.
20. Kushvah BS (2008) Linear Stability of equilibrium points in the generalized photogravitational chermnykh’s Problem. Astrophy. Space Sci 318: 41-50.
21. Singh J, Leke O (2014) Analytic and numerical treatment of motion of dust grain particles around triangular equilibrium points with post-AGB binary star and disc. Advances in Space Research 54(8): 1659-1677.
22. Vincent AE, Perdiou AE, 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 Frontiers in Astron. and Space Sci 9: 872459.
23. Singh J, Taura JJ, (2013) Motion in the generalized restricted three-body problem. Astrophys. Space Sci 343(1): 95-106.
24. Krasnov ML, Kiselyov AI, Makarenko GI (1983) A Book of Problems in Ordinary Differential Equations, MIR Publications, Moscow p: 255-291.
25. Wolfram S (2020) The Mathematica Book.5th Edition. Wolfram Media, Champaign pp: 1488.