Modelling the Evolution of the Two-Planetary Three-Body System of Variable Masses

. A classical non-stationary three-body problem with two bodies of variable mass moving around the third body on quasi-periodic orbits is considered. In addition to the Newtonian gravitational attraction, the bodies are acted on by the reactive forces arising due to anisotropic variation of the masses. We show that New-tonian’s formalism may be generalized to the case of variable masses and equations of motion are derived in terms of the osculating elements of aperiodic motion on quasi-conic sections. As equations of motion are not integrable the perturbative method is applied with the perturbing forces expanded into power series in terms of eccentricities and inclinations which are assumed to be small. Averaging these equations over the mean longitudes of the bodies in the absence of a mean-motion resonances, we obtain the differential equations describing the evolution of orbital parameters over long period of time. We solve the evolution equations numerically and demonstrate that the mass change modify essentially the system evolution.


Introduction
The classical three-body problem [26] is a well-known model of celestial mechanics in the framework of which a motion of three bodies P 0 , P 1 , P 2 of masses m 0 , m 1 , m 2 , respectively, under their mutual gravitational attraction is studied.The bodies are assumed to interact according to Newton's law of gravitation that is a good approximation in the case of spherically symmetric bodies.However, real heaven bodies may have different shape and mass distribution and their interaction is more complicated.To make model more realistic some authors take into account an oblateness of heaven bodies and some other perturbing forces like radiation pressure and quantum effects (see, for example, [1,2,13]).Such perturbations can modify motion of the bodies or change the regions of permissible motion but the problem is the stationary one as physical parameters of the system remain constant.
From the other side, real-life celestial bodies are non-stationary; their characteristics such as mass, size, shape, and internal structure, may vary with time (see, for example, [8,9,25]).Non-stationarity of the bodies may influence the dynamical evolution of their systems that makes a study of such systems highly relevant (see [3,11,16]).At the same time, non-stationarity complicates essentially the corresponding mathematical models of the bodies motion.Even in the case of classical two-body problem, a general solution of which is wellknown, dependence of mass on time makes the problem non-integrable; only in some special cases its exact analytical solution can be found (see survey of such models in [4,23,24]).
The bodies masses influence essentially on their interaction and motion and so it is especially interesting to investigate the dynamics of the many-body system with variable mass.One of the first works in this direction were done by T.B. Omarov [20] and J.D. Hadjidemetriou [10] who started investigation of the non-stationary two-body problem and showed that mass variability affects essentially on the dynamic evolution of the system.Later these investigations were generalized to the system of three bodies although works in this field are not numerous (see, for instance, [3,14,27,28]).
As the equations of motion are non-integrable the perturbation theory is usually used (see [5]).Its application involves quite cumbersome symbolic computation which can be best performed with computer algebra [21].Investigations of the three-body problem with variable masses, changing isotropically or anisotropically, were continued in a series of works [15,17,18,22], where equations of motion were obtained in terms of the second system of Poincaré elements (see [7]) in the framework of the Hamiltonian formalism.In order to obtain equations of motion accurate to linear terms of the orbital elements we need to compute the series expansion of the perturbing functions in terms of the orbital elements up to second order inclusive.It should be noted that it is a nontrivial task and we demonstrated that it may solved efficiently with the aid of the computer algebra (see [15,18,21,22]).
In the present work, we study the dynamical evolution of two-planetary system of three bodies when two planets P 1 , P 2 move around a central star P 0 in quasi-elliptic orbits such that their orbits do not intersect.In the first approx-imation, these orbits are determined by the exact solutions of the unperturbed equations of motion which can be obtained in analytic form for arbitrary laws of mass variation of the bodies (see [16]).Mutual attraction of the bodies P 1 , P 2 and reactive forces arising in the case of anisotropic mass variation enforce the orbital elements to change.In contrast to our previous works (see [15,18,22]), the differential equations determining the perturbed motion of the bodies are obtained in terms of the osculating elements in the framework of Newton's formalism what enables to write out expressions for the reactive forces and to obtain directly differential equations for the orbital elements (see [6]).In the case of small eccentricities and inclinations of the orbits the perturbing forces may be expanded in series in these parameters up to any desired order but here we consider only the first order terms what is sufficient to obtain the results corresponding to the accuracy of the observations.Averaging the equations of the perturbed motion over mean longitudes of the bodies P 1 , P 2 in the absence of a mean-motion resonances, we obtain the differential equations describing the evolution of orbital elements over long periods of time.These equations are solved numerically for different laws of the masses change.All relevant symbolic and numerical calculations are performed here with the aid of the computer algebra system Wolfram Mathematica [29].

