Investigation of the effect of albedo and oblateness on the circular restricted four variable bodies problem

Few-body problem attract many scientists for a long time in celestial mechanics and dynamical astronomy. The restricted three-body problem and four-body problem with many perturbations like different shapes of the primaries, resonance, variable mass of the primaries as well as infinitesimal body, Coriolis and centrifugal forces, Pointing-Robertson drag, solar radiation pressure and albedo effects etc., have been studied by many scientists. Simó [36] investigated the linear stability of the relative Lagrangian solutions in the four-body problem. Hadjidemetriou [17]studied the periodic orbits of fourth body in the restricted four-body problem with respect to rotating frame. Michalodmitrakis [27] generalized the restricted three-body problem to the restricted four-body problem and studied about the equilibrium points, regions of possible motion and periodic orbits. Kalvouridis et al. in the series of three papers [20,21 and 22] investigated the equilibrium points and stability in the restricted four-body problem under the effect of oblateness and radiation pressure. And also performed the zero-velocity surface and curves. Baltagiannis et al. [14] shown that equilibrium points depend on the mass of the primaries in the restricted four-body problem. Papadauris et al. [31,32] investigated the existence, locations, stability and periodic orbits of the equilibrium points in and out of the orbital plane in the photo-gravitational circular restricted four-body problem. Ansari [1] studied the periodic orbits in the restricted four-body problem around lagrangian points in three cases. In the first case, he has considered all three primaries as spherical in shape. In the second case, he has taken one of the three primaries as an oblate body. And in the third and last case, he has taken two of the primaries as oblate body and all the three primaries are source of radiation pressure. Falaye [16] investigated the stability of the equilibrium points in the restricted four-body problem under the effects of oblateness and solar radiation pressure and found that these equilibrium points are unstable. Arribas et al. [22] investigated the equilibria of the symmetric collinear restricted four-body problem where primaries are placed in a collinear central configuration with both masses and radiation pressure of the peripheral bodies are equal. Papadakis [30] performed the 21 families of simple 3D symmetric periodic orbits as well as the typical orbits of all symmetry type 3D orbits in the circular restricted four-body problem. After examined the stability, he illustrated the characteristic curves and stability diagrams of families of 3D periodic orbits. Asique et al. [10][11][12][13] studied the restricted four-body problem with different shapes of the primaries with solar radiation pressure. They have placed one of the primaries at the lagrangian points of the classical restricted three-body problem. They have illustrated the equilibrium points and zero-velocity curves for these models. Singh et al. [45,46] investigated in and out of plane equilibrium points in the circular restricted four-body problem with the effect of solar radiation pressure.

On the other-hand many scientists have studied about the albedo on these models. Anselmo et al. [8] performed the periodic perturbations of the satellite is the radiation pressure due to the sun-light reflected by the Earth. Rocco [34] evaluated the terrestrial albedo by using earth albedo model and orbital dynamics model and also calculated the irradiance at the satellite with radiant flux from each cells. Idrisi [18] investigated the existence and stability of the circular restricted three-body problem under the effect of albedo when smaller primary is an ellipsoid. Many scientists have investigated on these models with variable masses as Jeans [19], Meshcherskii [26], Shrivastava et al. [35], Lichtenegger [24], Singh et al. [37][38][39][40][41], Lukyanov [25], Zhang et al. [44], Abouelmagd and Mostafa [2], Mittal et al. [28], Ansari [3][4][5][6], etc. And also many scientists have explained the basins of attraction in these models as Douskos [15], Kumari and Kushvah [23], Paricio [33], Ansari [5,6], Zotos [45][46][47][48], etc.

This paper deals the motion of the infinitesimal body in the circular restricted four-body problem in which the masses of the primaries as well as the mass of the infinitesimal body vary and also one of the primaries is taken as source of radiation pressure due to which albedo effects are produced and another primary is considered as oblate body. We have studied this problem in various sections. In the statement of the problem and equations of motion section, we have constructed the equations of motion of the infinitesimal variable mass under the effect of these perturbations where the third primary is fixed at the triangular equilibrium point of the classical restricted three-body problem and also the variation of Jacobi Integral constant has been determined. In the computational work section, we have plotted numerically the equilibrium points, Poincaré surface of sections and basins of attraction in five cases by using Mathematica software. In the stability section, we have examined the stability of the equilibrium points under the effect of the perturbations. Finally, we have concluded the problem.

