Secular perturbations of quasi-elliptic orbits in the restricted three-body problem with variable masses

International Journal of Non-Linear Mechanics 73 (2015) 58–63

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Secular perturbations of quasi-elliptic orbits in the restricted three-body problem with variable masses$ A.N. Prokopenya a,b,n, M.Zh. Minglibayev c,d, B.A. Beketauov c a

Warsaw University of Life Sciences – SGGW, Nowoursynowska Str. 159, 02-776 Warsaw, Poland Collegium Mazovia Innovative Higher School, Sokolowska Str. 161, 08-110 Siedlce, Poland c Al-Farabi Kazakh National University, Al-Farabi Ave. 71, Almaty 050038, Kazakhstan d Fessenkov Astrophysical Institute, Observatoriya 23, Almaty 050020, Kazakhstan b

art ic l e i nf o

a b s t r a c t

Article history: Received 14 February 2014 Received in revised form 21 August 2014 Accepted 3 November 2014 Available online 13 November 2014

In this work we consider the satellite version of the restricted three-body problem when masses of the primary bodies P0, P1 vary isotropically with different rates, and their total mass changes according to the joint Meshcherskii law. Equations of motion of the body P2 of infinitesimal mass are obtained in terms of the osculating elements of the aperiodic quasi-conical motion about the body P0. Doubly averaging these equations and using the Hill approximation, we have obtained the differential equations determining the secular perturbations of the orbital elements and determined the domains of possible values of the system parameters for which their analytical solutions are expressed in terms of elementary or elliptic functions. The bodies P0, P1 mass variation laws for which the corresponding differential equations are integrable, have been found. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Restricted problem of three bodies Variable masses Aperiodic quasi-conical motion Secular perturbations

1. Introduction The restricted three-body problem is a well-known model of celestial mechanics, having a long history (see [1]), it was first considered by Euler (1772) and Jacobi (1836). Many celestial systems can be reduced in a first approximation to this model, and due to this it remains a popular subject of research for mathematicians, mechanicians and astronomers. Besides, a general solution of this problem has not been found, in spite of its apparent simplicity, and this stimulates further research of the model and its extensions and generalizations (see, for instance, [2, Chapter 4.7]). Remind that in the simplest case the restricted three-body problem describes a motion of the point P2 of negligible mass in the gravitational field of two massive points P0 and P1, moving in Keplerian orbits about their common center of mass. It is assumed that the masses of points P0 and P1 are given and their orbits are completely determined by the known solution of the two-body problem. Although this problem is not integrable, some of its exact

☆ This work was partly supported by the Grant 0688/GF of the Scientific-Technical Programs and projects of the Committee of Science of the Republic of Kazakhstan, 2012–2014. n Corresponding author at: Warsaw University of Life Sciences – SGGW, Nowoursynowska Str. 159, 02-776 Warsaw, Poland. E-mail addresses: [email protected] (A.N. Prokopenya), [email protected] (M.Zh. Minglibayev), [email protected] (B.A. Beketauov).

http://dx.doi.org/10.1016/j.ijnonlinmec.2014.11.007 0020-7462/& 2014 Elsevier Ltd. All rights reserved.

particular solutions can be found and their stability can be studied; such investigations contributed substantially to development of the theory of non-linear Hamiltonian systems and their stability (see, for example, [3]). Note that real-life celestial bodies are not stationary; their parameters such as mass, size, shape, and internal structure may change with time (see [4,5]). Taking into consideration possible changes of one parameter only, for example, the masses of points P0 and P1, considerably complicates the problem because an exact solution of the corresponding two-body problem can be found in an analytical form only for special cases. In the present paper we assume that the masses m0 ðtÞ and m1 ðtÞ of points P0 and P1, respectively, vary isotropically with different rates, but their total mass reduces according to the joint Meshcherskii law m00 þ m10 ¼ m0 ðtÞ þ m1 ðtÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi At 2 þ 2Bt þC  vðtÞ;

ð1Þ

where m00 ¼ m0 ðt 0 Þ, m10 ¼ m1 ðt 0 Þ, t0 is an initial instant of time, and A; B; C are parameters chosen in a way to satisfy the condition vðt 0 Þ ¼ 1 and v(t) to be an increasing function for t 4 t 0 . In this case equation of the points P0, P1 relative motion is reduced to ordinary equation of Keplerian motion for constant masses by means of a scale transformation of the spatial coordinates and time (see [6]). Such equation is integrable and one of its exact particular solution describes a uniform motion in a circle. Therefore, using the relative coordinates system with origin at point P0 and assuming, without losing generality, that trajectory of the point P1 is situated in the

A.N. Prokopenya et al. / International Journal of Non-Linear Mechanics 73 (2015) 58–63

2. Equations of motion in Delaunay's elements

coordinate plane XOY, one can describe the point P1 motion as XðtÞ ¼ vðtÞa1 cos θ1 ðtÞ;