Model description
Consider a system of three bodies of variable mass attracting each other according to Newton's law of universal gravitation.Denoting the position vectors of the bodies P 1 , P 2 relative to the primary P 0 by ⃗ r j = (x j , y j , z j ) and applying Newton's second law, the equations of motion may be written as (see [15,16,18]) Here G is the constant of gravitation, and the twice differentiable functions γ 1 (t) and γ 2 (t) are defined by where m 00 = m 0 (t 0 ), m j0 = m j (t 0 ) are the masses of the bodies P 0 , P 1 , P 2 , respectively, at the initial instant of time.The forces ⃗ F 1 and ⃗ F 2 on the righthand side of (2.1) can be represented by where and the reactive forces ⃗ Q 1 , ⃗ Q 2 are determined by the expressions (see [12]) The dot above a symbol in (2.2)-(2.4)denotes the total time derivative of the corresponding function, and ⃗ V j , (j = 0, 1, 2) are the relative velocities of the particles leaving the body P j or falling on it.
Note that in the case of constant masses when γ 1 (t) = 1, γ 2 (t) = 1 equations (2.1) reduce to the well-known equations determining relative motion of the bodies in the classical three-body problem.These equations are not integrable and are usually studied by methods of perturbation theory using an exact solution of the two-body problem as the first approximation (see, for example, [5,19]).Similar approach may be also applied in the case of variable masses but the corresponding two-body problem is integrable only for some special laws of mass change (see [23,24]).To obtain integrable two-body problem for arbitrary law of the mass variation we add the terms γj /γ j ⃗ r j in the left-hand side of Equations (2.1) and in expressions (2.2), (2.3) for the forces ⃗ F 1 , ⃗ F 2 .As a result, a general solution of equations obtained from (2.1) at ⃗ F 1 = 0, ⃗ F 2 = 0 could be written for arbitrary laws of mass variation of the bodies.
Actually, at ⃗ F j = 0, (j = 1, 2) two equations (2.1) become independent of each other and each of them has an exact solution that describes aperiodic motion of the body P j , (j = 1, 2) on a quasi-conic section (see [16]); it can be written as x j = γ j ρ j (cos(ω j + ν j ) cos Ω j − sin(ω j + ν j ) sin Ω j cos i j ) , y j = γ j ρ j (cos(ω j + ν j ) sin Ω j + sin(ω j + ν j ) cos Ω j cos i j ) , z j = γ j ρ j (sin(ω j + ν j ) sin i j ) , where ν j is the true anomaly and ρ j = a j (1 − e 2 j )/(1 + e j cos ν j ). (2.6) The constants a j , e j , i j , Ω j and ω j in (2.5), (2.6) are analogs of the well-known Kepler orbital elements and are determined from the initial conditions of motion (see [16]).The true anomaly ν j characterizes the position of the body on the orbit; introducing an analog of the eccentric anomaly E j by the relation we obtain the known Kepler equation where the mean anomaly M j is given by and κ j = G(m 00 + m j0 ), (j = 1, 2).The functions Φ j (t) have the form (2.10) By τ j in (2.9) we denote an analog of the time when the body P j passes through the pericenter.
It is readily seen that, for given orbital elements a j , e j , i j , Ω j , ω j , and τ j of each of the bodies P 1 and P 2 and the known functions γ 1 (t) and γ 2 (t), which depend on the laws of mass variation of all three bodies, Equations (2.7)-(2.10)make it possible to find the mean anomalies M j , the eccentric anomalies E j and true anomalies ν j as functions of time.As a result, solutions (2.5), (2.6) enable to compute the relative Cartesian coordinates of the bodies P 1 and P 2 at ⃗ F 1 = 0, ⃗ F 2 = 0 as functions of time and to describe their unperturbed motion.Using Equations (2.6)-(2.10),one can write the total time derivatives of the coordinates (2.5) in the form where (1 + e j cos ν j ) 2 γ 2 j (t) .