Letm₁(t), m₂(t), m₃(t) and m(t) be four variable masses (m₁(t) > m₂(t) >> m₃(t)). The configuration of the m₁(t), m₂(t) and m(t) is taken as the restricted three-body problem. The third primary m₃(t) is considered as so small that it is not influencing the motion of the primaries m₁(t) and m₂(t). Thus the center of the rotation of the system remains the same as for the bodies m₁(t) and m₂(t). The body m₃(t) is placed at one of the triangular equilibrium points of the R3BP. m(t) is moving in the space under the influence of these three primaries but not influencing them. In this way these four variable bodies formed the restricted circular four-body problem.

Here we have considered m₃(t) as oblate body with oblateness factor σ (the shape of the oblate body remain same though the mass is variable) and m₁(t) is a source of radiations(F_p) due to these radiations, m₂(t) and m₃(t) produce albedo F_A and F_a respectively (i.e. Albedo = (radiation reflected back into the space)/(incident radiation)). The infinitesimal body m(t) is moving in the space under the gravitational forces F₁, F₂ and F₃ of the primaries respectively. The total force on infinitesimal body will be F = F₁(1−p₁) + F₂(1−p₂) + F₃(1−p₃). Where 0≤p₁=F_p/F₁ < 1,0 ≤ p₂ = F_A/F₂ < p₁, 0 ≤ p₃ = F_a/F₃ < p₁. The line joining m₁(t) and m₂(t) is taken as x-axis and are moving around their common center of mass which is taken as origin O and the line perpendicular to the x-axis and passing through the origin is considered as y-axis, the line passing through the origin and perpendicular to the plane of motion of the primaries is taken as z-axis. Let us consider the synodic coordinate system, initially coincide with the inertial coordinate system, and revolving with angular velocity ω(t) about z-axis. Let the coordinates of m₁(t), m₂(t), m₃(t) and m(t) in the rotating frame be (x₁, 0, 0), (x₂, 0, 0), (x₃, y₃, 0) and (x, y, z) respectively (Fig. 1). Following the procedure given in Ansari [7], we can write the equations of motion of the infinitesimal variable mass in the rotating coordinate system when the variation of mass is non-isotropic with zero momentum as

Fig. 1

