universe Article Bifurcation Analysis and Periodic Solutions of the HD 191408 System with Triaxial and Radiative Perturbations Fabao Gao 1,2, * and Ruifang Wang 1 1 School of Mathematical Science, Yangzhou University, Yangzhou 225002, China; wangruifang16@sina.com 2 Departament de Matemàtiques, Edifici C, Facultat de Ciències, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain * Correspondence: gaofabao@gmail.com  Received: 31 December 2019; Accepted: 16 February 2020; Published: 22 February 2020  Abstract: The nonlinear orbital dynamics of a class of the perturbed restricted three-body problem is studied. The two primaries considered here refer to the binary system HD 191408. The third particle moves under the gravity of the binary system, whose triaxial rate and radiation factor are also considered. Based on the dynamic governing equation of the third particle in the binary HD 191408 system, the motion state manifold is given. By plotting bifurcation diagrams of the system, the effects of various perturbation factors on the dynamic behavior of the third particle are discussed in detail. In addition, the relationship between the geometric configuration and the Jacobian constant is discussed by analyzing the zero-velocity surface and zero-velocity curve of the system. Then, using the Poincaré–Lindsted method and numerical simulation, the second- and third-order periodic orbits of the third particle around the collinear libration point in two- and three-dimensional spaces are analytically and numerically presented. This paper complements the results by Singh et al. [Singh et al., AMC, 2018]. It contains not only higher-order analytical periodic solutions in the vicinity of the collinear equilibrium points but also conducts extensive numerical research on the bifurcation of the binary system. Keywords: binary HD 191408; perturbation; bifurcation and chaos; periodic orbit 1. Introduction As we know, approximately two-thirds of the stars are part of the multistellar system in our galaxy. Due to their diversity and the unpredictable characteristics of planets in our solar system, exploring the planets in the multistellar system is of great interest for many scientists. The planets in a multiplanetary system and celestial bodies constitute N-body problems. Among the N-body problems, the three-body problem has been a hot research field because of its obvious practical value, and the most widely used model is the classical restricted three-body problem (RTBP). There are two important solutions in RTBP: Equilibrium points and periodic solutions. Several analytical and numerical methods for searching periodic solutions of restricted three-body problems can be found in the review article Musielak and Quarles [1]. Regarding the existence of periodic solutions of a circular RTBP, Gao and Zhang [2] gave a rigorous proof and found that the periodic solutions were mainly affected by factors such as the initial values and the masses of the two primaries. When the restricted three-body problem is perturbed by some other factors (such as light pressure, oblate spheroidal primaries, radiation, and albedo, etc.), it will become the restricted three-body problem, which is closer to the natural system. Singh and his collaborators have done a lot of excellent research in this field. For example, Singh et al. [3,4] studied the existence of equilibrium points Universe 2020, 6, 35; doi:10.3390/universe6020035 www.mdpi.com/journal/universe Universe 2020, 6, 35 2 of 22 and the periodic orbits around the triangular equilibrium points of the perturbed RTBP, where the larger primary and smaller primary are considered triaxial and oblate spheroidal bodies, respectively. Moreover, Tsirogiannis et al. [5] and Singh et al. [6] considered a modification of the RTBP with radiation and oblateness and studied the periodic motions around the collinear equilibrium points. Furthermore, semi-analytical solutions around these points for both planar and spatial cases were obtained using the Lindstedt–Poincaré method by Singh et al. [6]. Kalantonis et al. [7] considered the asymptotic motion to collinear equilibrium points of the RTBP with the oblateness and they computed an asymptotic orbit by using a fourth-order local analysis, numerical integration, and standard differential corrections. For the case of a planar circular RTBP, where the first primary is an oblate spheroid, and the oblateness coefficients affect the character of the orbits, Zotos [8] reported three types of orbits: Bounded, escaping, and collisional orbits. When the primaries are triaxial rigid bodies, Elshaboury et al. [9] investigated the basic dynamical characteristics of the RTBP and obtained the equilibrium points as well as some simple symmetric periodic orbits. For the RTBP with radiation and triaxiality, Jain et al. [10] analyzed the effects of these perturbations on periodic orbits under different energy constants. Furthermore, Idridi and Ullah [11] discussed the effects of the radiation and albedo on the existence of the noncollinear libration points in the elliptic RTBP, in which the oblateness of the second primary was considered. In the past few decades, researchers have observed multiplanetary systems using space telescopes, observational data, and statistics tools (see Chen [12] for details). The binary system consisting of two stars moving around their common barycenter is particularly worth our attention. The periodic orbits are important keys to understand the motion of the third particle in the binary system. When two asteroids were approximated as triaxial ellipsoids, Hou et al. [13] studied the forced periodic orbits around the triangular libration points in a binary asteroid system influenced by the solar radiation pressure. Recently, a numerical method was proposed to search for three-dimensional periodic orbits by Shi et al. [14]. After applying this method to binary asteroid 1999 KW4, they found five kinds of periodic orbits of the binary system. By quantifying the detected set of planet masses and orbits, Howard and Fulton [15] efficiently made planet discovery and characterization. Berardo et al. [16] presented new observations of HIP 41378, and its possible orbital periods were obtained through observation. Singh et al. [17] studied the collinear equilibrium points and periodic motion around them in the RTBP for the binary HD 191408 system, where the two primaries are triaxial rigid bodies and emit radiation. Moreover, the effects of different parameters on the collinear equilibrium points were discussed. Das et al. [18] investigated the field of the radiating binary stellar system in the circular RTBP. Singh et al. [19] found three-dimensional periodic orbits around the collinear equilibrium points of the RTBP with oblate and radiating primaries. Singh and Umar [20] found that the positions of the third particle depended on the oblateness, radiation coefficients of the primaries, and the eccentricity of their orbits in the elliptic RTBP. They provided the numerical application of this problem in the stellar-oblate binary system. HD 191408 is a high-velocity star that belongs to the southern hemisphere main-sequence stars with debris disks, and García and Gómez [21] studied their optical aperture polarimetry. The study of high-resolution, high-signal-to-noise spectra of field stars of different metallicities becomes an effective technique to tackle various problems related to the chemical evolution of the galaxy, and Abia et al. [22] provided the atmospheric parameters, elemental abundance ratio, and signal-to-noise ratios for some stars, including the parameters of HD 191408. Karaali et al. [23] investigated the metallicity calibration of several dwarfs and metal-poor stars at different distances from the galactic plane, which contribute to the implications for the galactic formation and evolution. In addition, Perrin [24] analyzed the chemical composition for twelve K dwarfs, whose masses were estimated. However, few researchers have studied the periodic solutions of this system. Bifurcation refers to the motion of a system with suddenly changing parameters, e.g., the equilibrium state or the number and stability of the periodic motion, when the parameters change and become a certain value. Yumagulov et al. [25] studied the bifurcation in the planar elliptical RTBP, where Universe 2020, 6, 35 3 of 22 the parameters of the system around the triangular equilibrium points were discussed. Perdomo [26] observed the reduced periodic solutions of the spatial isosceles three-body problem, which contained a bifurcation point and provided an explanation for the existence of this point. Maciejewski and Rybicki [27] studied global bifurcation from the equilibrium points of the nonstationary periodic solutions in the RTBP. There are also some other excellent related literatures [28–31], mainly including the reports of a binary system with short-orbital-period, the spatial periodic orbits in various resonances, and the periodic motions near a high mass ratio binary star system. In this paper, we inherit the model of the perturbed restricted three-body problem of Singh et al. [17] and Jain et al. [10] and continue to study the binary HD 191408 system. In Section 2, the dynamical equations that involve the parameters of the third particle in the binary system are obtained. In Section 3, the bifurcation diagrams of the system’s state variables in terms of different parameters are illustrated and explained. In Section 4, the equilibrium points of the system are introduced by discussing the geometric configurations. In Section 5, the 2- and 3-dimensional periodic orbits around the collinear equilibrium points are computed using the Lindstedt–Poincaré method. Section 6 summarizes the study. 2. Dynamical Equations For the sake of convenience, we first introduce the normal rotating dimensionless coordinate system oxyz, of which the origin is located in the centroid of two primary bodies with masses m1 = 1 − µ and m2 = µ, where µ = m2 /(m1 + m2 ) ∈ (0, 1). Assuming that the coordinate of the third particle is (x, y, z), the two primaries are located at (µ, 0, 0) and (µ − 1, 0, 0), respectively. According to Singh et al. [17] and Jain et al. [10], the equations of motion of the third particle are .. . x − 2αn y = Ωx , .. . y + 2αnx = Ω y , (1) .. z = Ωz , where α represents the modification of perturbation on the Coriolis force and the potential function Ω(x, y, z) admits the following form n2 β 2 (1−µ)q1 µq (1−µ)(2s1 −s2 )q1 µ(2s3 −s4 )q2 Ω(x, y, z) = 2 (x + y ) + 2 r1 + r22 + 2r1 3 + 2r2 3 3(1−µ)(s1 −s2 ) y2 q1 3µ(s3 −s4 ) y2 q2 3(1−µ)s1 z2 q1 3µs3 z2 q2 (2) − 2r1 5 − 2r2 5 − 2r1 5 − 2r 5 , 2 β denotes the modification of perturbation on the centrifugal force and q1 and q2 are the radiation factors q q of the primaries. In addition, r1 = (x − µ)2 + y2 + z2 and r2 = (x + 1 − µ)2 + y2 + z2 denote the distances between the third particle and the first and second primaries, respectively. On account of the triaxiality of the primaries, the mean perturbed motion n is defined by r 3 n= 1 + [(2s1 − s2 ) + (2s3 − s4 )], (3) 2 (a1 2 − a3 2 ) (a2 2 − a3 2 ) (a1 02 − a3 02 ) (a2 02 − a3 02 ) s1 = , s 2 = , s3 = , s4 = , 5R2 5R2 5R2 5R2 where ai (i = 1, 2, 3) is the semi-axis of the first primary, ai 0 (i = 1, 2, 3) is the semi-axis of the second primary, and R is the dimensional distance between the two primaries. The physical parameters of the binary HD 191408 system are shown in Table 1. Universe 2020,6,6,3535 Universe2020, 4 4ofof2422 Table 1. Physical parameters of the binary HD 191408 system [17]. Table 1. Physical parameters of the binary HD 191408 system [17]. Parameters Parameters µ s µ1 s s1 2 s s2 3 s3 4 s s4 1q q1 2 q q α α β β 2 Values 0.1881 0.008 0.004 0.002 0.001 0.407761 0.991476 1.0003 1.0002 Values 0.1881 0.008 0.004 0.002 0.001 0.407761 0.991476 1.0003 1.0002 The first Jacobi-type function of the system (1) is of (the The first Jacobi-type function2Ω , z ) − (1) x, ysystem x 2 is ( + y 2 + z 2 = C, ) (4) where 𝑉𝑉 2 = 𝑥𝑥̇ 2 + 𝑦𝑦̇ 2 + 𝑧𝑧̇ 2 is the motion velocity y, z) − of x the + ythird = C , and C is the Jacobi constant. + z particle  . 2 . 2 . 2 2Ω(x, (4) 3. Analysis of. the 2 Bifurcation .2 .2 and Chaos where V 2 = x + y + z is the motion velocity of the third particle and C is the Jacobi constant. In this section, the bifurcation diagrams of the state variables in terms of different parameters in 3. Analysis the binary HD of 191408 the Bifurcation system are andillustrated, Chaos and the effects of each parameter on the system are analyzed. In this section, the bifurcation diagrams of the state variables in terms of different parameters in the For binary HD the given 191408 systemvalues of Tableand1,thewe are illustrated, effectsselected of each the iterative parameter on theinitial systemvalue 𝑥𝑥0 = are analyzed. (−0.08, 0.001, 0.2, 0, −1.6, 0) and limit parameter 𝜇𝜇 in the interval (0, 1). By applying For the given values of Table 1, we selected the iterative initial value x0 = (−0.08, 0.001, 0.2, 0, −1.6, 0) the ode45 numerical integration and limit parameter algorithm µ in [32](0,, we the interval obtain 1). By applying the results the ode45 in numerical Figure 1; integration namely, thealgorithm bifurcation [32], diagrams of three state variables, frames of 𝜇𝜇𝜇𝜇𝜇𝜇, 𝜇𝜇𝜇𝜇𝜇𝜇, and 𝜇𝜇𝜇𝜇𝜇𝜇, in terms of we obtain the results in Figure 1; namely, the bifurcation diagrams of three state variables, frames of mass parameter 𝜇𝜇 in the binary HDand µxy, µyz, 191408 µxz,system. in termsFigure of mass 1aparameter shows theµeffectin theofbinary 𝜇𝜇 in the HD 𝜇𝜇𝜇𝜇𝜇𝜇 191408frame; for aFigure system. given1avalue shows of the 𝜇𝜇, the bifurcation effect of µ in thediagram µxy frame; in for thea 𝜇𝜇𝜇𝜇𝜇𝜇 givenframe value in terms of µ, of 𝜇𝜇 is obtained. the bifurcation diagram Figure 1a frame in the µxy also shows in terms theof bifurcation µ is obtained. diagrams Figure 1a of also variables shows the𝑥𝑥 and 𝑦𝑦 in terms bifurcation diagrams of 𝜇𝜇of, variables i.e., the projection from of x and y in terms the µ, three- i.e., the dimensional projection from diagram to the two-dimensional the three-dimensional diagram plane. to the Thetwo-dimensional projections haveplane. similarThestructure, which projections haveis inevitable in qualitative analysis by comparing the projections. Similarly, the similar structure, which is inevitable in qualitative analysis by comparing the projections. Similarly, bifurcation diagrams in the the𝜇𝜇𝜇𝜇𝜇𝜇 and 𝜇𝜇𝜇𝜇𝜇𝜇diagrams bifurcation frames in interms the µyzof and 𝜇𝜇 are µxzshown framesininFigureterms1b,c.of µ are shown in Figure 1b,c. (a) Figure 1. Cont. Universe 2020, 6, 35 5 of 24 Universe 2020, 6, 35 5 of 22 (b) (c) Figure 1. (a–c). Bifurcation diagrams of the µxy, µyz, and µxz frames in terms of mass parameter Figure 1. (a,b,c). Bifurcation diagrams of the 𝜇𝜇𝑥𝑥𝑥𝑥, 𝜇𝜇𝑦𝑦𝑦𝑦, and 𝜇𝜇𝑥𝑥𝑥𝑥 frames in terms of mass parameter µ, respectively. 𝜇𝜇, respectively. Universe 2020, 6, 35 6 of 22 Universe 2020, 6, 35 6 of 24 For Foreach eachdifferent differentinitial initialvalue, value,system system(1) (1)will willadmit admitdifferent differentbifurcation bifurcationdiagrams, diagrams,but butfor foraa given given initial value, we can analyze the influence of the change of perturbation parameters onthe initial value, we can analyze the influence of the change of perturbation parameters on the dynamic behavior of the third particle by displaying a set of bifurcation dynamic behavior of the third particle by displaying a set of bifurcation diagrams. For example, we diagrams. For example, we can canfind findclearly clearlyfrom fromFigure Figure11that thatchaotic chaoticmotion motionexists existsininthe thethree threedirections directionsof ofspace spaceforforthethethird third particle particlewhen when thethe mass parameter µ𝜇𝜇 ∈∈((0, mass parameter 0.5]. 0, 0.5 Especially when ]. Especially whenµ𝜇𝜇==0.5 0.5ororso, so,that thatis is to to say, say, when when the the masses of two primaries are approximately equal, the dynamic masses of two primaries are approximately equal, the dynamic behavior of the third particle changes behavior of the third particle changes greatly. greatly. However, However, its dynamicits dynamic behaviorbehavior decreases decreases significantlysignificantly when µ ∈when (0.5, 𝜇𝜇1)∈. (0.5, 1). Similarly, weselectselectx 𝑥𝑥and and restrict the triaxial coefficient 𝑠𝑠𝑖𝑖 (𝑖𝑖3,=41, 2, 3, to the[−0.5, the4)interval interval Similarly, we 0 0 restrict the triaxial coefficient si (i = 1, 2, ) to 0.5]. [−0.5, 0.5]. The bifurcation diagrams of the state variables with respect The bifurcation diagrams of the state variables with respect to the triaxial coefficient si (i = 1, 2, 3, 4) to the triaxial coefficient 𝑖𝑖 (𝑖𝑖 shown 𝑠𝑠are are shown = 1, 2, 3,in4)Figures 2 andin3.Figures Figure 22aand shows3. Figure that the2a triaxial shows that the triaxial coefficient s1 of coefficient 𝑠𝑠1 of the the first primary has first primary has little effect on the motion amplitude of the third particle little effect on the motion amplitude of the third particle in the y direction, although it will cause the in the 𝑦𝑦 direction, although itoverall will cause the overallofdisplacement displacement the third particle of theinthird this particle direction. in It this hasdirection. It has a great a great influence on theinfluence movement on the of themovement of the third third particle in theparticle in theespecially x direction, 𝑥𝑥 direction, whenespecially s1 takeswhen 𝑠𝑠1 takes the the threshold threshold value around value ±0.32, around the movement of the third particle in the x direction will be the opposite. Figure 2b shows that 2b ±0.32, the movement of the third particle in the 𝑥𝑥 direction will be the opposite. Figure the shows triaxialthat the triaxial coefficient s2 ofcoefficient 𝑠𝑠2 of the the first primary willfirst causeprimary will cause the motion the motion amplitude amplitude of the third particleof to theexpand third particle gradually to expand in the ygradually direction,inand thewhen𝑦𝑦 direction, s2 takesand the when threshold 𝑠𝑠2 takes valuethe of threshold −0.07 in the value of −0.07 the x direction, in the 𝑥𝑥 direction, the movement of the third particle in this direction movement of the third particle in this direction will also be in the opposite direction. Figure 2c shows will also be in the opposite direction. Figure 2c shows that the combination that of the the combination triaxial coefficientsofs1the andtriaxial s2 actscoefficients on the third and 𝑠𝑠2 in 𝑠𝑠1 particle acts theonx the third direction. particle The result in the shows𝑥𝑥 direction. The result that the effect on theshowsdynamic thatbehavior the effect of on thethe thirddynamic particlebehavior is greaterofwhenthe third their particle is greater when their values values pass the threshold s = −0.32 and s = −0.07. pass the threshold 𝑠𝑠1 = −0.32 and 𝑠𝑠2 = −0.07. 1 2 (a) Figure 2. Cont. Universe 2020, 6, 35 7 of 24 Universe 2020, 6, 35 7 of 22 (b) (c) (a–c).Bifurcation Figure2.2.(a,b,c). Figure Bifurcation diagrams diagrams in in𝑠𝑠1 s𝑥𝑥𝑥𝑥, 1 xy,𝑠𝑠 s𝑥𝑥𝑥𝑥, 2 2 xy,and and𝑠𝑠 s𝑠𝑠1 s2𝑥𝑥x frames 1 2 frames with withrespect respecttotothe thetriaxial triaxial coefficient s ( i = 1, 2 ) , respectively. coefficient 𝑠𝑠𝑖𝑖 (𝑖𝑖 = 1, 2 ), respectively. i direction. In the 𝑥𝑥 direction, when the value of 𝑠𝑠3 passes the threshold value of −0.2, the movement of the third particle will change greatly. Figure 3b shows that the effect of 𝑠𝑠4 on the movement of the third particle in 𝑥𝑥 and 𝑦𝑦 directions is extremely complex, and obvious chaos has appeared. Figure 3c shows that under the joint action of 𝑠𝑠3 and 𝑠𝑠4 , with the increase of 𝑠𝑠3 and the decrease of 𝑠𝑠4 , the dynamic behavior of the third particle in the 𝑥𝑥 direction tends to be stable. Universe 2020, 6, 35 8 of 22 Universe 2020, 6, 35 9 of 24 (a) (b) Figure 3. Cont. Universe 2020, 6, 35 9 of 22 (b) (c) Figure 3. (a–c). Bifurcation diagrams in s3 xy, s3 xy, and s3 s4 x frames with respect to the triaxial Figure 3. (a,b,c). Bifurcation diagrams in 𝑠𝑠3 𝑥𝑥𝑥𝑥, 𝑠𝑠3 𝑥𝑥𝑥𝑥, and 𝑠𝑠3 𝑠𝑠4 𝑥𝑥 frames with respect to the triaxial coefficient si (i = 3, 4), respectively. coefficient 𝑠𝑠𝑖𝑖 (𝑖𝑖 = 3, 4), respectively. Figure 3 reflects the effect of the three semi-axes perturbations of the second primary on the dynamics of the third particle. Figure 3a shows that s3 has little effect on the movement amplitude of the third particle in the y direction, except that the third particle will have an overall shift in this direction. In the x direction, when the value of s3 passes the threshold value of −0.2, the movement of the third particle will change greatly. Figure 3b shows that the effect of s4 on the movement of the third particle in x and y directions is extremely complex, and obvious chaos has appeared. Figure 3c shows that under the joint action of s3 and s4 , with the increase of s3 and the decrease of s4 , the dynamic behavior of the third particle in the x direction tends to be stable. Next we select x0 and restrict the radiation coefficient qi (i = 1, 2) to the interval [0, 1]. Figure 4 reflects the influence of the radiation factors qi (i = 1, 2) of the primaries on the dynamic behavior of the third particle. Figure 4a shows that the effect on the third particle in the x and y directions is the same basically, both of which are gradually reduced from large-scale movements as the radiation factor q1 increases. It can be found from Figure 4b that the change of radiation factor q2 has little effect on the dynamic behavior of the third particle, but under the combined effect of radiation factors q1 and q2 (see Figure 4c), the third particle exhibits periodic motion in the x direction when q1 ∈ [0.7, 0.8]. the same basically, both of which are gradually reduced from large-scale movements as the radiation factor 𝑞𝑞1 increases. It can be found from Figure 4b that the change of radiation factor 𝑞𝑞2 has little effect on the dynamic behavior of the third particle, but under the combined effect of radiation factors 𝑞𝑞1 and 𝑞𝑞2 (see Figure 4c), the third particle exhibits periodic motion in the 𝑥𝑥 direction when 𝑞𝑞1 ∈ Universe 2020, 6, 35 10 of 22 [0.7, 0.8]. Universe 2020, 6, 35 11 of 24 (a) (b) Figure 4. Cont. Universe 2020, 6, 35 (b) 11 of 22 (c) 4. (a–c). FigureFigure Bifurcation diagrams in q1 xy,𝑞𝑞 q𝑥𝑥𝑥𝑥, 4. (a,b,c). Bifurcation diagrams in 𝑞𝑞1 𝑥𝑥𝑥𝑥, 2 xy, and q q x frames with respect to the radiation and 𝑞𝑞1 𝑞𝑞21𝑥𝑥 2frames with respect to the radiation 2 coefficient q ( coefficient i i = 1, 2 ) , respectively. 𝑞𝑞𝑖𝑖 (𝑖𝑖 = 1, 2), respectively. Universe 2020, 6, 35 12 of 24 For the values of the Coriolis and centrifugal forces, which are restricted in the interval [1, 1.1]. For the The bifurcation values of diagram inthe Coriolis Figure 5aand centrifugal shows that the forces, which are restricted modification in the interval in the Coriolis force α[1,has 1.1]. a similar The bifurcation diagram in Figure 5a shows that the modification in the Coriolis force effect on the third particle in the x and y directions, both of which tend to a certain range gradually 𝛼𝛼 has a similar effect on the third particle in the 𝑥𝑥 and 𝑦𝑦 directions, both of which tend to a certain range gradually with the increase of α from large-scale motion. From Figure 5b, we find that the movement of the with the increase of 𝛼𝛼 from large-scale motion. From Figure 5b, we find that the movement of the third particle in theinx the third particle and𝑥𝑥 yanddirections appears 𝑦𝑦 directions divergent appears divergentasas the modification the modification in the in the force force β centrifugal centrifugal 𝛽𝛽 increases; increases; thatmovement that is, the is, the movement amplitude amplitude in in bothdirections both directions increases. increases.However, However, the third theparticle third particle shows periodic shows periodic motion motion in the in x the 𝑥𝑥 direction direction with with thetheincrease increase of and 𝛽𝛽. α and of 𝛼𝛼 β. (a) Figure 5. Cont. Universe 2020, 6, 35 12 of 22 (a) Universe 2020, 6, 35 13 of 24 (b) (c) Figure 5. (a–c). Bifurcation diagrams in αxy, βxy, and αβx frames with respect to the Coriolis α and the Figure 5. (a,b,c). Bifurcation diagrams in 𝛼𝛼𝛼𝛼𝛼𝛼, 𝛽𝛽𝛽𝛽𝛽𝛽, and 𝛼𝛼𝛼𝛼𝛼𝛼 frames with respect to the Coriolis 𝛼𝛼 centrifugal force β, respectively. and the centrifugal force 𝛽𝛽, respectively. 4. Equilibrium Points For the Jacobi-type function (4) of the system, when the motion velocity of the third particle is zero, the three-dimensional space of the system with the change of Jacobi constant C is shown in Figures 6–10. When C decreases from C = 3, five equilibrium points A, B, D, E, and F can be obtained. A smaller C value corresponds to the larger permissible regions of motion of the third particle. The specific analysis is as follows: When C = 3, the zero-velocity surface of the third particle is shown in Figure 6. The third particle For the Jacobi-type function (4) of the system, when the motion velocity of the third particle is zero, the three-dimensional space of the system with the change of Jacobi constant C is shown in Figures 6–10. When C decreases from C = 3, five equilibrium points A, B, D, E, and F can be obtained. A smaller C value corresponds to the larger permissible regions of motion of the third particle. Universe 2020, 6, 35 13 of 22 The specific analysis is as follows: When C = 3, the zero-velocity surface of the third particle is shown in Figure 6. The third particle 4. Equilibrium can only skim overPointstwo primaries under the action of gravity but cannot pass through the forbidden area around them. For the When Cfunction Jacobi-type = 2.914, (4) theofforbidden the system,area of one when of the primaries the motion velocity ofand the the outer third forbidden particle is one will intersect at point A (see Figure 7). The third particle can fly to the outer space through the zero, the three-dimensional space of the system with the change of Jacobi constant C is shown in channel A. In Figures fact, 6–10. this Cpoint When A is the decreases from C =equilibrium first point.points A, B, D, E, and F can be obtained. 3, five equilibrium A smaller C value corresponds to the larger permissible regions of motion of the third particle. The specific analysis is as follows: Universe 2020, 6, 35 14 of 24 Universe 2020, 6, 35 14 of 24 Figure 6. Zero-velocity surface of the third particle C = 3. Figure 6. Zero-velocity surface of the third particle C = 3. Figure 7. Zero-velocity surface of the third particle C = 2.914. When C is 2.1943, as shown in Figure 8, the prohibited region of the third particle decreases. The “Channel B” appears, through which the third particle can fly from the permissible regions of one primary to another. Meanwhile, the second equilibrium point can be obtained. When C drops to 1.9777 (see Figure 9), “Channel D” appears, where the third particle can fly into another outer space. Thus, the third equilibrium point appears. Figure 7.7.Zero-velocity Figure Zero-velocity surface ofthe surface of thethird thirdparticle particleC =C2.914. = 2.914. When C is 2.1943, as shown in Figure 8, the prohibited region of the third particle decreases. The “Channel B” appears, through which the third particle can fly from the permissible regions of one primary to another. Meanwhile, the second equilibrium point can be obtained. When C drops to 1.9777 (see Figure 9), “Channel D” appears, where the third particle can fly into another outer space. Thus, the third equilibrium point appears. Figure 8. Zero-velocity surface of the third particle C = 2.1943. Figure 8. Zero-velocity surface of the third particle C = 2.1943. Figure 8. Zero-velocity surface of the third particle C = 2.1943. Universe 2020, 6, 35 14 of 22 Universe 2020, 6, 35 15 of 24 When C = 1 (see Figure 10), the last two equilibrium points E and F appear. With the decrease of C, the third particle can remove the primaries to fly into outer space. Figure 9.9.Zero-velocity Figure Zero-velocity surface ofthe surface of thethird thirdparticle particleC =C1.9777. = 1.9777. Figure 10. Zero-velocity surface of the third particle C = 1. Figure 10. Zero-velocity surface of the third particle C = 1. When C = 3, the zero-velocity surface of the third particle is shown in Figure 6. The third particle According can only skimtoover the two discussion primariesofunder zero velocity the actionsurface, when of gravity C takespass but cannot 2.914, 2.1943, through theand 1.9777, the forbidden area around them. When C = 2.914, the forbidden area of one of the primaries and related zero velocity curves are shown in Figure 11. Therefore, three collinear equilibrium the outer forbidden points are one will intersect at point A (see Figure 7). The third particle can fly to the outer space through the 𝐿𝐿1 (−0.34630003, 0), 𝐿𝐿2 (0.84584372, 0), and 𝐿𝐿3 (−1,23037223, 0), respectively. channel A. In fact, this point A is the first equilibrium point. When C is 2.1943, as shown in Figure 8, the prohibited region of the third particle decreases. The “Channel B” appears, through which the third particle can fly from the permissible regions of one primary to another. Meanwhile, the second equilibrium point can be obtained. When C drops to 1.9777 (see Figure 9), “Channel D” appears, where the third particle can fly into another outer space. Thus, the third equilibrium point appears. When C = 1 (see Figure 10), the last two equilibrium points E and F appear. With the decrease of C, the third particle can remove the primaries to fly into outer space. According to the discussion of zero velocity surface, when C takes 2.914, 2.1943, and 1.9777, the related zero velocity curves are shown in Figure 11. Therefore, three collinear equilibrium points are L1 (−0.34630003, 0), L2 (0.84584372, 0), and L3 (−1, 23037223, 0), respectively. Figure 11. Zero-velocity curves of the third particle for C = 2.914, 2.1943, and 1.9777, respectively. 5. Periodic Orbits Around the Collinear Equilibrium Points 5.1. Expansion of The 2D Equations of Motion In the plane, the motion equations of the third particle are (see [17]) x1 − 2α nx2 =  Ω x1 , (5) x2 + 2α nx1 =  Ω x2 , Figure 10. Zero-velocity surface of the third particle C = 1. According to the discussion of zero velocity surface, when C takes 2.914, 2.1943, and 1.9777, the related zero velocity curves are shown in Figure 11. Therefore, three collinear equilibrium points Universe 2020, 6, 35 15 of 22 are 𝐿𝐿1 (−0.34630003, 0), 𝐿𝐿2 (0.84584372, 0), and 𝐿𝐿3 (−1,23037223, 0), respectively. Figure 11. Zero-velocity curves of the third particle for C = 2.914, 2.1943, and 1.9777, respectively. Figure 11. Zero-velocity curves of the third particle for C = 2.914, 2.1943, and 1.9777, respectively. 5. Periodic Orbits Around the Collinear Equilibrium Points 5. Periodic Orbits Around the Collinear Equilibrium Points 5.1. Expansion of The 2D Equations of Motion 5.1. In Expansion the plane, of theThe 2D Equations motion of Motion equations of the third particle are (see [17]) In the plane, the motion equations.. of the third . particle are (see [17]) x1 − 2αnx2 = Ωx1 , − 2xα (5) x12αn  nxΩ 2x2= , Ωx , .. . x2 + 1 = 1 (5) x2 +in2α  where x1 = x, x2 = y and the potential function nx1 = Equation (5)x2is, Ω n2 β (1−µ)q µq (1−µ)(2s −s )q µ(2s3 −s4 )q2 Ω= 2 (x1 2 + x2 2 ) + ra 1 + rb2 + 1 2 1 2ra 3 + 2rb 3 = where x1 x= , x2 y and the potential − function 3(1−µ)( s1 −s2 )x2 qin 2 1 Equation − 3µ(s3 −s4 )x2(5) 2 q is 2 . (6) 2r 5 a 2r 5 b n2 β 2 2 (1 − µ )q1 µ q2 (1 − µ )(2 s1 − s2 )q1 µ (2 s3 − s4 )q2 = Deriving the ( x1 with Ω derivation + x2respect ) + to x1 and+ x2 on+the right side + (6), we of Equation obtain 2 ra rb 2ra 3 2rb 3 (1−µ)q1 (x1 −µ) µq2 (x1 +1−µ) 3(1−µ)(2s1 −s2 )q1 (x1 −µ) (6) Ωx1 3(1 = n−2 βx µ )( 1− s − s2 )rax32 2 q1 −3 ( sr3b 3− s4 ) −µ x2 2 q2 2r 5 −(2s3 −s4 )q2 (x1 +1−µ 3µ 5 ) 1 15(1−µ−)(s1 −s2 )x2 2 q1 (5x1 −µ) .15µ(s3 −sa 4 )x2 2 q2 (x1 +1−µ) 2rb 5 2ra 2ra 7 2rb − + + 2rb 7 , ( 1−µ ) q x µq x 3 ( 1−µ )( 4s −3s ) q x 3µ ( 4s −3s ) q x Ωx2 = n2 βx2 − ra 3 1 2 − r 23 2 − 1 2ra 5 2 1 2 − 3 2rb 5 4 2 2 b 15(1−µ)(s1 −s2 )q1 x2 3 15µ(s3 −s4 )q2 x2 3 + 2ra 7 + 2rb 7 , q q where ra = (x1 − µ) + x2 rb = (x1 + 1 − µ)2 + x2 2 are the distances between the third particle 2 2, and the primaries. By substituting the transformation of x1 = xLi + ξ and x2 = η into Equation (5), we obtain the following equations of motion of the third particle .. . ξ − 2αnη = Ωξ , .. . (7) η + 2αnξ = Ωη . Expand the right side of Equation (7) to the second order and get .. . ξ − 2αnη = K1 ξ + K2 ξ2 + K3 η2 , .. . (8) η + 2αnξ = L1 η + L2 ξη, and to the third order, as follows .. . ξ − 2αnη = K1 ξ + K2 ξ2 + K3 η2 + K4 ξ3 + K5 ξη2 , .. . (9) η + 2αnξ = L1 η + L2 ξη + L3 ξ2 η + L4 η3 , where the coefficients as shown in Appendix A. Universe 2020, 6, 35 16 of 22 5.2. Periodic Orbits in the Plane Now, we examine the periodic solutions of system (5). Suppose that the periodic solutions of system (8) has the following form in the powers of orbital parameter e (|e|  1) ξ = ξ1 e + ξ2 e2 , (10) η = η1 e + η2 e2 , Universe 2020, 6, 35 17 of 24 By substituting Equation (10) into (8) and the series expansion of periodic solution with respect to ewhere (ξ0 , 0 ) denotes the position is as follows of the collinear equilibrium points in the previous section. ξ = ξ0 + cos(ωt)e + [$21 + $23 cos(2ωt)]e2 , (11) Similarly, we suppose that the η =periodic v12 sin(ωt solutions )e + v24of (2ωt)e2(9) sinsystem , in powers of 𝑒𝑒 are where (ξ0 , 0) denotes the position of the collinear ξ =ξ1e +equilibrium ξ 2 e 2 + ξ3e3 points , in the previous section. Similarly, we suppose that the periodic solutions of system (9) in powers of e are (12) η =η1e + η2 e 2 + η3e3 , ξ = ξ1 e + ξ2 e2 + ξ3 e3 , (12) η = η1 e + η2 e2 + η3 e3 , Thus, the periodic solution of system (5) into series expansions of 𝑒𝑒 up to the third-order terms is ξ 0 + cos(ωt )e + [ϖ 21 + ϖ 23 cos(2ωt ) ] e + ϖ 34 cos(3ωt )e , Thus, the periodic solution of system (5) into series expansions of e up to the third-order terms is ξ= 2 3 (13) ηξ = vξ120 + cosω(tωt )e)e++v24[$sin(2 21 + $ ω23t )2cos e +ω $t34)ecos 2 3, = sin( e 2 +(2ωt v36 )]sin(3 3 (3ωt)e , (13) η = v12 sin(ωt)e + v24 sin(2ωt)e + v36 sin(3ωt)e , 3 where the coefficients are shown in Appendix B. where the coefficients are shown in Appendix B. For the the collinear collinear equilibrium equilibrium point point L𝐿𝐿1((−0.34630003, 0) around which the periodic orbits are For 1 −0.34630003, 0) around which the periodic orbits are obtained,as obtained, asshown shown ininFigure Figure12, 12,the thered redcurve curveisisup upto tothe thesecond-order second-orderterms, terms,and andthe theblue bluecurve curveisis about the third-order ones. By comparing these two orbits, we find that these two orbits are about the third-order ones. By comparing these two orbits, we find that these two orbits are close closeto to each other. each other. Figure 12. Periodic orbits around the collinear equilibrium point L1 up to second-order terms (red curve) and third-order terms (blue curve). Figure 12. Periodic orbits around the collinear equilibrium point 𝐿𝐿1 up to second-order terms (red curve) and third-order terms (blue curve). 5.3. Expansions of the Three-Dimensional Equations of Motion Substituting the transformations of x xLi + ξ , y = η and z = ζ into Equation (1), we = obtain the motion equations of the third particle in the (ξ ,η , ζ ) coordinate system Universe 2020, 6, 35 17 of 22 5.3. Expansions of the Three-Dimensional Equations of Motion Substituting the transformations of x = xLi + ξ, y = η and z = ζ into Equation (1), we obtain the motion equations of the third particle in the (ξ, η, ζ) coordinate system .. . ξ − 2αnη = Ωξ , .. . η + 2αnξ = Ωη , (14) .. ζ = Ωζ . Using Taylor expansion, the RHSs of Equation (14) are expanded to the second-order terms (see [17]); thus, the motion equations of the third particle become .. . ξ − 2αnη = K1 ξ + K2 ξ2 + K3 η2 + K6 ζ2 , .. . η + 2αnξ = L1 η + L2 ξη, (15) ζ = M1 ζ + M2 ξζ, Here, we use Taylor expansion, and the RHSs of Equation (14) are expanded to the third-order terms. Then, the motion equations of the third particle become .. . ξ − 2αnη = K1 ξ + K2 ξ2 + K3 η2 + K6 ζ2 + K7 ξ3 + K8 ξη2 + K9 ξζ2 , .. . η + 2αnξ = L1 η + L2 ξη + L5 η3 + L6 ξ2 η + L7 ηζ2 , (16) ζ = M1 ζ + M2 ξζ + M3 ζ3 + M4 ξ2 ζ + M5 η2 ζ, where the coefficients are shown in Appendix A. 5.4. Periodic Orbits in the Spatial Space Using the method of successive approximations up to the second-order terms, the periodic solutions of system (14) are in the following form (see [17]) ξ = ξ2 e2 , η = η 2 e2 , (17) ζ = ζ1 e, So, the periodic solutions in the form of e up to the second-order terms can be written as ξ = ξ0 + [$02 + $03 cos(2ωt)]e2 , η = v05 sin(2ωt)e2 , (18) ζ = sin(ωt)e. Similarly, using the same method up to the third-order terms, the periodic solutions of the above system are ξ = ξ1 e + ξ2 e2 + ξ3 e3 , η = η1 e + η2 e2 + η3 e3 , (19) ζ = ζ1 e + ζ2 e + ζ3 e , 2 3 The rationality of doing this is that the more theoretically expanded, the closer it is to the exact solution. However, the problem is that the more we expand, the greater the challenge we will face in obtaining the analytical solution, so we try our best to solve the approximate analytical solution by semi-analytical method. The periodic solutions in the form of e up to the third-order terms can be obtained ξ = ξ0 + cos(ωt)e + [$20 + $24 cos(2ωt)]e2 + $32 cos(3ωt)e3 , η = v12 sin(ωt)e + v21 sin(2ωt)e2 + v34 sin(3ωt)e3 , (20) ζ = sin(ωt)e + v32 sin(2ωt)e2 + v33 sin(3ωt)e3 . Universe 2020, 6, 35 19 of 24 ξ=ξ 0 + cos(ωt )e + [ϖ 20 + ϖ 24 cos(2ωt ) ] e 2 + ϖ 32 cos(3ωt )e3 , Universe 2020, 6, 35 η=v12 sin(ωt )e + v21 sin(2ωt )e 2 + v34 sin(3ωt )e3 , 18(20) of 22 ζ = sin(ωt )e + v32 sin(2ωt )e 2 + v33 sin(3ωt )e3 . where the coefficients can be found in Appendix B. whereThetheperiodic coefficients canaround orbits be foundtheincollinear Appendix B. equilibrium point L1 are plotted in the spatial space as The periodic orbits around the collinear equilibrium point 𝐿𝐿1 are plotted shown in Figure 13. Figure 13a is a periodic orbit up to the second-order in the terms, andspatial Figurespace 13b isasa shown periodic orbit up to the third-order terms. Because we selected different solutions in terms of e,isi.e., in Figure 13. Figure 13a is a periodic orbit up to the second-order terms, and Figure 13b a periodic orbit up to the third-order terms. Because we selected different solutions in terms effects of various of factors of deformation, orbital vibration, and motion, there are different formsof 𝑒𝑒, i.e., effects of various of factors of deformation, orbital vibration, and motion, there are different forms of of orbits. orbits. (a) (b) Figure Figure13. 13.(a,b) Periodic (a,b) orbits Periodic around orbits 𝐿𝐿1 in around L1the in spatial spacespace the spatial up to up second and third to second andorder thirdterms, order respectively. terms, respectively. Universe 2020, 6, 35 19 of 22 6. Conclusions In this paper, we focus on constructing the approximate analytical periodic solutions of the binary HD 191408 system by using Lindstedt-Poincaré method. The obtained second- and third-order periodic solutions in the plane and three-dimensional space generalized the corresponding ones in [17]. In addition, we perform extensive numerical research on the bifurcation of the system to discuss the effects of nine perturbations on the third particle’s dynamic behavior. These nine parameters include mass ratio µ, triaxial coefficients si (i = 1, 2, 3, 4), radiation factors q1 and q2 , the modification in the Coriolis force α, as well as the modification in the centrifugal force β. The results show that when the mass ratio parameter µ < 0.5, it has a greater effect on the third particle’s dynamic behavior than when µ > 0.5, and it has the strongest effect on the dynamic behavior of the third particle when µ = 0.5. Furthermore, the triaxial coefficients s1 and s2 have a more significant impact on the system than s3 and s4 . This mainly reflects that the dynamic behavior of the third particle changes greatly when the system is under the joint action of s1 and s2 , but it gradually stabilizes under the combined action of s3 and s4 . Furthermore, compared with the corresponding radiation factor q2 and the modification in the Coriolis force α, we also find that radiation factor q1 and the modification in the centrifugal force β have a greater impact on the dynamic behavior of the third particle, but with the increase of q1 and q2 , the dynamic behavior of the third particle stabilizes gradually, and the same happens when α and β increase together. It is hoped that the above results will help us to understand the dynamic evolution of binary system. Author Contributions: Formal analysis, R.W.; software, F.G.; writing—original draft, R.W.; writing—review and editing, F.G. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the National Natural Science Foundation of China (NSFC) though grant No.11672259, the China Scholarship Council through grant No. 201908320086. Acknowledgments: The authors thank the anonymous reviewers whose comments and suggestions helped improve and clarify this manuscript. Conflicts of Interest: The authors declare that there is no competing interests. Appendix A 2(1−µ)q1 2µq2 6(1−µ)q1 (2s1 −s2 ) 6µq (2s −s ) K1 = n2 β + r01 3  + r02 3 + + 2 r 53 4 ,  r01 5  02  5(2s −s ) 5(2s −s ) K2 = −3(1 − µ)q1 r 1 4 + r 1 6 2 V1 − 3µq2 r 1 4 + r 3 6 4 V2 , 01 01 02 02     3(1−µ)q1 1 5(4s1 −3s2 ) 3µq2 1 5(4s3 −3s4 ) K3 = 2 r01 4 + 2r01 6 V 1 + 2 r02 4 + 2r02 6 V2 ,     15(2s −s ) 15(2s −s ) K4 = 4(1 − µ)q1 r 1 5 + 2r 1 7 2 +4µq2 r 1 5 + 2r 3 7 4 , 01 01 02 02     15 ( 4s −3s ) 15 (4s −3s ) K5 = −6(1 − µ)q1 r 1 5 + 4r1 7 2 − 6µq2 r 1 5 + 4r3 7 4 , 01 01 02 02     3(1−µ)q1 1 5(4s1 −s2 ) 3µq2 1 5(4s3 −s4 ) K6 = 2 r01 4 + 2r01 6 V 1 + 2 r02 4 + 2r02 6 V2 ,     15(2s −s ) 15(2s −s ) K7 = 4(1 − µ)q1 r 1 5 + 2r 1 7 2 + 4µq2 r 1 5 + 2r 3 7 4 , 01 01 02 02     15(4s −3s ) 15(4s −3s ) K8 = −6(1 − µ)q1 r 1 5 + 4r1 7 2 − 6µq2 r 1 5 + 4r3 7 4 , 01 01 02 02     15(4s −s ) 15(4s −s ) K9 = −6(1 − µ)q1 r 1 5 + 4r 1 7 2 − 6µq2 r 1 5 + 4r 3 7 4 , 01 01 02 02 Universe 2020, 6, 35 20 of 22 (1−µ)q1 µq2 3(1−µ)q1 (4s1 −3s2 ) 3µq2 (4s3 −3s4 ) L1 = n2 β − r01 3 − r02 3 − − , 2r01 5  2r02 5  5(4s1 −3s2 ) 5 ( 4s −3s ) L2 = 3(1 − µ)q1 1 r01 4 + 2r 6 V1 + 3µq2 r 1 4 + 2r3 6 4 V2 , 01 02 02     15(4s −3s ) 15(4s −3s ) L3 = −6(1 − µ)q1 r 1 5 + 2r1 7 2 − 6µq2 r 1 5 + 2r3 7 4 , 01 01 02 02     3(1−µ)q1 1 15(6s1 −5s2 ) 3µq2 1 15(6s3 −5s4 ) L4 = 2 r01 5 + r01 7 + 2 r 5+ r02 7 , 02     3(1−µ)q1 1 5(6s −5s ) 3µq 5(6s −5s ) L5 = 2 r01 5 + 2r1 7 2 + 2 2 r 1 5 + 2r3 7 4 , 01 02 02     15(4s1 −3s2 ) 15(4s −3s ) L6 = −6(1 − µ)q1 r 5 + 4r 7 1 1 − 6µq2 r 5 + 4r3 7 4 , 01 01 02 02     3(1−µ)q1 1 15(2s1 −s2 ) 3µq2 1 15(2s3 −s4 ) L7 = 2 r 5 + 2r 7 + 2 r 5 + 2r 7 , 01 01 02 02 V1 = sgn(ξ − µ), V2 = sgn(ξ + 1 − µ), r01 = ξ − µ , r02 = ξ + 1 − µ , (1−µ)q1 µq 3(1−µ)q (4s −s ) 3µq (4s −s ) M1 = − − r 23 − r01 3 1 2r01 5 1 2 − 2 3 4 ,  02 2r02 5  5(4s −s ) 5(4s −s ) M2 = 3(1 − µ)q1 r 1 4 + 2r1 6 2 V1 + 3µq2 r 1 4 + 2r3 6 4 V2 , 01 01 02 02     3(1−µ)q1 1 5(6s1 −s2 ) 3µq2 1 5(6s3 −s4 ) M3 = 2 r01 5 + 2r 7 + 2 r 5 + 2r 7 , 01 02 02     15 ( 4s −s ) 15(4s −s ) M4 = −6(1 − µ)q1 r 1 5 + 4r 1 7 2 − 6µq2 r 1 5 + 4r 3 7 4 , 01 01 02 02     3(1−µ)q1 1 15(2s1 −s2 ) 3µq2 1 15(2s3 −s4 ) M5 = 2 r 5 + r 7 + 2 r 5+ r 7 . 01 01 02 02 Appendix B −2αnωK4 v05 = , 2(K1 L1 + 4(K1 + L1 − 4α2 n2 )ω2 + 16ω4 ) v12 = − L2αnω2, 1 +ω (K1 +4ω2 )L2 v12 +4αnω(K2 +K3 v12 2 +K4 ) v21 = − 2(K L +4(K +L −4α2 n2 )ω2 +16ω4 ) , 1 1 1 1 4αnωK2 −v12 [K1 L2 +4ω(αnK3 v12 +L2 ω)] v24 = 2K1 L1 +8(K1 +L1 −4α2 n2 )ω2 +32ω4 , M2 v32 = − 6ω2 , −2M2 (v32 +$21 )+M3 +M5 v12 2 −M4 v33 = 40ω2 , n  h i v34 = 1 4(K1 L1 +9(K1 +L1 −4α2 n2 )ω2 +81ω4 ) − K1 + 9ω2 2L2 (v21 + $21 v12 ) − L3 v12 3 + (L4 − L5 )v12 +6αnω[4(K2 $21 − K3 v12 v21 − K4 v32 ) + K5 − K6 − K7 ] , −6αnω(K3 v12 v24 − 14 K4 + 14 K5 v12 2 )− 12 (L2 v24 +L2 v12 $23 + 21 L3 v12 + 12 L4 v12 3 )(9ω2 +K1 ) v36 = K L +9(K +L −4α2 n2 )ω2 +81ω4 , 1 1 1 1 Universe 2020, 6, 35 21 of 22 K $02 = − 2K41 , K4 (L1 +4ω2 ) $03 = 2(K1 L1 +4(K1 +L1 −4α2 n2 )ω2 +16ω4 ) , K2 +K3 v12 2 +K4 $20 = − 2K1 , K1 +K3 v12 2 $21 = − 2K1 , −K2 (L1 +4ω2 )+v12 [4αnωL2 +K3 v12 (L1 +4ω2 )] $23 = 2K1 L1 +8(K1 +L1 −4α2 n2 )ω2 +32ω4 , 4αnωL2 v12 +(K2 +K3 v12 2 +K4 )(L1 +4ω2 ) $24 = 2(K1 L1 +4(K1 +L1 −4α2 n2 )ω2 +16ω4 ) , n h i $32 = 1 6αnω 2L2 (v21 + $21 v12 ) − L3 v12 3 + (L4 − L5 )v12 4(K1 L1 +9(K1 +L1 −4α2 n2 )ω2 +81ω4 )  o −[4(K2 $21 − K3 v12 v21 − K4 v32 ) + K5 − K6 − K7 ] L1 + 9ω2 ,      − 9ω2 + K1 K3 v12 v24 − 14 K4 + 14 K5 v12 2 + 3αnω L2 v24 + L2 v12 $23 + 12 L3 v12 + 12 L4 v12 3 $34 = . K1 L1 + 9(K1 + L1 − 4α2 n2 )ω2 + 81ω4 References 1. Musielak, Z.E.; Quarles, B. The three-body problem. Rep. Prog. Phys. 2014, 77, 30. [CrossRef] 2. Gao, F.B.; Zhang, W. A study on periodic solutions for the circular restricted three-body problem. Astron. J. 2014, 148, 116. [CrossRef] 3. Singh, J.; Begha, J.M. Periodic orbits in the generalized perturbed restricted three-body problem. Astrophys. Space Sci. 2011, 332, 319–324. [CrossRef] 4. Singh, J.; Begha, J.M. Stability of equilibrium points in the generalized perturbed restricted three-body problem. Astrophys. Space Sci. 2011, 331, 511–519. [CrossRef] 5. Tsirogiannis, G.A.; Douskos, C.N.; Perdios, E.A. Computation of the Lyapunov orbits in the photogravitational RTBP with oblateness. Astrophys. Space Sci. 2006, 305, 389–398. [CrossRef] 6. Singh, J.; Kalantonis, V.S.; Gyegwe, J.M.; Perdiou, A.E. Periodic motions around the collinear equilibrium points of the restricted three-body problem where the primary is a triaxial rigid body and secondary is an oblate spheroid. Astrophys. J. Suppl. 2016, 227, 1–13. [CrossRef] 7. Kalantonis, V.S.; Douskos, C.N.; Perdios, E.A. Numerical determination of homoclinic and heteroclinic orbits at collinear equilibria in the restricted three-body problem with oblateness. Celest. Mech. Dyn. Astron. 2006, 94, 135–153. [CrossRef] 8. Zotos, E.E. How does the oblateness coefficient influence the nature of orbits in the restricted three-body problem? Astrophys. Space Sci. 2015, 358, 1–33. [CrossRef] 9. Elshaboury, S.M.; Abouelmagd, E.I.; Kalantonis, V.S.; Perdios, E.A. The planar restricted three-body problem when both primaries are triaxial rigid bodies: Equilibrium points and periodic orbits. Astrophys. Space Sci. 2016, 361, 315. [CrossRef] 10. Jain, P.; Aggarwal, R.; Mittal, A.; Abdullah. Periodic orbits in the photogravitational restricted problem when the primaries are triaxial rigid bodies. Int. J. Astron. Astrophys. 2016, 6, 111–121. [CrossRef] 11. Idrisi, M.J.; Ullah, M.S. Non-collinear libration points in ER3BP with albedo effect and oblateness. J. Astrophys. Astron. 2018, 39, 28. [CrossRef] 12. Chen, J.J. A Bayesian Approach to the Understanding of Exoplanet Populations and the Origin of Life. Ph.D. Thesis, Columbia University, New York, NY, USA, 2018. 13. Hou, X.Y.; Xin, X.S.; Feng, J.L. Forced motions around triangular libration points by solar radiation pressure in a binary asteroid system. Astrodynamics 2019. [CrossRef] 14. Shi, Y.; Wang, Y.; Xu, S. Global search for periodic orbits in the irregular gravity field of a binary asteroid system. Acta Astronaut. 2019, 163, 11–23. [CrossRef] Universe 2020, 6, 35 22 of 22 15. Howard, A.; Fulton, B. Limits on planetary companions from doppler surveys of nearby stars. Publ. Astron. Soc. Pac. 2016, 128, 114401. [CrossRef] 16. Berardo, D.; Crossfield, I.J.M.; Werner, M.; Petigura, E.A.; Christiansen, J.; Ciardi, D.; Dressing, C.; Fulton, B.; Gorjian, V.; Greene, T.; et al. Revisiting the HIP 41378 System with K2 and Spitzer. Astron. J. 2019, 157, 185. [CrossRef] 17. Singh, J.; Perdiou, A.; Gyegwe, J.M.; Perdios, E.A. 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. Math. Comput. 2018, 325, 358–374. [CrossRef] 18. Das, M.K.; Narang, P.; Mahajan, S.; Yuasa, M. Effect of radiation on the stability of equilibrium points in the binary stellar systems: RW-Monocerotis, Krüger 60. Astrophys. Space Sci. 2008, 314, 261–274. [CrossRef] 19. Singh, J.; Perdiou, A.E.; Gyegwe, J.M.; Kalantonis, V.S. Periodic orbits around the collinear equilibrium points for binary Sirius, Procyon, Luhman 16, α-Centuari and Luyten 726-8 systems: The spatial case. J. Phys. Commun. 2017, 1, 025008. [CrossRef] 20. Singh, J.; Umar, A. On “out of plane” equilibrium points in the elliptic restricted three-body problem with radiating and oblate primaries. Astrophys. Space Sci. 2013, 344, 13–19. [CrossRef] 21. García, L.; Gómez, M. Optical polarization of solar type stars with debris disks. Rev. Mex. De Astron. Y Astrofísica 2015, 51, 1–10. 22. Abia, C.; Rebolo, R.; Beckma, J.E.; Crivellari, L. Abundances of light metals and Ni in a sample of disc stars. Astron. Astrophys. 1988, 206, 100–107. 23. Karaali, S.; Bilir, S.; Karatas, Y.; Ak, S.G. New metallicity calibration down to [Fe/H] = −2.75 dex. Publ. Astron. Soc. Aust. 2003, 20, 165–172. [CrossRef] 24. Perrin, M.N. Mass estimation of twelve K dwarfs. Symp. Int. Astron. Union 1976, 72, 167–172. [CrossRef] 25. Yumagulov, M.G.; Belikova, O.N.; Isanbaeva, N.R. Bifurcation near boundaries of regions of stability of libration points in the three-body problem. Astron. Rep. 2018, 62, 144–153. [CrossRef] 26. Perdomo, O.M. A bifurcation in the family of periodic orbits for the spatial isosceles 3 body problem. Qual. Theory Dyn. Syst. 2018, 17, 411–428. [CrossRef] 27. Maciejewski, A.; Rybicki, S. Global bifurcations of periodic solutions of the restricted three body problem. Celest. Mech. Dyn. Astron. 2004, 88, 293–324. [CrossRef] 28. Burdge, K.; Coughlin, M.W.; Fuller, J.; Kupfer, T.; Bellm, E.C.; Bildsten, L.; Graham, M.J.; Kaplan, D.L.; Van Roestel, J.; Dekany, R.G.; et al. General relativistic orbital decay in a seven-minute-orbital-period eclipsing binary system. Nature 2019, 571, 528–531. [CrossRef] 29. Antoniadou, K.; Libert, A.S. Spatial resonant periodic orbits in the restricted three-body problem. Mon. Not. R. Astron. Soc. 2018, 483, 2923–2940. [CrossRef] 30. Antoniadou, K.I.; Voyatzis, G. Resonant periodic orbits in the exoplanetary systems. Astrophys. Space Sci. 2014, 349, 657–676. [CrossRef] 31. Bosanac, N.; Howell, K.; Fischbach, E. Stability of orbits near large mass ratio binary systems. Celest. Mech. Dyn. Astron. 2015, 122, 27–52. [CrossRef] 32. Shampine, L.F.; Reichelt, M.W. The Matlab ode suite. SIAM J. Sci. Comput. 1997, 18, 1–22. [CrossRef] © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).