**1. Introduction**

In photogravitational celestial mechanics, along with the forces of Newtonian attraction *Fg*, the light pressure is taken into account *Fp*, coming from the radiating body (star) [1]. In some cases, the luminous flux is so intensive that the force *Fp* competes with gravity *Fg*, and can be even greater than that.

For a particular particle, the magnitude of the light pressure force depends not only on the power of the radiation source (star), but also on the cross-sectional area, the mass and the reflectivity of the particle. To determine the connection between the parameters of the star and the particle, a coefficient *Q* is introduced, called the particle mass reduction coefficient. For a particular particle, *Q* has a constant value that characterizes its susceptibility to radiation. The relationship between the parameters of the star [2] and the particle gives the reduction coefficient *Q*

$$Q = 1 - (1 - \varepsilon)A \frac{E}{fM} \tag{1}$$

(*f* is the gravitational parameter of the star, *E* and *M* is the mass and power of the star, *A* is a windage of the particle, determined by the ratio of the cross–sectional area to its mass, *ε* is the coefficient of light reflection). Sufficiently large and dense particles with small values of the parameters *A* and *ε* are most affected by the gravitational force of the star, therefore, *Q* > 0. For the smallest particles with high windage and reflection coefficient, the action of light is greater than gravity (*Q* < 0).

**Citation:** Tureshbaev, A.; Omarova, U.; Myrzayev, R. On the Stability of Collinear Libration Points in the Three-Body Problem with Two Radiating Masses. *Eng. Proc.* **2023**, *33*, 65. https://doi.org/10.3390/ engproc2023033065

Academic Editors: Askhat Diveev, Ivan Zelinka, Arutun Avetisyan and Alexander Ilin

Published: 22 August 2023

