## 1 Introduction  
The classical two-body problem describes the motion of two particles \(P_{0}\), \(P_{1}\) of constant masses \(m_{0}\), \(m_{1}\), respectively, attracting each other according to [PERSON]'s law of universal gravitation. It was first formulated and solved by [PERSON] and has many applications in celestial mechanics and cosmic dynamics (see [1, 2, 3]). The main reason for its importance is a complete and general solution to this problem may be used as a first approximation in a wide variety of practical orbital motion problems. For example, a planet's motion around a star or a spacecraft's motion about a planet is described, to the first approximation, by the corresponding solution of the two-body problem. Although the attraction of other bodies disturbs this motion and forces its orbital parameters to change, the application of perturbation theory, which has been developed quite well, enables one to investigate the system's motion accurately (see, for example, [4, 5, 6]). Note that some characteristics of celestial bodies such as mass, shape, size, and internal structure may vary with time (see, for example, [7, 8, 9]). For instance, the mass of a star decreases considerably during its lifetime due to corpuscular and light emission, and this mass loss may play an important role in the evolution of the two-body system (see [10, 11, 12, 13, 14]). To describe the motion of celestial bodies more precisely and to make the model more realistic, one needs to take into account the variability of the system parameters. Similar to the case of constant masses, to solve the many-body problem with variable masses, it is expedient to start from the two-body problem. For simplicity, one can assume that each body of the system may eject mass isotropically in all directions; the velocity of ejection is high enough for the particles to escape from the system, and their interaction with the bodies may be ignored. Hence, the only effect of the loss of mass is the decrease in attraction between the two bodies. Then, the relative motion of one body around the other may be described by the following differential equation:  
\[\frac{d^{2}\vec{r}}{dt^{2}}+Gm(t)\frac{\vec{r}}{r^{3}}=0, \tag{1}\]  
where \(G\) is the gravitational constant, \(\vec{r}\) is the radius-vector of the body \(P_{1}\) relative to \(P_{0}\), and \(m(t)=m_{0}(t)+m_{1}(t)\) is the total mass of the system, which is a function of time. Equation (1) was proposed first by [PERSON] [15] and is known now as Gylden's problem. Later, [PERSON] [16] explained the physical meaning of this equation and showed its integrability for some specific mass variation laws. Although many different models of mass variation in the two-body problem have been proposed (see [17]), it turned out that an exact solution to Equation (1) can be found only for some specific functional dependence \(m(t)\) (see [18; 19; 20; 21; 22; 23]). Note that by choosing some realistic law of the total mass variation \(m(t)\) (see, for example, [24; 25]) and initial conditions, one can find the corresponding numerical solution to Equation (1). However, such a solution cannot be obtained for all possible initial conditions, and so it does not allow one to obtain some general results, for example, to investigate the stability of the system motion. In the case of constant masses, a general solution to Equation (1) is known and describes the motion of bodies on conic sections, the orbital elements of which depend on the total mass \(m(t)=const\) (see [2; 3]). Therefore, it appears reasonable that variability of masses may disturb this motion and enforce the orbital elements to change with time. To derive the differential equations determining the behavior of the orbital elements, one needs to reduce Equation (1) to a form that is appropriate for the application of perturbation theory. There are two approaches that are used most often in studying the many-body problem with variable masses.