Equations of perturbed motion
In the absence of forces (2.2), (2.3) the orbital elements a j , e j , i j , Ω j , ω j , and τ j of the bodies P 1 , P 2 do not change with time.However, mutual attraction and reactive forces (2.4) arising in the case of anisotropic mass variation of the bodies affect their motion and the orbital elements must necessarily vary with the time.Solving Equations (2.1) at ⃗ F 1 ̸ = 0, ⃗ F 2 ̸ = 0 numerically, one can find the perturbed coordinates of the bodies as functions of time but it will be equally effective to obtain the orbital elements as functions of the time.These functions may be used then to investigate the long-term evolution of orbital elements which is the most interesting for applications in celestial mechanics.
Taking into account the dependence of the orbital elements on time, the solutions to Equations (2.1) can be written in the general form x j = x j (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), y j = y j (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), z j = z j (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), where the functions in the right-hand sides are determined by expressions (2.5), in which true anomalies ν j are replaced by the mean anomalies M j and relations (2.9), (2.10) determine the explicit dependence of the mean anomalies M j on the parameters a j and τ j .Such representation of solutions is well-known in the theory of differential equations as the method of the variation of arbitrary constants.
Direct substitution of solutions (3.1) into Equations (2.1) gives six secondorder differential equations for 12 unknown functions a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (t), j = 1, 2. Since each such system has an infinite number of solutions, we should introduce six additional equations for the variables a j , e j , i j , Ω j , ω j , M j , j = 1, 2. Usually such equations are obtained from the condition that the rates of variation of the perturbed coordinates ẋj , ẏj , żj are equal to the partial derivatives of functions (3.1) with respect to time.It means that in perturbed motion both the coordinates and the velocity components at time t are given by the formulas (2.5), (2.11) expressed in terms of the time and the instantaneous orbital elements at t.Such instantaneous elements are known as osculating elements (for details, see [6]).
As a result, we obtain the following equations where j = 1, 2.
The time derivatives of the coordinates (3.1) can be formally written as ẋj = ẋj (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), ẏj = ẏj (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), żj = żj (a j (t), e j (t), i j (t), Ω j (t), ω j (t), M j (a j , τ j , t), t), where the functions in the right-hand sides are determined by expressions (2.11).Taking into account that solutions (2.5), (2.11) satisfy the equations of motion (2.1) in the absence of perturbations and substituting (3.3) into (2.1),we obtain the following equations for the functions ẋj , ẏj , żj where j = 1, 2. By solving Equations (3.2), (3.4), we can obtain explicit expressions for the derivatives of the orbital elements ȧj , ėj , ij , Ωj , ωj , and Ṁj for each body P 1 and P 2 .Note that carrying out the corresponding calculations and deriving the differential equations for the orbital elements requires quite cumbersome symbolic computations.Such computations can be performed efficiently using a computer algebra system such as Wolfram Mathematica (see [29]).By performing the corresponding calculations, we finally obtain the following system of differential equations for finding the dependence of the orbital elements on time: ) F nj (cos E j − e j ) cos ω j − 1 − e 2 j sin ω j sin E j , F nj sin i j (cos E j − e j ) sin ω j + 1 − e 2 j cos ω j sin E j , e j √ κ j (1 − e j cos E j ) 1 − e 2 j (−2 + e 2 j + e j cos E j ) sin E j F τ j The forces F rj , F τ j , and F nj on the right-hand sides of (3.5)-(3.6)are the radial, transversal and normal components of the forces ⃗ F 1 , ⃗ F 2 , respectively, determined by expressions (2.2), (2.3).As the reactive forces ⃗ Q 1 , ⃗ Q 2 defined by (2.4) are usually determined in the orbital systems of coordinates of the bodies P 1 , P 2 the forces ⃗ F 1 , ⃗ F 2 are also written in these systems of coordinates.The direction cosines of the unit vectors ⃗ e rj = (e xj , e yj , e zj ), ⃗ e τ j = (τ xj , τ yj , τ zj ), and ⃗ e nj = (n xj , n yj , n zj ) along the radial, transversal, and normal directions, respectively, can be easily written on the basis of solutions (2.5): e xj = cos(ω j + ν j ) cos Ω j − sin(ω j + ν j ) sin Ω j cos i j , e yj = cos(ω j + ν j ) sin Ω j + sin(ω j + ν j ) cos Ω j cos i j , e zj = sin(ω j + ν j ) sin i j , τ xj = − sin(ω j + ν j ) cos Ω j − cos(ω j + ν j ) sin Ω j cos i j , τ yj = − sin(ω j + ν j ) sin Ω j + cos(ω j + ν j ) cos Ω j cos i j , τ zj = cos(ω j + ν j ) sin i j , n xj = sin Ω j sin i j , n yj = − cos Ω j sin i J , n zj = cos i j , j = 1, 2. (3.8) Denoting the components of the relative velocities of particles leaving the bodies P 1 , and P 0 or falling on them along the radial, transversal, and normal directions in the orbital system of coordinates related to the body P 1 by V r1 , V r0 , V τ 1 , V τ 0 , V n1 , and V n0 and using (2.2)-(2.4),we obtain for the first body where the corresponding components of the reactive force ⃗ Q 1 are given by Similarly, denoting the components of the relative velocities of particles leaving body P 2 or falling on it along the radial, transversal, and normal directions in the orbital system of coordinates related to the body P 2 by V r2 , V τ 2 , V n2 , we obtain the radial, transversal, and normal components of the force ⃗ F 2 in the form where Note, that the relative velocities ⃗ V 0 in (3.11) of the particles leaving the body P 0 or falling on it are given in the orbital system of coordinates related to the body P 1 .If the relative velocities ⃗ V 0 , ⃗ V 1 , and ⃗ V 2 and laws of variation of body masses are given, Equations (3.5)-(3.11)completely determine the perturbed motion of the bodies P 1 , P 2 .