YðtÞ ¼ vðtÞa1 sin θ1 ðtÞ;

ZðtÞ ¼ 0;

ð2Þ

where the function θ1 ðtÞ is determined by differential equation dθ 1 _  θ1 ¼ dt

AC  B2 þ

K a31

!1=2

1 ; v2 ðtÞ

59

ð3Þ

positive parameter a1 is determined from the initial conditions, K ¼ f ðm00 þ m10 Þ, and f is the constant of gravitation. Note that in case of constant masses (A ¼ B ¼ 0; C ¼ 1) solution (2) determines a uniform motion of point P1 in a circular orbit of radius a1 with angular velocity ω1 ¼ ðK=a31 Þ1=2 . If solution of the two-body problem with variable masses is given, one can formulate the restricted three-body problem, similarly to the case of constant masses. As this problem is also not integrable, two main approaches to its investigation are usually applied. First, one can try to find exact particular solutions and investigate the effect of variations of the points masses on these solutions and their stability; different aspects of this problem are considered in a number of papers (see, for example, [7–12]). On the other hand, one can suppose that point P2 is a satellite of point P0, for instance. Its motion in a zero approximation (when m1 ¼ 0) can be described by the corresponding exact solution of the two-body problem but it is perturbed by the effect of attraction to point P1. As a result parameters of the point P2 orbit change with time, and the problem is to describe their evolution on the basis of the perturbation theory. Such a model was first successfully exploited in the theory of lunar motion [13,14]. It is usually called a satellite version of the restricted three-body problem and it is used in dynamics of the solar and stellar systems, and artificial satellites [4,15–18]. In the present paper we consider the satellite version of the restricted three-body problem when two points P0 and P1 form a binary system, losing the mass due to the corpuscular and photon radiation. We assume that the points masses vary isotropically with different rates with the only restriction that their total mass reduces according to the joint Meshcherskii law (1). We use the relative coordinates system with origin at point P0, where motion of point P1 is described by functions (2). Motion of point P2, being a satellite of point P0, is determined in a zero approximation by an exact solution of the two-body problem with variable masses that was obtained in [5]. But this motion is disturbed by the gravity of point P1. We are interested in the long-term evolution of the point P2 trajectory, based on an averaging scheme (see, for example, [17]), in the framework of Hill's approximation [13], when a distance between points P0 and P1 is considered to be much greater than a distance between P0 and P2. The main aim of the paper is to find possible values of the system parameters for which the differential equations, determining the secular perturbations of the point P2 orbital elements, are integrable and describe its quasi-elliptic motion. The corresponding solutions of the evolutionary equations are obtained in the analytical form. The paper is organized as follows. In Section 2 we obtain the equations of motion of point P2 in terms of the Delaunay elements. Then in Section 3 we introduce the Hamiltonian of the system in Hill's approximation and doubly average it over the mean anomalies of the points P1 and P2. As a result we obtain the system of differential equations, determining long-term evolution of the orbital elements. In Section 4 we analyze the averaged equations of motion to find the domains of the system parameters for which the solution, describing a quasi-elliptic motion of point P2, can be obtained. Then in Section 5 we determine the mass variation laws for which a general solution of the averaged equations of motion can be found in analytical form. And we conclude in Section 6.

Denoting the relative Cartesian coordinates of points P1, P2 as ! ! R ¼ ðX; Y; ZÞ, r ¼ ðx; y; zÞ, respectively, we can write the equations of motion of point P2 in the form (see [5]) ! r γ€ ! € !  0 r ¼ grad!W; r þ fm00 ð4Þ r γ 0 ðtÞr3 γ 0 ðtÞ where the force function W is given by ! !! γ€ 0 2 1 r  R r þ fm1  ; W¼ 2γ 0 ðtÞ Δ12 R3 

Δ12 ¼ ðX  xÞ2 þ ðY yÞ2 þ ðZ  zÞ2

1=2

;

 1=2  1=2 ; R ¼ X2 þ Y 2 þ Z2 ¼ a1 vðtÞ; r ¼ x2 þ y2 þ z2 and expressions (2) are taken into account. Here γ 0 ðtÞ is an increasing twice continuously differentiable function given by

γ 0 ðtÞ ¼

m0 ðt 0 Þ ; m0 ðtÞ

ð5Þ

which determines the dependence of point P0 mass on time. Note ! that the left-hand side of Eq. (4) contains the term γ€ 0 r =γ 0 that is also obtained in the right-hand side after differentiation of the force function W. This term has been added to ensure integrability of Eq. (4) in case of the unperturbed problem (W ¼0) of two bodies of variable mass for arbitrary function γ 0 ðtÞ (see [5]). Remind that our aim is to find a class of functions γ 0 ðtÞ, for which the evolutionary equations are integrable. In case of W¼0 Eq. (4) has an exact solution that can be obtained by integrating the corresponding Hamilton–Jacobi equation, for example (for details see [5]). This solution describes an aperiodic motion of point P2 on quasi-conical section and can be represented in the form x ¼ γ 0 ρ cos θ cos φ;