{m˙(t)m(t)(x˙−ω(t) y)+(x¨−ω˙(t) y−2ω(t) y˙−ω2(t)x)=−μ1(t)(x−x1)(1−p1)r13−μ2(t)(x−x2)(1−p2)r23−μ3(t)(x−x3)(1−p3)r33−3μ3(t)(x−x3)(1−p3)σ2r35,m˙(t)m(t)(y˙+ω(t) x)+(y¨+ω˙(t) x+2ω(t) x˙−ω2(t) y)=−μ1(t)y(1−p1)r13−μ2(t)y(1−p2)r23−μ3(t)y(1−p3)r33−3μ3(t)(y−y3)(1−p3)σ2r35,m˙(t)m(t)z˙+z¨=−μ1(t)z(1−p1)r13−μ2(t)z(1−p2)r23−μ3(t)z(1−p3)r33−3μ3(t)z(1−p3)σ2r35.$$ \begin{equation} \left\{\begin{array}{c} \frac{\dot{m}(t)}{m(t)} (\dot{x}-\omega(t)~ y)+(\ddot{x}-\dot{\omega}(t) ~y-2\omega(t) ~\dot{y}-\omega ^{2}(t) x) \\ \\ = -\frac{\mu _{1}(t) (x-x_{1} )(1-p_1)}{r_{1}^{3}} -\frac{\mu _{2}(t) (x-x_{2} )(1-p_2)}{r_{2}^{3}} -\frac{\mu _{3}(t) (x-x_{3} )(1-p_3)}{r_{3}^{3}} \\ \\-\frac{3\mu _{3}(t) (x-x_{3} )(1-p_3)\sigma}{2r_{3}^{5}},\\ \\ \frac{\dot{m}(t)}{m(t)} (\dot{y}+\omega(t) ~x)+(\ddot{y}+\dot{\omega}(t)~x + 2\omega(t)~\dot{x}-\omega ^{2}(t) ~y)\\ \\ = -\frac{\mu _{1}(t) y (1-p_1)}{r_{1}^{3}} -\frac{\mu _{2}(t) y (1-p_2)}{r_{2}^{3}}-\frac{\mu _{3}(t) y (1-p_3)}{r_{3}^{3}} \\ \\-\frac{3\mu _{3}(t) (y-y_{3} )(1-p_3)\sigma}{2r_{3}^{5}},\\ \\ \frac{\dot{m}(t)}{m(t)} \dot{z}+\ddot{z} =-\frac{\mu _{1}(t) z(1-p_1)}{r_{1}^{3}} -\frac{\mu _{2}(t) z(1-p_2)}{r_{2}^{3}} -\frac{\mu _{3}(t) z(1-p_3)}{r_{3}^{3}} \\ \\-\frac{3\mu _{3}(t) z(1-p_3)\sigma}{2r_{3}^{5}}. \end{array} \right. \end{equation} $$(1)

where, ri2=(x−xi)2+(y−yi)2+z2 $r_{i}^{2} =(x-x_{i} )^{2} +(y-y_i)^{2} +z^{2}$ are the distances from the primaries to the infinitesimal body in the rotating coordinate system respectively, μ_i(t) = G m_i(t), (i = 1, 2, 3).

Using Meshcherskii transformation

x=ξR(t),y=ηR(t),z=ζR(t),dtdτ=R2(t),ri=ρiR(t),ω(t)=ω0R2(t),xi=ξiR(t),yi=ηiR(t),μ(t)=μ1(t)+μ2(t)=μ0R(t),μi(t)=μi0R(t),m(t)=m0R(t),R(t)=αt2+2βt+γ $x=\xi R(t),\; y=\eta R(t),\; z=\zeta R(t),\frac{dt}{d\tau} =R^{2} (t),r_{i} =\rho _{i} R(t), \omega (t)=\frac{\omega _{0}}{R^{2} (t)} ,x_{i} =\xi _{i} R(t),y_{i} =\eta _{i} R(t),\mu(t)=\mu _{1}(t)+\mu _{2} (t)=\frac{\mu _{0}}{R(t)} , \mu _{i} (t)=\frac{\mu _{i0}}{R(t)},m(t) =\frac{m_{0}}{R(t)} ,R(t)=\sqrt{\alpha t^{2} +2\beta t+\gamma}$ ,

where α, β, γ, μ₀, μ_i0, m₀ are constants for i = 1, 2, 3.

The system (1) becomes

{ξ″−2ω0η′−(αt+β)ξ′=∂W∂ξ,η″+2ω0ξ′−(αt+β)η′=∂W∂η,ζ″−(αt+β)ζ′=∂W∂ζ.$$\begin{equation} \left\{\begin{array}{c} \xi''-2\omega _{0}\eta'-(\alpha t + \beta)\xi'=\frac{\partial W}{\partial \xi},\\ \\ \eta''+2\omega _{0}\xi'-(\alpha t + \beta)\eta'=\frac{\partial W}{\partial \eta }, \\ \\ \zeta''-(\alpha t + \beta)\zeta'=\frac{\partial W }{\partial\zeta }. \end{array}\right. \end{equation}$$(2)

where,

W=12((αt+β)2+ω02−(αγ−β2))(ξ2+η2)+12((αt+β)2−(αγ−β2))ζ2−(αt+β)ξη+μ10(1−p1)ρ1+μ20(1−p2)ρ2+μ30(1−p3)ρ3+μ30(1−p3)σ2ρ33,ρi2=(ξ−ξi)2+(η−ηi)2+ζ2.$$ \begin{equation}\begin{array}{l} W =\frac{1}{2} ((\alpha t + \beta)^{2} +\omega _{0}^{2} -(\alpha \gamma-\beta^2) )(\xi^{2} +\eta^{2})+\frac{1}{2} ((\alpha t + \beta)^{2} -(\alpha \gamma-\beta^2) )\zeta^{2}-(\alpha t + \beta)\xi\eta \\ + \frac{\mu _{10}(1-p_1)}{\rho _{1}} +\frac{\mu _{20} (1-p_2)}{\rho _{2}} +\frac{\mu _{30} (1-p_3)}{\rho _{3}} +\frac{\mu _{30} (1-p_3)\sigma}{2\rho _{3}^{3}}, \\ \rho _{i}^{2} =(\xi -\xi _{i} )^{2} +(\eta-\eta_i)^{2} +\zeta^{2}. \end{array} \end{equation} $$

Prime (′) is the differentiation w.r.to τ. Taking unit of mass, distance and time at initial time t₀ such that

μ0=1, ξ1+ξ2=1, G=1, ω0=1, αt0+β=α1(constant).$$\mu _{0} =1,~\xi_1+\xi_2 = 1,~G = 1,~\omega _{0} = 1,~\alpha t_{0} +\beta=\alpha _{1} (constant).$$

Introducing the new mass parameter as μ₁₀ = 1−υ, μ₂₀ = υ, μ₃₀ = α₂υ, where α₂ << 1 and υ is the ratio of the mass of the primary m₂ to the total mass of the primaries m₁ and m₂. Finally, the system (2) becomes

{ξ″−2η′−α1ξ′=∂U∂ξ,η″+2ξ′−α1η′=∂U∂η,ζ″−α1ζ′=∂U∂ζ.$$ \begin{equation} \left\{\begin{array}{c} \xi''-2\eta'-\alpha _{1} \xi'=\frac{\partial U}{\partial \xi}, \\ \\ \eta''+2\xi'-\alpha _{1} \eta'=\frac{\partial U}{\partial \eta}, \\ \\ \zeta''-\alpha _{1} \zeta'=\frac{\partial U}{\partial \zeta}. \end{array} \right. \end{equation} $$(3)

where, U=12(α12+k)(ξ2+η2+ζ2)−12ζ2−α1ξη+(1−υ)(1−p1)ρ1+υ(1−p2)ρ2+α2υ(1−p3)ρ3+α2υ(1−p3)σ2ρ33 $U =\frac{1}{2} (\alpha _{1}^{2}+k )(\xi^{2} +\eta^{2} +\zeta^{2})-\frac{1}{2}\zeta^{2}-\alpha _{1}\xi\eta + \frac{(1-\upsilon)(1-p_1)}{\rho _{1}} +\frac{\upsilon (1-p_2)}{\rho _{2}}+\frac{\alpha_2 \upsilon (1-p_3)}{\rho _{3}}+\frac{\alpha_2 \upsilon (1-p_3)\sigma}{2\rho _{3}^{3}}$ , ρi2=(ξ−ξi)2+(η−ηi)2+ζ2,αγ−β2=1−k $\rho_{i}^{2} =(\xi -\xi _{i})^{2} +(\eta-\eta_i)^{2} +\zeta^{2}, \alpha \gamma-\beta^{2}=1-k$ and (ξ₁, η₁) = (υ, 0),(ξ₂, η₂) = (1−υ, 0), (ξ3,η3)=(υ−12,32) $(\xi _{3},\eta_3) =(\upsilon-\frac{1}{2},\frac{\sqrt{3}}{2})$ .

If there are constant masses then there is a constant motion (i.e Moulton [29]), the Jacobi Integral Constant defined as

JIC=2U−2((ξ′)2+(η′)2+(ζ′)2),$$ \begin{equation} J_{IC}=2U-2((\xi')^{2}+(\eta')^{2}+(\zeta')^{2}), \end{equation} $$(4)

Multiplying in the first equation of (3) by ξ′, in the second equation of (3) by η′ and in the third equation of (3) by ζ′ and add and using equation (4), we get the variation of the Jacobi Integral Constant as

dJICdτ=−2α1((ξ′)2+(η′)2+(ζ′)2),$$ \begin{equation} \frac{dJ_{IC}}{d\tau}=-2\alpha_{1}((\xi')^{2}+(\eta')^{2}+(\zeta')^{2}), \end{equation} $$(5)

Where J_IC is the Jacobi Integral Constant.

In this section, we have drawn the locations of equilibrium points, the Poincaré surfaces of section and the basins of attraction for five different cases by using Mathematica software (in all the calculations υ = 0.019, α₂ = 0.01):

Third primary is placed at one of the Lagrangian points of the classical restricted three-body problem

(i.e. k = 1, α₁ = 0, p₁ = 0, p₂ = 0, p₃ = 0, σ = 0),

3.1

Equilibrium points during in-plane and out of plane motions

The solutions of U_ξ = 0, U_η = 0, and U_ζ = 0, will be the location of equilibrium points but the solutions of these equations represent the location of equilibrium points during in-plane motions when (ξ≠0, η≠0, ζ = 0), (Fig. 2) and represent the location of equilibrium points of the out of plane when (ξ≠0, η = 0, ζ≠0), (Fig. 3) and (ξ = 0, η≠0, ζ≠0), (Fig. 4). At the equilibrium points, all the derivatives of the co-ordinates with respect to the time will be zero. Therefore the infinitesimal body will stay at these points with zero velocity. Hence the singularities are known as stationary points. Where

Fig. 2

Uξ=(α12+k)ξ−α1η−(1−υ)(ξ−ξ1)(1−p1)ρ13−υ(ξ−ξ2)(1−p2)ρ23−υ3(ξ−ξ3)(1−p3)ρ33−3α2υ(ξ−ξ3)(1−p3)σ2ρ35,$$ \begin{equation} U_{\xi}=(\alpha _{1}^{2} +k)\xi -\alpha _{1}\eta -\frac{(1-\upsilon )(\xi -\xi_1 )(1-p_1)}{\rho _{1}^{3}}-\frac{\upsilon (\xi -\xi_2)(1-p_2)}{\rho_{2}^{3}}-\frac{\upsilon _{3} (\xi-\xi_{3} )(1-p_3)}{\rho_{3}^{3}}-\frac{3\alpha_2 \upsilon(\xi-\xi_{3} )(1-p_3)\sigma}{2\rho_{3}^{5}}, \end{equation} $$(6)

Uη=(α12+k)η−α1ξ−(1−p1)(1−υ)ηρ13−(1−p2)υηρ23−υ3(η−η3)(1−p3)ρ33−3α2υη(1−p3)σ2ρ35,$$ \begin{equation}U _{\eta} = (\alpha _{1}^{2} +k)\eta -\alpha _{1} \xi -\frac{(1-p_1)(1-\upsilon )\eta}{\rho _{1}^{3}}-\frac{(1-p_2)\upsilon \eta}{\rho _{2}^{3}}-\frac{\upsilon_{3}(\eta-\eta_{3} )(1-p_3)}{\rho_{3}^{3}}-\frac{3\alpha_2 \upsilon \eta (1-p_3)\sigma}{2\rho_{3}^{5}}, \end{equation} $$(7)

Uζ=(α12+k−1)ζ−(1−p1)(1−υ)ζρ13−(1−p2)υζρ23−υ3ζ(1−p3)ρ33−3α2υζ(1−p3)σ2ρ35.$$ \begin{equation}U_{\zeta} = (\alpha _{1}^{2} +k-1)\zeta -\frac{(1-p_1)(1-\upsilon )\zeta}{\rho _{1}^{3}} -\frac{(1-p_2)\upsilon \zeta}{\rho _{2}^{3}}-\frac{\upsilon_{3}\zeta (1-p_3)}{\rho_{3}^{3}}-\frac{3\alpha_2 \upsilon \zeta (1-p_3)\sigma}{2\rho_{3}^{5}}. \end{equation} $$(8)

3.1.1

During in-plane motion (i.e. ξ ≠ 0, η ≠ 0, ζ = 0) locations of lagrangian points

During in-plane motion, we have plotted graphs for the locations of the lagrangian points in five cases. We found six lagrangian points in which three points (L₁, L₂, L₃) are collinear and rest three points (L₄, L₅, L₆) are non-collinear when third primary is placed at one of the lagrangian points of the classical restricted circular three-body problem (Fig. 2(a) points with blue color). It is observed that the lagrangian points (L₁, L₄, L₅, L₆) and (L₂, L₃) lie left side and right side of the origin respectively. In the variable mass case, we found six lagrangian points (Fig. 2(b) points with green color) in which points (L₂, L₃, L₆) and (L₁, L₄, L₅) lie left side and right side of the origin respectively. Here (L₃, L₄) and (L₅, L₆) are symmetrical and rest points are non-symmetrical. In the rest three cases, we found eight lagrangian points with slight different locations (Fig. 2(c) points with red color, Fig. 2(d) points with black color, Fig. 2(e) points with magenta color). In all the figures black stars denote the locations of the primaries respectively.

3.1.2

During out of plane (i.e. ξ ≠ 0, η = 0, ζ ≠ 0 and ξ = 0, η ≠ 0, ζ ≠ 0) the locations of lagrangian points

During out of plane (i.e. ξ ≠ 0, η = 0, ζ ≠ 0), we found three lagrangian points on the ξ-axis in all the five cases (Fig. 3 (a) points in blue color, Fig. 3(b) points in green color, Fig. 3(c) points with red color, Fig. 3(d) points with black color, Fig. 3(e) points with magenta color). In all the cases, the points L₁ and (L₂, L₃) lie left and right side of the origin with slight different locations respectively. The black stars denote the locations of the primaries m₁ and m₂.

Fig. 3

On the other hand, during the out of plane (i.e. ξ = 0, η ≠ 0, ζ ≠ 0), we found three lagrangian points (Fig. 4) in all five cases. The points L₁ and L₃ lie left and right side of the origin respectively but the point L₂ lies at the origin. It is also observed from the Figure 4 that lagrangian points are moving away from the origin from the first case to the variable mass case and then toward the origin in the other cases.

Fig. 4

3.2

Poincaré Surface of Section

We also have drawn the Poincaré surface of sections for five cases in both the (ξ−ξ′)−plane (Fig.5 (i)) and the (η−η′)−plane (Fig.5 (ii), Fig. 5 (iii)). It is observed that in the ξ−ξ′-plane, the surfaces are shrinking and in the η−η′-plane, the surfaces first expanding and then shrinking.

Fig. 5

3.3

Basins of attraction

In this section, we have drawn the basins of attraction for the circular restricted four variable bodies problem in which we have taken one primary as solar radiation pressure due to which other two primaries produced albedo and also one of the primaries as an oblate body, by using Newton-Raphson iterative method for the five cases (i. Third primary is placed at one of the lagrangian points of the restricted three-body problem, ii. Variation of mass case, iii. Solar radiation pressure effect, iv. Albedo effect, v. Oblateness effect). The algorithm of our problem is given by

{ξn+1=ξn−(UξUηη−UηUξηUξξUηη−UξηUηξ)(ξn,ηn),ηn+1=ηn−(UηUξξ−UξUηξUξξUηη−UξηUηξ)(ξn,ηn).$$ \begin{equation} \left\{\begin{array}{c} \xi _{n+1} =\xi _{n} -\left(\frac{U _{\xi} U _{\eta \eta} -U _{\eta} U _{\xi \eta}}{U _{\xi \xi} U_{\eta \eta} -U _{\xi \eta} U _{\eta \xi}} \right)_{(\xi _{n} ,\eta _{n} )} ,\\ \eta _{n+1} =\eta _{n} -\left(\frac{U _{\eta} U _{\xi \xi} -U _{\xi} U _{\eta \xi}}{U _{\xi \xi} U _{\eta \eta} -U _{\xi \eta} U _{\eta \xi}} \right)_{(\xi _{n} ,\eta _{n} )} . \end{array}\right. \end{equation} $$(9)

Where ξ_n, η_n are the values of ξ and η coordinates of the n^th step of the Newton-Raphson iterative process. If the initial point converges rapidly to one of the lagrangian points then this point (ξ, η) is a member of the basin of attraction of the root. This process stops when the successive approximation converges to an attractor. We used color code for the classification of the lagrangian points on the (ξ, η)−plane. In the first case (Fig. 6), L₁, L₂, L₃ represent cyan color regions, L₄, L₅, L₆ represents light blue color regions. The basins of attraction corresponding to the lagrangian points L₁, L₂, L₃, L₄, L₅, L₆ extend to infinity. In the variable mass case (Fig. 7), L₁ represents cyan color region, L₂, L₃, L₅ represent light blue color regions, L₄ and L₆ represent light green color regions. The basins of attraction corresponding to the lagrangian point L₁ cover finite area but corresponding to the lagrangian points L₂, L₃, L₄, L₅, L₆ extend to infinity. While in the rest three cases (Fig. 8, Fig. 9, Fig. 10), L₁ represents cyan color region, L₂, L₃ represent light blue, L₅ represents blue color region and L₄, L₆, L₇, L₈ represent light green color regions. The basins of attraction corresponding to the all lagrangian points L₂, L₃, L₄, L₅, L₆, L₇, L₈ extend to infinity except L₁ which covers finite area. In this way a complete view of the basin structures created by the attractors. We can observe in detail from the zoomed part of all the figures in Fig. 6(b), Fig. 7(b), Fig. 8(b), Fig. 9(b), Fig. 10(b). The black points and black stars denote the location of the lagrangian points and the primaries respectively.

Fig. 6

We can examine the stability of the equilibrium points under the effect of Albedo and oblateness when all the masses are varying, by taking ξ = ξ₀ + ξ_d, η = η₀+ η_d, ζ = ζ₀ + ζ_d in system (3), we get

{ξd″−2ηd′−α1ξd′=ξdUξξ0+ηdUξη0+ζdUξζ0ηd″+2ξd′−α1ηd′=ξdUηξ0+ηdUηη0+ζdUηζ0ζd″−α1ζd′=ξdUζξ0+ηdUζη0+ζdUζζ0$$\begin{equation} \left\{\begin{array}{ccc} \xi_d''-2\; \eta_d'-\alpha _{1} \xi_d'&=&\xi_d~U _{\xi \xi }^{0} +\eta_d~U_{\xi \eta }^{0} +\zeta_d~U_{\xi \zeta }^{0} \\ \eta_d''+2\; \xi_d'-\alpha _{1} \eta_d'&=&\xi_d~U_{\eta \xi }^{0} +\eta_d~U_{\eta \eta }^{0} +\zeta_d~U_{\eta \zeta }^{0} \\ \zeta_d''-\alpha _{1} \zeta_d'&=&\xi_d~U_{\zeta \xi }^{0} +\eta_d~U_{\zeta \eta }^{0} +\zeta_d~ U_{\zeta \zeta }^{0} \end{array}\right. \end{equation}$$(10)

Where ξ_d, η_d and ζ_d are the small displacements of the infinitesimal body from the equilibrium point. The superscript zero denotes the value at the equilibrium point.

To solve system (10), let ξ_d = C₁e^λτ, η_d = C₂e^λτ, ζ_d = C₃e^λτ, where C₁, C₂ and C₃ are constant parameters. Substituting these values in system (10) and rearranging, we get:

{C1(λ2−α1λ−Uξξ0)−C2(2λ+Uξη0)−C3Uξζ0=0,C1(2λ−Uηξ0)+C2(λ2−α1λ−Uηη0)−C3Uηζ0=0,−C1Uζξ0−C2Uζη0+C3(λ2−α1λ−Uζζ0)=0.$$\begin{equation} \left\{\begin{array}{ccc} C_1(\lambda ^{2} -\alpha _{1} \lambda - U_{\xi \xi }^{0})-C_2(2\lambda + U_{\xi\eta }^{0} )-C_3 U_{\xi \zeta }^{0} =0,\\ C_1(2\lambda - U_{\eta \xi }^{0} )+C_2(\lambda ^{2}-\alpha _{1}\lambda - U_{\eta\eta }^{0} )-C_3 U_{\eta\zeta }^{0} =0,\\ -C_1 U_{\zeta \xi }^{0}-C_2 U_{\zeta\eta }^{0} +C_3(\lambda ^{2} -\alpha _{1} \lambda - U_{\zeta \zeta }^{0} )=0. \end{array}\right. \end{equation}$$(11)

The system (11), will have a non-trivial solution for C₁, C₂ and C₃ if

|λ2−α1λ−Uξξ0−(2λ+Uξη0)−Uξζ02λ−Uηξ0λ2−α1λ−Uηη0−Uηζ0−Uζξ0−Uζη0λ2−α1λ−Uζζ0|=0,$$\left|\begin{array}{ccc} {\lambda ^{2} -\alpha _{1} \lambda - U_{\xi \xi}^{0}} & {-(2\lambda + U_{\xi \eta}^{0} )} & {- U _{\xi \zeta}^{0}} \\ {2\; \lambda - U_{\eta \xi}^{0}} & {\lambda ^{2} -\alpha _{1} \lambda- U_{\eta\eta}^{0}}&{- U_{\eta \zeta}^{0}} \\ {- U_{\zeta\xi}^{0}} & {- U_{\zeta\eta}^{0}} & {\lambda ^{2} -\alpha_{1} \lambda - U_{\zeta\zeta}^{0}} \end{array}\right|=0,$$

which is equivalent to

λ6−3α1λ5+λ4(4+3α12−Uξξ0−Uηη0−Uζζ0)+α1λ3(−4−α12+2 Uξξ0+2 Uηη0+2 Uζζ0)+λ2(−(Uξη0)2−(Uξζ0)2+Uξξ0 Uηη0−(Uηζ0)2−4Uζζ0+Uξξ0Uζζ0+Uζζ0 Uηη0−α12 Uξξ0−α12 Uηη0−α12 Uζζ0)+α1λ((Uξη0)2+(Uξζ0)2−Uξξ0 Uηη0+(Uηζ0)2−Uξξ0 Uζζ0−Uηη0 Uζζ0)+(Uξζ0)2 Uηη0−2 Uξη0 Uξζ0 Uηζ0+Uξξ0 (Uηζ0)2+(Uξη0)2 Uζζ0−Uξξ0 Uηη0 Uζζ0=0,$$ \begin{equation} \begin{array}{l} {\lambda ^{6} -3\alpha _{1} \lambda ^{5} +\lambda ^{4} (4 +3\alpha _{1}^{2} - U_{\xi \xi} ^{0} - U_{\eta \eta}^{0} - U_{\zeta \zeta}^{0} )} {+\alpha _{1} \lambda ^{3} (-4 -\alpha _{1}^{2} +2~ U_{\xi \xi}^{0} +2~ U_{\eta \eta}^{0} +2~ U_{\zeta \zeta}^{0} )}\\\\ {+\lambda ^{2} (-(U_{\xi \eta} ^{0})^{2} -(U_{\xi \zeta} ^{0})^{2} + U _{\xi \xi}^{0} ~ U _{\eta \eta}^{0} -(U_{\eta \zeta}^{0})^{2} -4 U_{\zeta \zeta}^{0} + U_{\xi \xi}^{0} U_{\zeta \zeta}^{0} + U_{\zeta \zeta}^{0} ~U_{\eta \eta}^{0} -\alpha _{1}^{2}~ U_{\xi \xi}^{0} -\alpha _{1}^{2}~ U_{\eta\eta}^{0} -\alpha _{1}^{2} ~ U_{\zeta \zeta}^{0})}\\\\ +\alpha _{1} \lambda ((U_{\xi \eta}^{0})^{2} +(U_{\xi \zeta}^{0})^{2} -U_{\xi \xi}^{0} ~U_{\eta \eta}^{0} +(U_{\eta \zeta}^{0})^{2} -U_{\xi \xi}^{0}~ U_{\zeta \zeta}^{0}- U_{\eta \eta}^{0}~ U_{\zeta \zeta}^{0})\\\\ +(U_{\xi \zeta}^{0})^{2} ~ U_{\eta \eta}^{0} -2~U_{\xi \eta}^{0}~ U_{\xi \zeta}^{0}~U_{\eta \zeta}^{0} + U_{\xi \xi}^{0}~ (U_{\eta \zeta}^{0})^{2} +(U_{\xi \eta}^{0})^{2} ~U_{\zeta \zeta}^{0} -U_{\xi \xi}^{0}~ U_{\eta \eta}^{0}~ U _{\zeta \zeta}^{0} =0, \end{array} \end{equation} $$(12)

From the solution of the equation (12) we found that λ has complex values in all the cases and at least one of them has positive real value. Hence all the equilibrium points are unstable.

In this paper, we have investigated the effect of Albedo and oblateness of the primary in the circular restricted four-body problem with variable masses. We have determined the equations of motion, when the masses of the primaries as well as the infinitesimal body vary, which are different from the classical case by the variation parameters α₁, k, oblateness factor σ and the radiations effect p₁, p₂, p₃. Further the expression for the variation of Jacobi integral constant have been evaluated which is also depending on the variation parameter α₁. We have plotted equilibrium points, Poincaré surface of sections and basins of attraction in five cases (i. Case when third primary is placed at one of the lagrangian points of the classical circular restricted three-body problem, ii. Variation of mass, iii. Solar radiation pressure effect, iv. Albedo effect, v. Oblateness effect) by using Mathematica software. The equilibrium points during in-plane motion (Fig. 2), we found six equilibrium points in which three are collinear and three are non-collinear equilibrium points in the first case, in the variable mass case, we got one is collinear and rest five are non-collinear, while in the cases solar radiation pressure, albedo and oblateness, we found eight equilibrium points. During out-of-plane motions (i.e. ξ ≠ 0, η = 0, ζ ≠ 0)(Fig. 3) and (i.e. ξ = 0, η ≠ 0, ζ ≠ 0)(Fig. 4), we found three equilibrium points in all cases. The poincaré surface of sections have been determined in two phase spaces (ξ−ξ′ Fig. 5 (i) and η−η′ Fig. 5 (ii)). In ξ−ξ′-plane, the poincaré surface of sections are shrinking and η−η′-plane, the poincaré surface of sections are first expanding and then shrinking . The Newton-Raphson basins of attraction have been studied in all five cases (Fig. 6, Fig. 7, Fig. 8, Fig. 9, Fig. 10). We used color code for the classification of the equilibrium points on the (ξ−η)−plane. These points will be clearly visible in the zoomed part of all the figures near the lagrangian configuration (ie. Fig. 6 b, Fig. 7 b, Fig. 8 b, Fig. 9 b, Fig. 10 b). Finally, we have examined the stability of the equilibrium points in the circular restricted four variable bodies problem with oblateness and albedo effect and found that all the equilibrium points are unstable.

Fig. 7

Fig. 8

Fig. 9

Fig. 10