Small eccentricities and inclinations
It is quite obvious that exact solution to nonlinear differential equations (3.5)-(3.6)cannot be obtained and one can try to find only approximate solutions.
Note that in many problems of celestial mechanics, eccentricities and inclinations of body orbits are small (see [6,19]).Here we consider this practically important case of small eccentricities e j << 1 and inclinations i j << 1, (j = 1, 2) and expand the right-hand sides of Equations (3.5)-(3.6) in series in these parameters up to first-order terms.
First, we can find approximate solution to the Kepler equation (2.8) and represent the eccentric anomaly in the form of a converging series in e j (see [19]) E j = M j + e j sin M j + . . . .Using this solution, we obtain cos On substituting expansions (4.1) into solutions (2.5) and expanding the expressions obtained in series in small parameters, we obtain x j =a j γ j cos(M j +ω j +Ω j )+ ej 2 (cos(2M j +ω j +Ω j )−3 cos(ω j +Ω j )) , y j =a j γ j sin(M j +ω j +Ω j )+ ej 2 (sin(2M j +ω j +Ω j )−3 sin(ω j +Ω j )) , z j = a j γ j i j sin(M j + ω j ), j = 1, 2. (4.2) Using (4.2), we find The distance between the bodies P 1 and P 2 may be written then as where the mean longitude λ j = M j + ω j + Ω j has been introduced instead of the mean anomaly M j for the convenience of computations (see [19]), and Since ρ 0 is a periodic function of the variables λ 1 and λ 2 (see (4.4)), the expressions 1/ρ 0 , 1/ρ 3 0 , 1/ρ 5 0 which will appear in the expansion of r −3 12 in small parameters (see (3.9), (3.10)), may be replaced by the corresponding Fourier series where A k , B k , and C k are the Laplace coefficients satisfying the recurrences (see [7,19]) All the Laplace coefficients can be computed using the above recurrences and the following expressions for A 0 and A 1 : where the functions denote the complete elliptic integral of the first and second kinds, respectively, and the parameter α = a1γ1 a2γ2 < 1.The body P 2 is assumed to be an outer planet and the trajectory of body P 1 is located inside of the trajectory of body P 2 .
On substituting the expansions (4.1)-(4.5)into (3.9)-(3.11),we compute the expansions of the right-hand sides of Equations (3.5)-(3.6) in powers of eccentricities e 1 , e 2 and inclinations i 1 , i 2 .The coefficients of these expansions are periodic functions of mean longitudes λ 1 , λ 2 , and they are rational expressions the numerators of which include the trigonometric functions cos(kλ j ), sin(kλ j ), cos(kλ 1 ± nλ 2 ), and sin(kλ 1 ± nλ 2 ), (k, n = 1, 2, ...).These expressions are quite bulky and so we do not write them here.Since we are interested in the behaviour of the orbital elements on long time intervals, the terms on the right-hand sides of Equations (3.5)-(3.6)determining the short-term oscillations of the orbital elements can be eliminated by averaging the equations over the mean longitudes λ 1 and λ 2 (see [6,7,19]).We assume that the mean-motion resonances are absent in the system and the masses m 0 (t), m 1 (t), m 2 (t) of the bodies and velocities ⃗ Recall that averaging of the function W (λ 1 , λ 2 ) over the variables λ 1 and λ 2 and transition to the secular perturbations W (sec) is reduced to calculating the integral By substituting for the function W (λ 1 , λ 2 ) the right-hand sides of Equations (3.5)-(3.6) in which expansions in the small parameters are made and taking into account the expressions (3.9)-(3.11)for the forces and expansions (4.1)-(4.6),we obtain the following differential equations: Equations (4.7) determine the secular perturbations of the orbital elements of the bodies P 1 and P 2 .We do not write the averaged Equation (3.6) here because due to the integration of Equations (3.5)-(3.6)with respect to the mean longitudes an information about the location of the bodies in the orbits is lost and we can analyze only slow changes of the orbital parameters a j , e j , i j , Ω j , and ω j in time.