y ¼ γ 0 ρ cos θ sin φ;

z ¼ γ 0 ρ sin θ;

ð6Þ

where ρ; θ; φ are the quasi-spherical coordinates. We use here a term ”quasi” to emphasize that the unperturbed trajectories determined by functions (6) are conic sections only in the case of constant masses, i.e., when γ 0 ðtÞ ¼ 1. The conic is described by the equation

ρ¼

að1 e2 Þ ; 1 þe cos ν

ð7Þ

and it is an ellipse for the eccentricity eo 1, a is the semi-major axis of the ellipse. The true anomaly ν as a function of time is determined by the following equation: pffiffiffiffiffiffi Z ν K0 dν ¼ ðΦðtÞ  ΦðτÞÞ; ð8Þ 2 3=2 ð1 þe cos ν Þ a ð1  e2 Þ3=2 0 where

ΦðtÞ ¼

Z

t 0

dt

γ 20 ðtÞ

;

K 0 ¼ fm00 ;

and τ is the time of perihelion passage. It should be noted that the presence of a time-dependent scaling coefficient γ 0 ðtÞ in (6) results in the deformation of the conical section and aperiodic motion of point P2. For that reason, the solution to the problem of two bodies of variable mass is said to describe aperiodic motion on quasiconical sections (see [5]). Similarly to the case of constant masses, we can define three additional parameters, namely, the orbital inclination i, the longitude of the ascending node Ω, and the argument of the perihelion ω, determining orientation of the orbital plane and position of the conic in this plane (see, for example, [2, Chapter 2.5]). Thus, we

60

A.N. Prokopenya et al. / International Journal of Non-Linear Mechanics 73 (2015) 58–63

obtain six arbitrary constants i, Ω, e, a, ω, τ, known from the classical two-body problem as Keplerian orbital elements. They are determined by the initial conditions of motion and together with (6)–(8) they completely determine the solution of equation (4) in case of W¼ 0. To determine the unperturbed coordinates of point P2 as explicit functions of time, we need to obtain the function νðtÞ from transcendental equation (8). To simplify this problem it is convenient to introduce the eccentric anomaly w instead of true anomaly ν, using the following equation: rffiffiffiffiffiffiffiffiffiffi ν 1þe w tan : ð9Þ tan ¼ 1e 2 2

motion can be written in the standard Hamiltonian form as

Then the integral in the left-hand side of (8) can be evaluated, and the formula takes the form pffiffiffiffiffiffi K0 w  e sin w ¼ l  3=2 ðΦðtÞ  ΦðτÞÞ; ð10Þ a

Let us assume that the condition r 5 R is fulfilled and the Hill approximation may be used (see [13]). It means that the ratio r/R may be considered as a small parameter. Therefore, we may expand the force function W in (13) into a power series in terms of r/R and neglect all terms of order ðr=RÞ3 and higher. Then the Hamiltonian (13) takes the form

where parameter l may be called the mean anomaly, similar to the case of constant masses. If the function γ 0 ðtÞ is given, then the expression for l(t) in the right-hand side of (10) can be calculated, and thus the implicit function w ¼ wðtÞ can be obtained. Using (7) and (9), we can represent the unperturbed Cartesian coordinates (6) as x ¼ γ 0 aðð cos w  eÞð cos ω cos Ω  sin ω sin Ω cos iÞ pffiffiffiffiffiffiffiffiffiffiffiffi  1  e2 ð sin ω cos Ω sin w þ cos ω sin Ω sin w cos iÞÞ; y ¼ γ 0 aðð cos w  eÞð cos ω sin Ω þ sin ω cos Ω cos iÞ pffiffiffiffiffiffiffiffiffiffiffiffi  1  e2 ð sin ω sin Ω sin w  cos ω cos Ω sin w cos iÞÞ; z ¼ γ 0 aðð cos w  eÞ sin ω sin i þ

pffiffiffiffiffiffiffiffiffiffiffiffi 1  e2 cos ω sin w sin iÞ:

ð11Þ

Note that six constants i, Ω, e, a, ω, ΦðτÞ, arising during integration of the unperturbed equation (4), coincide with the classical Keplerian orbital elements in case of m0 ¼ const or γ 0 ¼ 1. So it is reasonable to consider them as the analogues of orbital elements (see [5]). Interaction with massive point P1 results in disturbing the point P2 trajectory determined by (11). It may be treated in such a way that solution of equation (4) has form (11) but the orbital parameters become functions of time. To derive the differential equations determining the time evolution of orbital elements in the simplest form, it is convenient to rewrite Eq. (4) in the Hamiltonian form and to change to the special set of canonical variables known as Delaunay's elements (see [5,19]). Three pairs of the corresponding canonical conjugate coordinates and momenta ðl; LÞ, ðg; GÞ, ðh; HÞ are related to the analogues of the Keplerian orbital elements as pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi K0 l ¼ 3=2 ðΦðtÞ  ΦðτÞÞ; L ¼ K 0 a; a g ¼ ω;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G ¼ K 0 að1 e2 Þ;