**Copyright:** © 2023 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 (https:// creativecommons.org/licenses/by/ 4.0/).

The photogravitational three-body problem introduced by V.V. Radzievskiy [3] has become a good dynamic model for studying the motion of microparticles in binary star systems.

In the elliptical version of the problem (when the orbits of a stellar pair are elliptical), we write the equations of motion of a particle in a rectangular system rotating together with the stars in the form

$$
\ddot{\mathbf{x}} - 2\dot{\mathbf{y}} = \frac{\partial W}{\partial \mathbf{x}},
\ddot{\mathbf{y}} + 2\dot{\mathbf{x}} = \frac{\partial W}{\partial \mathbf{y}},
\ddot{z} = \frac{\partial W}{\partial z} \tag{2}
$$

where *x*, *y*, *z* are dimensionless coordinates related to the distance *r* = *p*/(1 + *e* cos *v*) − 1 between the stars (*p* and *e* are the focal parameter and the eccentricity of the relative orbital motion of the stellar pair) that are the rectangular coordinates of the particle and *W* is the force function of the system, equal to

$$\begin{aligned} \mathcal{W} &= \left(1 + e \cos v\right)^{-1} \left[\left(x^2 + y^2 - z^2 e \cos v\right)/2 + Q\_1(1 - \mu)/R\_1 + Q\_2 \mu/R\_2\right] \\\\ \mathcal{R}\_1 &= \left[\left(x + \mu\right)^2 + y^2 + z^2\right]^{1/2}, \mathcal{R}\_2 = \left[\left(x - 1 + \mu\right)^2 + y^2 + z^2\right]^{1/2} \end{aligned} \tag{3}$$

Here, *μ* and 1 − *μ* are the dimensionless masses of the stars, and *Q*<sup>1</sup> and *Q*<sup>2</sup> are the reduction coefficients of their mass, which represent the ratio of the difference between the gravitational and repulsive forces to the gravitational force. In terms of physical meaning, numerical values *Q*<sup>1</sup> and *Q*<sup>2</sup> do not exceed 1. For the classical problem *Q*<sup>1</sup> = 1, *Q*<sup>2</sup> = 1 (no radiation of bodies) [4].

Libration points—constant solutions of the adopted system of dynamic equations represent relative equilibria in a circular problem and periodic motions in an elliptic problem. They are found from the system of equations

$$\frac{\partial \mathcal{W}}{\partial \mathbf{x}} = 0, \frac{\partial \mathcal{W}}{\partial y} = 0, \frac{\partial \mathcal{W}}{\partial z} = 0,\tag{4}$$

CLP are located on a straight line connecting the main bodies, and for them *y* = 0, *z* = 0. Their positions on the abscissa axis are determined from the first equation of the system (4). The triangular libration points were carefully studied in [5]. The stability of triangular points in strictly nonlinear formulation were considered in [6,7]. In [8,9], the nonlinear analysis of stable coplanar libration points that are not on the plane of orbital motion of main bodies was completed.

#### **2. Collinear Libration Points and Their Stability in a Plane Problem**

The CLP coordinate is determined from the following equation

$$f(\mathbf{x}) = \mathbf{x} - Q\_1 \frac{(1 - \mu)(\mathbf{x} + \mu)}{|\mathbf{x} + \mu|^3} - Q\_2 \frac{\mu(\mathbf{x} + \mu - 1)}{|\mathbf{x} + \mu - 1|^3} = 0 \tag{5}$$

Equation (5) contains three mutually independent parameters *μ*, *Q*1, *Q*2. Therefore, the problem can consider a three-parameter family of CLPs. Let us consider the case when both components of the binary star radiate and move in circular orbits. From Equation (5), we obtain

$$f'(\mathbf{x}) = 1 + 2a(\mathbf{x}),\\a = Q\_1 \frac{1 - \mu}{|\mathbf{x} + \mu|^3} + Q\_2 \frac{\mu}{|\mathbf{x} + \mu - 1|^3}$$

Thus, there is a parameter *a* [10], which is included in the characteristic equation

$$
\lambda^4 + (2 - a)\lambda^2 - (1 - a)(1 + 2a) = 0,\tag{6}
$$

its solutions are equal to

$$
\lambda\_a^2 = \frac{1}{2}(a-2 \pm \sqrt{(9a-8)a}, a = 1,2)
$$

When changing the parameter values *a* in the intervals 8/9 < *a* ≤ 1 and −0.5 < *a* ≤ 0, roots *λα* are purely imaginary, and at the boundaries there are multiples. Therefore, the parameter values indicated above by the inequalities correspond to the stability region of libration points in the first approximation. In [11], they constructed diagrams of CLP stability.

In the problem, resonances of the third and fourth orders are found; for third-order resonances, the resonant values of the parameter *a* have the form *a*∗ <sup>±</sup> <sup>=</sup> 41/108 <sup>±</sup> <sup>5</sup> <sup>√</sup>145/108; for fourth-order resonances, *<sup>a</sup>*<sup>±</sup> = (<sup>68</sup> <sup>+</sup> <sup>60</sup>√5)/209. As expected, the resonance of the third order leads to the instability of the CLP [12].

In [13], it is shown that at the resonance of the fourth order, the CLPs are stable by Lyapunov.

Below, we consider the stability of the CLP in the spatial problem. The Birkhoff normal form is used, and the Arnold–Moser theorem [14] is applied.

#### **3. Equations of Motion and Expansion of the Hamilton Function**

Particle motion *P*(*x*, *y*, *z*) is given by the canonical equations

$$\frac{d\overline{q\_i}}{dt} = \frac{\partial H}{\partial \overline{p\_i}}, \frac{d\overline{p\_i}}{dt} = -\frac{\partial H}{\partial \overline{q\_i}}, (i = 1, 2, 3) \tag{7}$$

where *qi* are the Cartesian coordinates of the particle *P*(*x*, *y*, *z*), *pi* are the corresponding canonical momenta and *H*(*x*, *y*, *z*, *p*1, *p*2, *p*3) is the Hamilton analytic function with respect to coordinates and momenta, which, in our case, has the form

$$H = \frac{1}{2} (\overline{p}\_1^2 + \overline{p}\_2^2 + \overline{p}\_3^2) + (\overline{p}\_1 y - \overline{p}\_2 x) - Q\_1 (1 - \mu) / R\_1 - Q\_2 \mu / R\_2$$

$$R\_\alpha = \sqrt{(\mathbf{x} - \mathbf{x}\_\alpha)^2} + y^2 + z^2 \text{ ( $\alpha = 1, 2$ )}$$

Here, *Q*<sup>1</sup> and *Q*<sup>2</sup> are coefficients of reduction of the masses of the main bodies, which, in the case of CLP, can take both positive and negative values [10].

We study the stability of the CTL under the assumption that the orbit of the main bodies is circular and the particle *P* of infinitely small mass at the initial moment of time experiences initial perturbations that take it out of the plane of rotation of the main bodies *S*<sup>1</sup> and *S*2.

We introduce perturbations into Equation (1) according to the formulas

$$\begin{aligned} \mathbf{x} &= \mathbf{x}^\* + q\_1, y = q\_2, z = q\_3, \overline{p}\_1 = \overline{p}\_1^\* + p\_1, \overline{p}\_2 = p\_2, \overline{p}\_3 = p\_3 \\\ p\_1^\* &= \mathbf{x}^\*, p\_2^\* = y^\* = p\_3^\* = z\_0^\* = 0, \end{aligned} \tag{9}$$

where

$$\alpha^\* = 0.5(Q\_1^{\frac{2}{3}} - Q\_2^{\frac{2}{5}} - 1) - \mu,\\ \overline{p}\_1^\* = \mp 0.5\sqrt{2(Q\_1^{\frac{2}{3}} + Q\_2^{\frac{2}{5}}) - (Q\_1^{\frac{2}{5}} - Q\_2^{\frac{2}{5}}) - 1}.\tag{10}$$

Expanding the Hamilton function in a row according to the degrees of perturbations and in the neighborhood of the considered collinear point, taken as the origin, we obtain

$$H = H\_2 + H\_3 + H\_4 + \dots \tag{11}$$

Here, *Hm* are homogeneous polynomials of degree *m*(*m* = 2, 3, 4, ...) with respect to generalized coordinates and impulses *pi*, so

$$H\_m = \sum\_{v+l=m} h\_{v\_1 v\_2 v\_3 l\_1 l\_2 l\_3} \cdot q\_1^{v\_1} q\_2^{v\_2} q\_3^{v\_3} p\_1^{l\_1} p\_2^{l\_2} p\_3^{l\_3} \tag{12}$$

Then, in Equation (11), the forms *H*2, *H*<sup>3</sup> and *H*<sup>4</sup> taking into account (9) will take the following form:

$$H\_2 = \frac{1}{2}(p\_1^2 + p\_2^2 + p\_3^2) + p\_1q\_2 - p\_2q\_1 + h\_{200}q\_1^2 + h\_{020}q\_2^2 + h\_{002}q\_3^2 + h\_{110}q\_1q\_2 \tag{13}$$

$$+ h\_{101}q\_1q\_3 + h\_{001}q\_2q\_3$$

$$H\_3 = h\_{300}q\_1^3 + h\_{030}q\_2^3 + h\_{300}q\_3^3 + h\_{210}q\_1^2q\_2 + h\_{201}q\_1^2q\_3 + h\_{120}q\_1q\_2^2 + h\_{021q\_2^2q\_3} \tag{14}$$

$$+ h\_{102}q\_1q\_3^2 + h\_{012}q\_2q\_3^2 + h\_{111}q\_1q\_2q\_3.$$

$$\begin{aligned} H\_4 &= h\_{400}q\_1^4 + h\_{040}q\_2^4 + h\_{004}q\_3^4 + h\_{310}q\_1^3q\_2 + h\_{130}q\_1q\_2^3 \\ &+ h\_{103}q\_1q\_3^3 + h\_{301}q\_1^3q\_3 + h\_{031}q\_2^3q\_3 + h\_{013}q\_3^3q\_2 + h\_{211}q\_1^2q\_2q\_3 \\ &+ h\_{121}q\_1q\_2^2q\_3 + h\_{112}q\_1q\_2q\_3^2 + h\_{220}q\_1^2q\_2^2 + h\_{202}q\_1^2q\_3^2 + h\_{022}q\_2^2q\_3^2 \end{aligned} \tag{15}$$

where

$$\begin{aligned} h\_{200} &= -8a, h\_{020} = 4a, h\_{002} = 4a, h\_{101} = 0, h\_{101} = 0, h\_{011} = 0 \\ h\_{300} &= 16b, h\_{120} = -16b, h\_{102} = -16b, h\_{030} = 0, h\_{003} = 0, h\_{210} = 0, \\ h\_{201} &= 0, h\_{021} = 0, h\_{012} = 0, h\_{111} = 0, \\ h\_{400} &= -32c, h\_{040} = -12c, h\_{004} = -12c, h\_{220} = 32c, h\_{202} = 32c, \\ h\_{022} &= -8c, h\_{310} = 0, h\_{130} = 0, h\_{103} = 0, h\_{301} = 0, h\_{031} = 0, \\ h\_{013 \to 0}, h\_{211} &= 0, h\_{121} = 0, h\_{112} = 0. \end{aligned} \tag{16}$$

$$\begin{aligned} a &= \frac{Q\_1(1-\mu)}{|Q\_1^{2/3} - Q\_2^{2/3} + 1|^3} + \frac{Q\_2\mu}{|Q\_1^{2/3} - Q\_2^{2/3} - 1|^3}, \\ b &= \frac{Q\_1(1-\mu)(Q\_1^{2/3} - Q\_2^{2/3} + 1)}{|Q\_1^{2/3} - Q\_2^{2/3} + 1|^5} + \frac{Q\_2\mu(Q\_1^{2/3} - Q\_2^{2/3} - 1)}{|Q\_1^{2/3} - Q\_2^{2/3} - 1|^5} \\ c &= \frac{Q\_1(1-\mu)}{|Q\_1^{2/3} - Q\_2^{2/3} + 1|^5} + \frac{Q\_2\mu}{|Q\_1^{2/3} - Q\_2^{2/3} - 1|^5} \end{aligned} \tag{17}$$

#### **4. Stability of CLP in a Spatial Problem**

The question of the stability of the investigated spatial CLPs can be considered as a stability problem of equilibrium positions *qi* = *pi* = 0(*i* = 1, 2, 3) of an autonomous Hamiltonian system with three degrees of freedom. As can be seen from (13), here we have the case when *H*<sup>2</sup> is not a sign-definite function, and the characteristic equation of the system has no roots with a nonzero real part. Hence, the stability of the complete system does not follow from the stability of a linear system.

Expanding the Hamilton function into a power series *qi*, *pi* in the vicinity of the considered equilibrium position, first the Hamiltonian *H*<sup>2</sup> is transformed to the normal form in the form

$$K2 = \omega\_1 r\_1 - \omega\_2 r\_2 + \omega\_3 r\_3. \tag{18}$$

The structure of the normal form depends on the type of resonance relation

$$
\omega\_1 r\_1 - \omega\_2 r\_2 + \omega\_3 r\_3 = 0 (|k\_1| + k\_2 + |k\_3| \le 4),
\tag{19}
$$

where the frequencies of the principal oscillations for the libration points are equal to

$$
\omega\_1 = \sqrt{(2 - a + \sqrt{(9a - 8)a}) / 2}, \\
\omega\_2 = \sqrt{(2 - a - \sqrt{(9a - 8)a}) / 2}, \\
\omega\_3 = \sqrt{a}. \tag{20}
$$

As can be seen from the last expression, for the frequency of spatial oscillations, the parameter *a* can take only positive values. Therefore, resonances containing the frequency of spatial oscillations can be realized only in a limited part of the region (8/9 < *a* ≤ 1 and −1/2 < *a* ≤ 0) for necessary stability conditions of the system. Let us investigate the stability of CLP at two-frequency resonances. For CLP, the following two-frequency resonances turned out to be possible:

$$
\omega\_1 = 2\omega\_2, \omega\_1 = 3\omega\_2, 2\omega\_1 = \omega\_3, 3\omega\_1 = \omega\_3, 2\omega\_2 = \omega\_3, 3\omega\_2 = \omega\_3
$$

Resonances *ω*<sup>1</sup> = 2*ω*<sup>2</sup> and *ω*<sup>1</sup> = 3*ω*2, discovered in the plane problem, were studied in [12,13]. In the spatial photogravitational problem, resonances of the third and fourth orders turned out to be possible

$$2\omega\_1 = \omega\_3, 3\omega\_1 = \omega\_3, 2\omega\_2 = \omega\_3, 3\omega\_2 = \omega\_3$$

which respectively correspond to the values of the parameter *a* defined as

$$a = 4(1 + 2\sqrt{7})/27,\\ a = 4(-1 + \sqrt{10})/9,\\ a = (63 + \sqrt{53217})/304,\\ a = (63 + \sqrt{53217})/304$$

Note that the last two resonances, 3*ω*<sup>1</sup> = *ω*<sup>3</sup> and 3*ω*<sup>2</sup> = *ω*3, match. To construct resonance curves (in the stability region in the linear approximation of the system) for the corresponding specific resonance value of the coefficient *a*, a curve is constructed, which is determined by the expression

$$\frac{Q\_1(1-\mu)}{|Q\_1^{2/3}-Q\_2^{2/3}+1|^3} + \frac{Q\_2\mu}{|Q\_1^{2/3}-Q\_2^{2/3}-1|^3} = a$$

At resonance 2*ω*<sup>1</sup> = *ω*<sup>3</sup> (which does not involve the frequency of plane oscillations), which corresponds to the value of the parameter *a* = 4(1 + 2 <sup>√</sup>7)/27, the normalized Hamiltonian takes the form [14]

$$H = 2\omega\_1 r\_1 - \omega\_1 r\_3 + A(\omega\_1, \omega\_3) r\_3 \sqrt{r\_1} \sin(\varphi\_1 + 2\varphi\_3) + O(r\_1 + r\_3)^2,\tag{21}$$

where *A*(*ω*1, *ω*3) = − *ω*1(*x*<sup>2</sup> <sup>1002</sup> + *<sup>y</sup>*<sup>2</sup> <sup>1002</sup>) and the coefficients *x*<sup>1002</sup> and *y*<sup>1002</sup> look like

$$\mathcal{X}\_{1002} = -\frac{\omega\_1 h\_{0111}}{2\omega\_1} - \frac{h\_{1002}}{2} + \frac{h\_{1200}}{2\omega\_1^2}, \\ \mathcal{Y}\_{1002} = -\frac{\omega\_1 h\_{012}}{2} - \frac{\omega\_1 h\_{0210}}{2\omega\_1^2} + \frac{h\_{1101}}{2\omega\_1}.$$

which, for collinear points, are equal to

$$
\alpha\_{1002} = -\frac{\hbar\_{1200}}{2\omega\_1^2}, y\_{1002} = 0\tag{22}
$$

hence, the expression

$$A(\omega\_1, \omega\_3) = -\sqrt{\omega\_1(\mathbf{x}\_{1002}^2 + y\_{1002}^2)} = -\sqrt{\omega\_1}\mathbf{x}\_{1002}$$

is not equal to zero anywhere; therefore, according to the Arnold–Moser theorem, at a third-order resonance from the stability region, in the first approximation, we can say that the CTLs are unstable. If there is a fourth-order resonance in the system, corresponding to the value of the parameter *<sup>a</sup>* = (<sup>63</sup> <sup>+</sup> <sup>√</sup>53217)/304, using the Birkhoff transformation in the original Hamiltonian, we annihilate the terms of the third degree. The Hamiltonian normalized in this case in polar coordinates will take the following form [13]:

$$\begin{aligned} H &= 3\omega\_1 r\_1 - \omega\_1 r\_3 + c\_{20} r\_1^2 + c\_{11} r\_1 r\_3 + c\_{02} r\_3^2 + \\ B(\omega\_1, \omega\_3) r\_3 \sqrt{r\_1 r\_3} \cos(\varphi\_1 + 3\varphi\_3) + O(r\_1 + r\_3)^{5/2} \end{aligned} \tag{23}$$

where

$$B(\omega\_1, \omega\_3) = \frac{1}{3}\omega\_3\sqrt{3(\mathbf{x}\_{1003}^2 + y\_{1003}^2)}$$

It is important to notice that in the classical problem for a fixed value *μ*, the coefficients *B*(*ω*1, *ω*3), *c*200, *c*<sup>100</sup> and *c*<sup>020</sup> take constant values (which simplifies the investigation of the problem). In this problem, the same coefficients do not remain constant and are functions of arbitrary coefficients of the coefficients *Q*<sup>1</sup> and *Q*2, which makes the task much more difficult. These are denoted by the coefficients of the Hamiltonian (23)

$$N\_1 = c\_{200} + 3c\_{110} + 9c\_{020}, \\ N\_2 = 3\sqrt{3}B(\omega\_1, \omega\_3)\_r$$

where

$$B(\omega\_1, \omega\_3) = \frac{1}{3}\omega\_3\sqrt{3(\varkappa\_{1003}^2 + \jmath\_{1003}^2)}$$

is defined by expressions

$$\begin{split} x\_{1003} &= \frac{1}{2}\omega\_1 h\_{0013} + \frac{1}{2\omega\_3^3}h\_{1300} - \frac{1}{2\omega\_3}h\_{102} - \frac{\omega\_1}{2\omega\_3^2}h\_{0211} \\ -\frac{9}{5}(\mathbf{x}\_{0120} \ast \mathbf{x}\_{0012} + y\_{0120} \ast y\_{0012}) - \frac{1}{\omega\_3}(\mathbf{x}\_{1002}y\_{1011} + \mathbf{x}\_{1011}y\_{1002}) \\ &+ \frac{4}{\omega\_3^2}(\mathbf{x}\_{1002}\mathbf{x}\_{0201} + y\_{1002}y\_{0201}) + \frac{3}{2}(\mathbf{x}\_{0003}\mathbf{x}\_{0111} + y\_{0003}y\_{0111}), \\ y\_{1003} &= -\frac{\omega\_1}{2\omega\_3}h\_{0112} + \frac{1}{2}h\_{1003} + \frac{1}{2\omega\_3^2}h\_{1201} + \frac{\omega\_1}{2\omega\_3^3}h\_{0310} \\ -\frac{9}{5}(\mathbf{x}\_{0120}\mathbf{x}\_{0012} + \mathbf{x}\_{0012}\mathbf{x}\_{0120}) - \frac{1}{\omega\_3}(y\_{1011}y\_{1002} - \mathbf{x}\_{1011}\mathbf{x}\_{1002}) \\ &+ \frac{4}{\omega\_3^2}(\mathbf{x}\_{0201}y\_{1002} + \mathbf{x}\_{1002}y\_{0201}) + \frac{3}{2}(\mathbf{x}\_{0111}\mathbf{x}\_{0003} + \mathbf{x}\_{0003}y\_{0111}), \end{split}$$

where coefficients

*h*0013, *h*1300, *h*1102, *h*0211, *h*0112, *h*1003, *h*1201, *h*0310, *x*0120, *y*0120, *x*1011, *y*1002, *y*1011, *x*0012, *y*0111, *x*0201, *y*0201, *x*0003, *x*<sup>0003</sup>

given for CLP above (16) take values equal to zero, therefore, they are identically equal to zero *x*<sup>1003</sup> and *y*<sup>1003</sup> . Then, the equality takes place

$$N\_2 = 3\sqrt{3}B(\omega\_1, \omega\_3) = 0.$$

Let us now define the value *N*<sup>1</sup> = *c*<sup>200</sup> + 3*c*<sup>110</sup> + 9*c*020. Here, the coefficients *c*200, *c*110*c*020, which are invariants of the Hamilton function (11) with respect to canonical transformations, depend on the coefficients *hv*1*v*2*l*1*l*<sup>2</sup> that are homogeneous polynomials (12) of degree *m*(*m* = 3, 4), which are equal to

$$\begin{split} c\_{200} &= \frac{3}{2\omega\_1^2} h\_{4000} - \frac{27}{8} \omega\_2^2 y\_{0030}^2 - \frac{3}{2} \mathbf{x}\_{1020}^2 \\ c\_{110} &= \frac{1}{\omega\_1 \omega\_2} h\_{2200} - \frac{2}{3} \mathbf{x}\_{1002}^2 + \frac{3}{10} \omega\_2^2 y\_{0012}^2 + 2 \mathbf{x}\_{0111} \mathbf{x}\_{1020} \\ c\_{020} &= \frac{3}{2\omega\_2^2} h\_{0400} - \frac{1}{6} \mathbf{x}\_{0111}^2 - \frac{3}{40} \omega\_2^2 y\_{0012}^2 \end{split} \tag{24}$$

where

$$\begin{split} y\_{0030} &= -\frac{1}{\omega\_1^2} h\_{3000}, \mathbf{x}\_{1020} = -\frac{3}{2\omega\_1^2} h\_{3000}, \mathbf{x}\_{1002} = \frac{1}{2\omega\_2^2} h\_{1200} \\ y\_{0012} &= \frac{1}{\omega\_1 \omega\_2^2} h\_{1200}, \mathbf{x}\_{0111} = \frac{1}{\omega\_1 \omega\_2} h\_{1200} \end{split} \tag{25}$$

Substituting from (16) the values *h*<sup>3000</sup> = 16*b*, *h*<sup>1200</sup> = −16*b* in (24), we have

$$\begin{split} c\_{200} &= -\frac{48}{\omega\_1^2}c - \frac{864\omega\_2^2}{\omega\_1^2}(1+\frac{1}{\omega\_1^2})^2b^2 - 96(1-\frac{3}{\omega\_1^2})^2b^2, \\ c\_{110} &= \frac{32}{\omega\_1\omega\_2}c - \frac{96}{\omega\_2^4}b\_2 + \frac{384}{5\omega\_1^2\omega\_2^2}b^2 - \frac{256}{\omega\_1\omega\_2}(1-\frac{3}{\omega\_1^2})b^2, \\ c\_{020} &= -\frac{18}{\omega\_2^2}c - \frac{32}{3\omega\_2^4}b^2 - \frac{96}{5\omega\_1^2\omega\_2^2b^2} \end{split} \tag{26}$$

where *a*, *b* and *c* are parameters that depend on reduction factors *Q*<sup>1</sup> and *Q*<sup>2</sup> and dimensionless mass parameter *μ*. As the calculations showed, the modulus of the expression *N*<sup>1</sup> = *c*<sup>200</sup> + 3*c*<sup>110</sup> + 9*c*<sup>020</sup> is always different from zero. Consequently, the inequality holds everywhere |*N*1| > *N*<sup>2</sup> = 0, which, according to [13], guarantees the existence of Lyapunov stability. In a similar way, it is proved that at a resonance of the third order 2*ω*<sup>2</sup> = *ω*3, CLPs are unstable, and at a fourth-order resonance 3*ω*<sup>2</sup> = *ω*<sup>3</sup> they are stable by Lyapunov. Below are the regions of the stability of the linear system (colored), in which the resonance curves of the fourth order are indicated 2*ω*<sup>2</sup> = *ω*<sup>3</sup> for two values of the mass parameter *μ*.

At *μ* = 0.001, the resonance curve is located closer to the middle of the stability region (Figure 1a). When the mass *μ* increases up to 0.01 (Figure 1b), the region becomes slightly smaller and the resonance curve becomes closer to the boundary of the stability region and becomes less noticeable than when *μ* = 0.001. Apparently, this fact confirms the above conclusion that resonances containing the frequency of spatial oscillations can be realized only in a limited part (8/9 < *a*), being the area of necessary conditions for the stability of CLP.

**Figure 1.** (**a**) Resonance curve of the fourth at *μ* = 0.001; (**b**) Resonance curve of the fourth at *μ* = 0.01.

Note that in the classical problem for a fixed value *μ*, the coefficients *c*200, *c*110, *c*<sup>020</sup> are constants (which much simplifies the investigation of the problem). However, in this problem, the same coefficients are not constants but functions of the coefficients *Q*<sup>1</sup> and *Q*2, which makes the task much more difficult.

If a *ω<sup>i</sup>* does not satisfy condition (19), then after applying the Birkhoff transformation, the Hamiltonian of the perturbed motion in polar coordinates normalized to the fourth order inclusive has the form

$$H^\* = K\_2(r\_{1\prime}, r\_{2\prime}r\_3) + K\_4(r\_{1\prime}r\_{2\prime}r\_3) \tag{27}$$

Here, *K*<sup>4</sup> is defined by the expression

$$K\_4 = c\_{200}r\_1^2 + c\_{110}r\_1r\_3 + c\_0 11r\_2r\_3 + c\_{002}r\_3^2 \tag{28}$$

Now, we use Arnold's results on the stability of Hamiltonian systems for most of the initial conditions [13]. It is known that the instability found in the plane problem remains such in the spatial problem.

Assuming that there are no resonances in the system 2*ω*<sup>1</sup> = *ω*3, *ω*<sup>1</sup> = 2*ω*2, *ω*<sup>1</sup> = 3*ω*2, 3*ω*<sup>1</sup> = *ω*3, 2*ω*<sup>2</sup> = *ω*3, 3*ω*<sup>2</sup> = *ω*3, consider a fourth-order determinant

$$D\_4 = \det \begin{vmatrix} \frac{\partial^2 K\_4}{\partial r\_i \partial r\_j} & \frac{\partial K\_2}{\partial r\_i} \\ \frac{\partial K\_2}{\partial r\_j} & 0 \end{vmatrix} \tag{29}$$

Expanding the determinant (29), we have

$$\begin{aligned} D\_4 &= \omega\_1^2 (c\_{011}^2 - 4c\_{020}c\_{002} + \omega\_2^2 (c\_{101^2} - 4c\_{200}c\_{002}) + \\ &+ \omega\_3^2 (c\_{110}^2 - 4c\_{200}c\_{020}) + 2\omega\_1 \omega\_2 (c\_{101}c\_{011} - 2c\_{002}c\_{110}) - \\ &- 2\omega\_1 \omega\_3 (c\_{001}c\_{110} - 2c\_{020}c\_{101}) + 2\omega\_2 \omega\_3 (c\_{110}c\_{101} - 2c\_{200}c\_{011}) \end{aligned} \tag{30}$$

The balance position *qi* = *pi* = 0 is stable for most initial conditions (by the Lebesgue measure) when the determinant *D*<sup>4</sup> = 0. After using numerical analysis, we check the validity of the inequality *D*<sup>4</sup> = 0. We see that in the spatial photogravitational three-body problem, the collinear libration points are stable for most initial conditions (by the Lebesgue measure) for all *a* (except for the values corresponding to internal resonances of the third 2*ω*<sup>1</sup> = *ω*<sup>3</sup> and 2*ω*<sup>2</sup> = *ω*<sup>3</sup> and fourth *ω*<sup>1</sup> = 3*ω*2, 21 = 3*ω*3, 3*ω*<sup>2</sup> = *ω*<sup>3</sup> orders) from the stability region in the linear approximation.

The presence of stability in the system for most of the initial conditions means, with a probability close to unity, that the KTLs are stable in the spatial problem.

As shown by numerical calculations, CLPs are formally stable for almost all values of the parameters from the stability region in the linear approximation. The exceptions are, in addition to the values of the parameters corresponding to the studied resonance, perhaps those values *μ*, *Q*1, *Q*<sup>2</sup> from the stability region, at which resonances above the fourth order are realized.

The presence of formal stability means that Lyapunov instability is not detected over a practically very long time interval. This suggests that the particles will stay near the stable libration points for quite a long time.

**Author Contributions:** Conceptualization, A.T.; methodology, A.T. and U.O.; software, R.M. and U.O.; alidation, A.T.; writing—rough preparation, R.M. and U.O.; writing—review and editing, A.T.; All authors have read and agreed with the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors have declared that there are no conflict of interest.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