Numerical solutions to evolution equations
Although the averaged Equations (4.7) are approximation of (3.5)- (3.6) accurate to the first order in eccentricities and inclinations their general solution cannot be found in symbolic form.In order to investigate an influence of masses change on the dynamic evolution of the system we can choose some realistic values for the system parameters and solve Equations (4.7) numerically.To simplify the calculations it is convenient to use the dimensionless variables.For example, we use initial values of the semi-major axis a 10 = a 1 (t 0 ) and the mass m 00 of body P 0 as units of distance and mass, respectively, and define dimensionless distance a * j , mass m * j and time t * by , j = 0, 1, 2.
The masses variation are described by the Eddington-Jeans law where for the bodies P 0 , P 1 , and P 2 we choose, respectively To be able to test the model we consider the Sun, Jupiter, and Saturn as bodies P 0 , P 1 and P 2 , respectively, and choose the following initial values for orbital elements (see [19]):  In the case of constant masses of the bodies the system (4.7)describes the secular perturbations of the orbital elements in the framework of the classical three-body problem and its solutions correspond to the known results (see [6,7,19]).Taking into account the isotropic masses variation according to the Eddington-Jeans law when reactive forces do not arise results in only some quantitative changes of solutions to (4.7) (see Figure 1).The semi-major axes a 1 and a 2 remain constant while the period of oscillations of the eccentricities, inclinations and longitudes of the ascending nodes increase as the masses decrease.If only one component of the relative velocity V r0 of the particles leaving the most massive body P 0 along the radial direction becomes greater than zero (V r0 = 1) dependance of the eccentricities e 1 , e 2 and arguments of pericenter ω 1 , ω 2 on time changes (see Figure 2).However, the corresponding component of the reactive force does not influence the other orbital elements.
Solving the system (4.7) in the case of V n0 = 1, V r1 = −1, V τ 2 = 1, when reactive forces along the radial, transversal and normal directions arise demonstrates noticeable changes in evolution of the orbital elements (see Figure 3).Period of the eccentricity oscillations decreases in comparison to the case of absence of the reactive forces while the inclinations undergo additional oscillations with greater period.Note that numerical solutions to the system (4.7) and visualization of the results are performed with the aid of the system Wolfram Mathematica.