ð12Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h ¼ Ω; H ¼ K 0 að1  e2 Þ cos i:

H¼ 

2γ 20 L2

 W;

dL ∂H ¼ ; dt ∂l

dg ∂H ¼ ; dt ∂G

dG ∂H ¼ ; dt ∂g

dh ∂H ¼ ; dt ∂H

dH ∂H ¼ : dt ∂h

ð14Þ

3. Averaged equations in Hill's approximation

γ€ þ 0 ðx2 þ y2 þ z2 Þ 2γ 20 L2 2γ 0  K1  2 2 þ x þ y2 þ z2  3ðx2 cos 2 θ1 þ y2 sin θ1 þ xy sin ð2θ1 ÞÞ ; 2γ 1 v3 a31

H¼ 

K 20

ð15Þ where K 1 ¼ fm10 is a constant, γ 1 ¼ m1 ðt 0 Þ=m1 ðtÞ is a doubly differentiable continuous function of time, determining a variation of the point P1 mass, and coordinates x; y; z are determined in (11). Taking into account Eq. (12), we can rewrite the Hamiltonian (15) in terms of the Delaunay elements but the corresponding expression turns out to be quite cumbersome and we do not write it here. It should be noted that in the Hill approximation the ratio a=a1 is a small value. Therefore, the third term in the Hamiltonian (15), being proportional to K1, contains a small coefficient ða=a1 Þ2 . We assume that the second derivative γ€ 0 is a small value, as well, what means that the masses m0 ðtÞ and m1 ðtÞ vary very slowly and they do not change substantially during one revolution of either point P1 or P2. Then Eqs. (3) and (14) with the Hamiltonian (15) show that the variables L, g, G, h, H will change slowly with time, undergoing small short-term perturbations, while the mean anomalies l and θ1 are increasing functions. As we are interested in the long-term evolution of the point P2 orbit under an influence of massive point P1, one may disregard the short-term perturbations of orbital elements by means of averaging of Eq. (14) over the mean anomalies of points P1 and P2 (see [17]). Instead of averaging Eq. (14), it is expedient to average the Hamiltonian (15) and to substitute it into (14) afterwards. The averaged Hamiltonian is determined as Z 2π Z 2π Z 2π Z 2π 1 1 H¼ 2 H dl dθ1 ¼ 2 Hð1  e cos wÞ dw dθ1 ; 4π 0 4π 0 0 0 and is given by H¼

K 20

þ 2γ 20 L2

γ 0 γ€ 0 a2 2



   K 1 γ 20 a2 3 3 2 e 1 þ e2 þ 1 þ 2 2 2γ 1 a31 v3

 3K 1 γ 20 a2  2  2 þ 2 cos 2 i þ e2 ð3 þ 3 cos 2 iþ 5 cos ð2ωÞ sin iÞ : 3 3 16γ 1 a1 v

The Hamiltonian in the Delaunay elements is represented in the form K 20

dl ∂H ¼ ; dt ∂L

ð13Þ

where the function W is defined in (4) and components of vectors ! ! R and r are given by (2) and (11), respectively. Then equations of

ð16Þ Now let us assume that the function differential equation:

γ€ 0 ¼  α where

K 1 γ 0 ðtÞ ; a31 γ 1 ðtÞv3 ðtÞ

γ0 satisfies the following ð17Þ

α is a parameter. In this case all terms in the Hamiltonian

A.N. Prokopenya et al. / International Journal of Non-Linear Mechanics 73 (2015) 58–63

61

(16), except for the first one, will have the same coefficient depending on time that can be eliminated by means of introducing a new independent variable, instead of time t. As a result, the averaged equations of motion will become integrable and their solution can be found in analytical form. Taking into account (12) and (17), we can rewrite the averaged Hamiltonian (16) in the form " ! ! K 1 γ 20 L4 1  α K 20 3G2 3 H2 1 þ 5  H¼ þ  4 2 2γ 20 L2 2γ 1 a31 v3 K 20 L2 G2 ! !!# 3 G2 3H 2 H2 1 2 3 þ 2 þ 5 cos ð2ωÞ 1  2 : ð18Þ  8 L G G

Using n as a new independent variable, we can rewrite the evolutionary differential equations (20)–(23) in the form

Obviously, the averaged Hamiltonian (18) does not depend on the mean anomaly l. According to (14), it means that L¼ const and the semi-major axis a ¼ L2 =K 0 is constant, as well. The first differential equation of system (14), determining evolution of the mean anomaly l, looks as pffiffiffiffiffiffi K 1 γ 20 a2 K0 _l ¼ pffiffiffiffiffiffiffiffiffi ð  4αð7 þ 3e2 Þ þ 7 15 cos ð2ωÞ þ γ 20 a3=2 8γ 1 a31 v3 K 0 a

 dΩ 2 cos i  ¼ pffiffiffiffiffiffiffiffiffiffiffiffi 5e2 cos ð2ωÞ  2  3e2 : 2 dn 1e

þ 3 cos 2 ið  7 þ5 cos ð2ωÞÞ þ 3e2 ð1  3 cos 2 i 5 cos ð2ωÞ sin 2 iÞÞ:

pffiffiffiffiffiffiffiffiffiffiffiffi de ¼ 10e 1  e2 sin 2 i sin ð2ωÞ; dn

ð25Þ

di 10e2 ¼  pffiffiffiffiffiffiffiffiffiffiffiffi sin i cos i sin ð2ωÞ; dn 1  e2

ð26Þ

dω 2  ¼ pffiffiffiffiffiffiffiffiffiffiffiffi 5 cos 2 i  1 þ e2 þ 4αð1  e2 Þ dn 1  e2  2 2  5e cos ð2ωÞ þ 5 sin i cos ð2ωÞ ;

ð27Þ ð28Þ

Note that the system of differential equations (25)–(28) looks similarly to the corresponding equations, describing evolution of satellites of Uranus (see [20]). But due to dependence of the points masses on time equation (27) contains additional term 4αð1  e2 Þ in the parentheses in the right-hand side and additional parameter α. Therefore, the system dynamics should be more complicated, although it can be investigated in a similar way as in [20].

ð19Þ As the mean anomaly l determines position of the point P2 on its orbit and does not influence on the orbit itself, further we will not consider Eq. (19), it is sufficient to note that it can be integrated as soon as other orbital parameters have been found. On substituting the averaged Hamiltonian (18) into the last four equations of (14) and taking into account (12), after quite standard calculations we obtain the following system of the evolutionary differential equations: pffiffiffiffiffiffiffiffiffiffiffiffi 15K 1 γ 20 a2 pffiffiffiffiffiffiffiffiffi e 1  e2 sin 2 i sin ð2ωÞ; e_ ¼ 3 3 K0a 8γ 1 a1 v i_ ¼ 

15K 1 γ 20 a2 e2 pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi sin i cos i sin ð2ωÞ; 3 3 K 0 a 1  e2 8γ 1 a1 v

3K 1 γ 20 a2 1  pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 5 cos 2 i 1 þ e2 þ4αð1  e2 Þ 8γ 1 a31 v3 K 0 a 1  e2  2 5e2 cos ð2ωÞ þ 5 sin i cos ð2ωÞ ;

ð20Þ

 3K 1 γ 20 a2 cos i  pffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 5e2 cos ð2ωÞ  2  3e2 : 8γ 1 a31 v3 K 0 a 1  e2

dn ¼

16γ

γ

dt:

ð1  e2 Þ cos 2 i ¼ c1 ¼ const;  e2

ð29Þ

 2 N  sin 2 i sin 2 ωÞ ¼ c2 ¼ const; 5

ð30Þ

ð21Þ

ð22Þ

sin 2 i ¼

ð23Þ

where z ¼ e2 . Taking them into account, we can eliminate the variables i and ω in Eq. (25) and obtain a differential equation with respect to z(n) that has a form

Now we can introduce a new dimensionless variable n according to the relationship 3K 1 20 ðtÞa2 pffiffiffiffiffiffiffiffiffi 3 3 1 ðtÞa1 v ðtÞ K 0 a

One can readily check that the system of three equations (25)– (27) has two independent integrals of motion

where new parameter N ¼ 1 þ α has been introduced for convenience. As we are interested in quasi-elliptic motion of point P2, eccentricity of its orbit should be less than 1 and, hence, the first integral c1 must belong to the interval 0 r c1 o 1. Consequently, for given c1, we have 0 r e2 o 1  c1 . Using the integrals (29) and (30), we obtain the following expressions:

ω_ ¼

_ ¼ Ω

4. Analysis of evolutionary equations

ð24Þ

1  z c1 ; 1z

sin 2 ω ¼

ð1  zÞð2Nz  5c2 Þ ; 5zð1  z  c1 Þ

pffiffiffiffiffiffiffiffiffiffi dz ¼ 8 sgnð sin ð2ω0 ÞÞ Q ðzÞ; dn

ð31Þ

where the function sgnðxÞ determines the sign of sin ð2ω0 Þ at the initial instant of time (ω0 ¼ ωðt 0 Þ), and the third degree polynomial

Fig. 1. Domains of possible values of the integrals c1, c2 for N o 0 (left) and N Z 0 (right) are bounded by bold curves.

62

A.N. Prokopenya et al. / International Journal of Non-Linear Mechanics 73 (2015) 58–63

Q(z) is given by 2

Q ðzÞ ¼ ð2Nz  5c2 Þð5c2 þ zð5  2N  5c1  5c2 Þ  z ð5  2NÞÞ;

ð32Þ

Analysis of Eq. (31) shows that under the conditions 0 r c1 r 1;