Conclusions
In this paper, we investigated a non-stationary three-body problem for bodies of variable masses that attract each other according to Newton's law of gravitation taking into account the reactive forces arising due to anisotropic variation of the bodies masses.The original equations of motion of the bodies in the relative system of coordinates are obtained in the framework of Newton's formalism, which makes it possible to write the reactive forces on the basis of Meshcherskii equation.Using the exact solutions of the non-stationary twobody problem (see [16]) and applying the method of variation of constants, we derived differential equations of the perturbed motion in terms of osculating elements of the aperiodic motion along quasi-conical section.It should be emphasized that the obtained Equations (3.5)-(3.6)are valid for any laws of the mass variation of the bodies and completely determine the perturbed motion of the bodies P 1 , P 2 .
In the case of small eccentricities and inclinations of orbits, we have expanded the right-hand sides of Equations (3.5)-(3.6) in power series in terms of the orbital elements up to the first order.As the coefficients of e 1 , e 2 and i 1 , i 2 in the obtained expressions are periodic functions of the mean longitudes λ 1 , λ 2 , we replaced them by the corresponding Fourier series.Finally, we have shown that the right-hand sides of differential equations (3.5)-(3.6)contain the terms describing behaviour of the orbital elements on long time intervals and quite cumbersome terms determining the short-term oscillations of the orbital elements.Assuming that the mean-motion resonances are absent in the system and averaging the equations over the mean longitudes λ 1 , λ 2 , we derived differential equations determining the secular perturbations of the orbital elements.Note that the equations obtained describe the perturbed motion of the bodies in the general case when the masses of all three bodies vary anisotropically, and reactive forces occur.
To test the model, we have solved the averaged Equations (4.7) numerically for some realistic values of the system parameters and some laws of the masses variations, the obtained results are presented on Figures 1-3.Comparison with the case of constant masses which is well-known (see, for example, [19]) demonstrates that masses variation can significantly affect the evolution of orbital parameters.In the next paper we plan to use the model proposed to investigate numerically some real two-planetary systems of three non-stationary bodies and to investigate an influence of the masses variation on their evolution.
Note that all symbolic and numerical calculations were carried out using Wolfram Mathematica (see [29]).

. 4 )
Using (2.7), expansions (4.1)-(4.3),and the expressions for the direction cosines (3.7)-(3.8),we find the scalar products of the units vectors appearing in the expressions for the components of the forces ⃗ F 1 and ⃗ F 2 along the radial, transversal, and normal directions (see (3.9)-(3.11)).They are given by Equations (3.5)-(3.6)change very slowly with time and the procedure of averaging does not change them.
Long-term evolution of the eccentricities e 1 , e 2 and inclinations i 1 , i 2 (solid curves -constant masses, dashed curves -isotropic mass changes, dotted curvesnon-isotropic mass changes,