0 r z r 1  c1 ;

ð38Þ

0 r sin 2 ω r 1;

where

a quasi-elliptic motion takes place when either z¼ const and Q ðzÞ ¼ 0 or z varies between two different positive zeros zi, zk of polynomial Q(z), provided that the polynomial is non-negative inside the interval z A ½zi ; zk . The corresponding values of integrals c1, c2 must belong to the domains shown in Fig. 1, which are bounded by the lines c1 ¼ 0;

c1 ¼ 1 

5c2 ; 2N

ð33Þ

and the curve

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 2N 2N  2 ð c2 Þ 1  : c1 ¼ 1 c2  5 5

ð34Þ

In contrast to the case of constant masses (see [20]), a shape and size of the domain of possible values of integrals c1, c2, corresponding to quasi-elliptic motion, depends on parameter N. Eq. (30) shows that c2 r0 for N o 0 and c2 Z 0 for N 4 5=2. In both cases c2 ¼ 0 only if z ¼0 (or e¼ 0) and it will be the only solution of equation (31). In case of 0 r N r 5=2 the polynomial Q(z) has three roots, namely, z1 ¼ 1  c1 =ð1  2N=5Þ and z2;3 ¼ 0. The root z1 is negative for 1  2N=5 o c1 r1 (the corresponding points are shown in the right-hand side of Fig. 1 as a dashed bold line), and Eq. (31) has the only solution z ¼0 again. But for 0 r c1 o 1  2N=5 the root z1 becomes smaller than 1 and the polynomial (32), taking a form Q ðzÞ ¼ 2Nz2 ð5  2N  5c1  zð5  2NÞÞ;

ð35Þ

is non-negative for 0 r z rz1 . Substituting the polynomial (35) into Eq. (31), one can readily see that the differential equation is integrated in elementary functions and its solution is determined by the following equation: pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi b  b  kz p ln ffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ¼ 8 10N sgnð sin ð2ω0 ÞÞn þ B0 ; ð36Þ b þ b  kz where

pffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b  b  kz0 B0 ¼ lnpffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; b þ b  kz0 z0 is the initial value of z at n ¼0 (t ¼ t 0 ), and b ¼ 1 c1 2N=5, k ¼ 1 2N=5. The case of line c1 ¼ 1 5c2 =ð2NÞ (see Fig. 1) is analyzed in a similar way. The polynomial Q(z) takes the form     5 1 ; ð37Þ Q ðzÞ ¼ ð2Nz  5c2 Þ2 z 1  2N and has three roots z1 ¼

1 5 1 2N

;

z2;3 ¼

corresponding solution is given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi z2;3  z1 þ z  z1 2 z2;3 p ffiffiffiffiffiffiffiffiffiffiffi ffi lnpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 16 N  1 sgnð sin ð2ω0 ÞÞn þ B0 ; z1 z2;3  z1  z  z1

5c2 : 2N

For N Z 5=2 we have 0 r z2;3 r 1 but z1 41 and so Eq. (31) has only a stationary solution z ¼ 5c2 =ð2NÞ ¼ 1  c1 , describing quasi-elliptic motion of point P2. In case of 0 r N o 5=2 and c2 4 0 the root z1 becomes negative and Eq. (31) has only a stationary solution z ¼ 1  c1 again. But for N o 0, when the line (33) touches the curve (34) (see Fig. 1) and parameters c1, c2 must satisfy the conditions 0 r c1 r 1=ð1  2N=5Þ, 2N=5 rc2 r2Nð1  c1 Þ=5, we obtain 0 oz1 r z2;3 r 1 and polynomial (37) is non-negative for z A ½z1 ; z2;3 . Then Eq. (31) is integrated in elementary functions and the

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi z2;3  z1 þ z0 z1 B0 ¼ lnpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi; z2;3  z1  z0 z1 and z0 is the initial value of z. On the curve (34) the polynomial Q(z) takes the form sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!2   5c2 ð  c2 Þ z Q ðzÞ ¼ 2Nð2N 5Þ z  ; 1  2N=5 2N

ð39Þ

where we have taken into account that c2 o 0 and N o 5=2. The corresponding roots are given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð  c2 Þ 5c2 ; z3 ¼ : z1;2 ¼ 1 2N=5 2N One can readily check that for 0 rN o 5=2 we have 0 r z1;2 r 1 and z3 o 0 and, hence, Eq. (31) has only a stationary solution z ¼ z1;2 . This solution remains also for N o0 and  1 þ 2N=5 r c2 o 2N=5 when the root z3 becomes greater than 1. But for N o 0 and 2N=5 rc2 ocn2 , where cn2 corresponds to the point of contact of the curve (34) and the line (33), we obtain 1 o z1;2 o z3 r 1; 1  5=ð2NÞ and polynomial Q(z) is non-negative for z A ½z1;2 ; z3 . For these values of parameters equation (31) is integrated in elementary functions as in two previous cases. The last boundary of the domain of possible values of c1, c2 is the line c1 ¼ 0, corresponding to the case i ¼ π =2. The polynomial Q(z) takes the form Q ðzÞ ¼ ð2Nz 5c2 Þð1 zÞð5c2 þ zð5  2NÞÞ;

ð40Þ

and has three different roots z1 ¼ 1;

z2 ¼

5c2 ; 2N  5

z3 ¼

5c2 : 2N

ð41Þ

One can readily see that for N ¼ 5ð1 þ c2 Þ=2 we have a double root z1 ¼ z2 ¼ 1 and 0 r z3 o1. In this case Eq. (31) is integrated in elementary functions as in previous cases. When all roots are different, depending on the value of parameter N, one can distinguish several cases when two of the roots (41) belong to the interval 0 r zj r zk r1 and Q ðzÞ Z 0 for zj r z r zk . They are (a) (b) (c) (d) (e)

N 4 5ð1 þ c2 Þ=2, 0 rc2 r2N=5 : z3 rz r z2 ; 5=2 oN o 5ð1 þ c2 Þ=2, 0 rc2 r2N=5: z3 r z r z1 ; 0 rN o 5=2, 0 rc2 r 2N=5 : z3 rz r z1 ; 0 rN o 5=2,  1 þ 2N=5 r c2 o 0 : z2 r z r z1 ; N o0,  1 þ 2N=5 r c2 o0 : z2 r z rz1 ¼ 1.

All such cases are analyzed in a similar way and Eq. (31) is integrated in the elliptic quadratures. As an example, let us consider the case (d). We have z3 o 0 and Q ðzÞ ¼ ðz1  zÞðz  z2 Þðz  z3 Þ Z0 for z A ½z2 ; z1 . Then Eq. (31) is integrated in quadrature and its solution is given by Z z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dz pffiffiffiffiffiffiffiffiffiffi: ð42Þ 8 2Nð5  2NÞ sgnð sin ð2ω0 ÞÞn ¼ Q ðzÞ z0 An integral in the right-hand side of (42) is calculated in terms of the elliptic functions and the result may be represented as (for

A.N. Prokopenya et al. / International Journal of Non-Linear Mechanics 73 (2015) 58–63

Making a substitution γ 0 ðtÞ ¼ vðtÞξðτÞ in differential equation (24), determining the variable n as function of time, we obtain

details see [20]) zðuÞ ¼

z2  z3 κ sn u ; 1  κ 2 sn2 u 2

2

ð43Þ

where

κ2 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ¼ 8 2Nð5  2NÞn þ u0 ;

u0 ¼ sgnð sin ð2ω0 ÞÞFðφ0 ; κ 2 Þ;

sin 2 φ0 ¼

ðz1  z3 Þðz0  z2 Þ : ðz1  z2 Þðz0  z3 Þ

5. Possible laws of mass variation Remind that we have obtained the evolutionary equations (20)–(23) under assumption that function γ 0 ðtÞ satisfies the differential equation (17). Taking into account Eq. (1) and the relationships m0 ðtÞ ¼ m00 =γ 0 ðtÞ, m1 ðtÞ ¼ m10 =γ 1 ðtÞ, we can rewrite Eq. (17) in the form Kγ K γ€ 0 ¼  α 3 04 þ α 3 03 : a1 v a1 v

ð44Þ

Making a substitution γ 0 ðtÞ ¼ vðtÞξðτÞ, we can reduce (44) to the form ! 2 d ξ K K 2 þ AC  B þ α ξðτÞ ¼ α 30 ; ð45Þ dτ 2 a1 a31 where τ is a new independent variable determined by the differential equation dt 2

At þ 2Bt þ C

:

ð46Þ

One can readily see that Eq. (45) is integrable and its general solution is given by

ξðτÞ ¼ α

K0

σ 2 a31

þ Φ1 cos ðστÞ þ Φ2 sin ðστÞ;

ð47Þ

where Φ1, Φ2 are two arbitrary constants, determined by the initial conditions, and K a1

σ 2 ¼ AC  B2 þ α 3 : General solution of equation (46) is also easily obtained and is 1

At þ B

τðtÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi arctanpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ τ0 ; AC  B2

AC  B2

ð49Þ

0

where the function ξðτÞ is given by (47). After integration of this equation, what can be easily done, we obtain an explicit expression for the function nðτÞ. Thus, the problem of integration of the evolutionary equations in the restricted three-body problem with variable mass has been completely solved.

Here snu and Fðφ0 ; κ 2 Þ are the Jacobi elliptic sine and the incomplete elliptic integral of the first kind, respectively. Note that in the internal points of the domains shown in Fig. 1 we have similar situation, when three roots of the polynomial Q(z) are different and there is an interval ½zj ; zk   ½0; 1, where Q ðzÞ Z 0. In all such cases Eq. (31) is integrated in a similar way and the results are expressed in terms of the elliptic functions. As soon as the function z(n) is found, one can write the solutions of differential equations (26)–(28) and investigate an influence of the system parameters variation on the quasi-elliptic orbital elements of point P2. But to find the orbital elements as explicit functions of time we still need to know the laws of mass variation γ 0 ðtÞ and γ 1 ðtÞ.

dτ ¼

dn 3a3=2 2 ¼ ðK ξ ðτÞ K 0 ξðτÞÞ; dτ 16a3 K 1=2 1

z1 z2 ; z1 z3

63

ð48Þ

where the constant τ0 is determined from the condition τðt 0 Þ ¼ 0. Thus, Eqs. (1), (47) and (48) completely determine a general solution of equation (44) and possible mass variation laws of the points P0, P1. One should remember only that parameters A, B, C should be chosen in such a way that the functions v(t), γ 0 ðtÞ and γ 1 ðtÞ must decrease with time.

6. Conclusion We have considered the satellite version of the restricted threebody problem in the case when masses of two bodies vary isotropically with different rates, while their total mass reduces according to the joint Meshcherskii law. We have obtained the evolutionary equations of the massless point P2, describing the secular perturbations of the orbital elements, in the Hill approximation, and investigated their integrability. We have shown that these equations are integrated in terms of the elementary and elliptic functions and their solutions describe quasi-elliptic motion of the point P2 if initial conditions of motion are chosen in such a way that two integrals of motion c1, c2 belong to the domains shown in Fig. 1. We have also found possible mass variation laws admitting integration of the evolutionary equations in analytic form. Using the results obtained, one can visualize phase trajectories of the point P2 and investigate their dependence on system parameters and different mass variation laws, this work will be done in our next paper. References [1] V.G. Szebehely, Theory of Orbits. The Restricted Problem of Three Bodies, Academic Press, New York, London, 1967. [2] D. Boccaletti, G. Pucacco, Theory of Orbits. vol. 1: Integrable Systems and Nonperturbative Methods, Springer-Verlag, Berlin, Heidelberg, 1996. [3] A.P. Markeev, The Libration Points in Celestial Mechanics and Cosmodynamics, Nauka, Moscow, 1978. [4] A.A. Bekov, T.B. Omarov, The theory of orbits in non-stationary stellar systems, Astron. Astrophys. Trans. 22 (2003) 145–153. [5] M.Zh. Minglibayev, Dinamika gravitiruyushchikh tel s peremennymi massami i razmerami (Dynamics of Gravitating Bodies of Variable Masses and Sizes), LAP Lambert Academic Publications, Library of the Al.-Farabi Kazakh National University, Kazahstan, 2012. [6] L.M. Berkovič, Gylden–Meščerski problem, Celest. Mech. 24 (1981) 407–429. [7] J. Singh, B. Ishwar, Effect of perturbations on the location of equilibrium points in the restricted problem of three bodies with variable masses, Celest. Mech. 32 (1984) 297–305. [8] L.G. Luk'yanov, Particular solutions in the restricted three-body problem with arbitrary varying masses, Sov. Astron. 33 (2) (1989) 194–197. [9] L.G. Luk'yanov, Stability of libration points in the restricted three-body problem with varying masses, Sov. Astron. 34 (1) (1990) 87–89. [10] A.A. Bekov, Particular solutions in the restricted, collinear three-body problem with variable masses, Sov. Astron. 35 (1) (1991) 103–106. [11] J. Singh, Nonlinear stability of equilibrium points in the restricted three-body problem with variable mass, Astrophys. Space Sci. 314 (2008) 281–289. [12] P.S. Letelier, T.A. da Silva, Solutions to the restricted three-body problem with variable masses, Astrophys. Space Sci. 332 (2011) 325–329. [13] G.W. Hill, Researches in the lunar theory, Am. J. Math. 1 (1878) 129–147. [14] E.W. Brown, C.A. Shook, Planetary Theory, Dover Publications, New York, 1964. [15] M.L. Lidov, Evolution of the orbits of artificial satellites of planets as affected by gravitational perturbations from external bodies, AIAA J. Russian Suppl. (1963) 1985–2002. [16] Y. Kozai, The motion of a lunar orbiter, J. Astron. Soc. Jpn. 15 (1963) 301–312. [17] M.L. Lidov, M.A. Vashkov'yak, On quasi-satellite orbits in a restricted elliptic three-body problem, Astron. Lett. 20 (5) (1994) 676–690. [18] A.A. Bekov, On the dynamics of non-stationary binary stellar systems, in: Order and Chaos in Stellar and Planetary Systems, ASP Conf. Series, vol. 316, 2004, pp. 366–370. [19] A.N. Prokopenya, M.Zh. Minglibayev, G.M. Mayemerova, Symbolic calculations in studying the problem of three bodies with variable masses, Progr. Comput. Softw. 40 (2) (2014) 79–85. [20] M.A. Vashkov'yak, Evolution of orbits of distant satellites of Uranus, Astron. Lett. 25 (7) (1999) 476–481.