# Analytical solution of the two-body problem against the cosmic vacuum background

N. Emelyanov, M. Yu. Kovalyov

► To cite this version:

N. Emelyanov, M. Yu. Kovalyov. Analytical solution of the two-body problem against the cosmic vacuum background. Monthly Notices of the Royal Astronomical Society, 2013, 429 (4), pp.3477-3491. 10.1093/mnras/sts619. hal-02542165

HAL Id: hal-02542165

https://hal.science/hal-02542165

Submitted on 28 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical solution of the two-body problem against the cosmic vacuum background

N. V. Emelyanov<sup>1,2★</sup> and M. Yu. Kovalyov<sup>3</sup>

Accepted 2012 December 12. Received 2012 October 23

## ABSTRACT

Motions in the two-body problem against the background of the cosmic vacuum are considered. The potential function of the problem is taken in the form $U = \mu/r + H^2r^2/2$, where $\mu$ is the gravitational parameter, $H$ is the effective Hubble constant, which represents cosmic vacuum and $r$ is a distance from attracting centre. Precise analytical solution of the equations of motion is obtained. The formulae for calculating a moving particle's coordinates for any instance of time with given initial conditions are derived. The chain of the formulae contains a transcendent equation analogous to the Kepler's equation in the two-body problem. The formulae are expressed using incomplete elliptic integrals of the first and third kind.

It is shown that finite motions take place only within a circle with a radius $\bar{r} = (\mu/H^2)^{1/3}$ centred in an attracting body. In this case, the minimum distance cannot exceed the value of $r^* = (\mu/4H^2)^{1/3}$. Circular motions can take place in the area of finite motions. The circular motions with the radii in the range $(r^*, \bar{r})$ are orbitally unstable. This means that the slightest variations in the initial conditions turn the circular motion either into finite motion with a great amplitude of distance's variation or into infinite motion. The latter can exist for any value of minimum distance.

Numerical integration of the equations of motion has been carried out for a number of instances. It is shown that, if the equations are solved for the rectangular coordinates, numerical integration can give erroneous solutions if the trajectory approaches the sphere with the radius $\bar{r}$. This is a result of the accuracy's loss in the case when two close values are subtracted. A change of variables in the equations of motion is found which gives no loss of accuracy so that numerical integration can be carried out with the same accuracy for any trajectory of motions.

Comparison of analytical solution with the numerical integration is made for a number of sets of initial conditions. It gives complete coincidence of solutions at the time intervals covering several dozens of revolutions of the particle around the attracting centre.

**Key words:** gravitation–galaxies: evolution–galaxies: general–galaxies: kinematics and dynamics–cosmology: theory.

## 1 INTRODUCTION

Nowadays, it is an obvious fact that all the bodies of the Universe are subjected to the mechanical influence of some additional force not related to the action of other bodies. This influence was predicted by *Einstein* and was confirmed by astronomical observations of the last few decades (Riess et al. 1998; Perlmutter et al. 1999). The physical cause of the phenomenon was called cosmic vacuum. It is this physical agent that gives rise to the acceleration of the observed

expansion of the Universe. The first models of the expansion of the Universe were elaborated by de Sitter (1917). Friedmann (1922, 1924) took, in addition to the cosmic vacuum, into consideration non-relativistic matter. But the both treated the Universe as a whole and the motions of individual galaxies were not considered. Probably it is the reason why, since de Sitter and Friedmann, it is supposed that the motions of galaxies in the cosmic vacuum field take place only in radial directions, tacitly assuming that there are no motions in lateral directions. This point of view is also encouraged by the absence of observable motions of galaxies in the picture plane. But it does not mean that there are no such motions.

<footnote n="1"><sup>1</sup>Sternberg Astronomical Institute, M. V. Lomonosov Moscow State University, 13 Universitetskij prospect, 119992 Moscow, Russia</footnote>

<footnote n="2"><sup>2</sup>Institut de Mécanique Céleste et de Calcul des Éphémérides, Observatoire de Paris, UMR 8028 du CNRS, 77 avenue Denfert-Rochereau, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Faculty of Physics, M. V. Lomonosov Moscow State University, Leninskie Gory, Moscow 119991, Russia</footnote>

<footnote>* E-mail: emelia@sai.msu.ru</footnote>
In fact, the study of motions of separate galaxies has begun only in recent years. For instance, the papers by Chernin (2001) and Dolgachev, Domozhilova & Chernin (2003) study the motions of the Local Group galaxies in the cosmic vacuum field influenced by the gravitational forces using the model of planar rectilinear restricted three-body problem.

Strict description of any motions in the Universe influenced by the gravitational forces and cosmic vacuum is possible only using the equations of Einstein's general relativity theory. But there is a possibility of studying the influence of cosmic vacuum using only the equations of Newtonian mechanics by adding to them terms which approximately describe the influence of cosmic vacuum upon all the bodies of the Universe. This kind of methodology is used in celestial mechanics to take into consideration relativistic effects in the equations of motion of celestial bodies. Obviously, the model is approximate but, while its accuracy is no worse than that of observations, the model can be used to solve both straight and reverse practical problems of the mechanics of the bodies of the Universe. The matching of the modelled motions with those obtained from observations can serve as criterium for the model's accuracy.

Lyk'yanov (2009) studied all possible motions of galaxies including those in lateral directions in the repulsion field created by cosmic vacuum assuming that there are no other forces. In other words, the problem of the motions of many bodies in cosmic vacuum field was considered. Strict analytical solution of this problem was obtained.

The next step in modelling the motions of the bodies of the Universe within the frames of the simplified model of cosmic vacuum is studying the motion of two bodies influenced by both their gravitational attraction and cosmic vacuum. Such an approach has already been applied by Nandra, Lasenby & Hobson (2012) where some features of the motions were obtained using both Newtonian mechanics and general relativity theory. However, the exact analytical solution of the equations of motion in such a statement of the problem has not yet been derived.

In this paper, we deduce the exact analytical solution of the two-body problem against the cosmic vacuum background where the influence is determined by the force depending solely upon the distance between bodies. The dependence of the force on the distance was taken in the same form as that used by Lyk'yanov (2009) and Nandra et al. (2012). In this way, we get a mechanical model involving two material points (particles) with finite masses whose motions are influenced by their mutual Newtonian attraction and by an additional force acting along the line connecting the particles, constantly directed opposite to the origin of the reference frame with a value of $H^2r$, where $H$ is an effective Hubble constant defined only by the cosmic vacuum's density and $r$ is the distance to the origin of coordinates. The origin is placed in one of the moving particles. Since the model involves forces depending solely on the distance, it is easy to prove that the particles move in a constant plane. Let us place on this plane a rectangular system of coordinates $(x, y)$ centred in one of the moving particles. Then, the motion of another particle in this system of coordinates is described by the solution of the following system of differential equations:

$$ \frac{d^2x}{dt^2} = -\frac{\mu}{r^3}x + H^2x, \quad \frac{d^2y}{dt^2} = -\frac{\mu}{r^3}x + H^2y, \quad (1) $$

where $\mu$ is a gravitational parameter defined by the formula

$$ \mu = G(M_1 + M_2), $$

where $G$ is the universal gravitational constant, $M_1$ and $M_2$ are the masses of the particles.

Now we can define the problem more exactly. It is necessary to obtain the general solution of the equations of motion, i.e. the coordinates $x, y$ of the moving particle, as explicit functions of time and arbitrary constants. These functions can be given by a chain of formulae where transcendent equations are permissible analogous to the Kepler equation in the classical two-body problem.

We leave outside of this paper the question of how accurately such equations of Newtonian mechanics describe the motions of two gravitating particles in the presence of the cosmic vacuum. Moreover, it is obvious that the real Universe has more objects than just two material points and cosmic vacuum, and these objects are not material points since galaxies are pretty extensive formations. Nevertheless, we neglect the presence of other attracting bodies as well as their shapes and dimensions. It is left for further researches to estimate the errors resulting from the adopted assumptions. We suppose that the model in consideration may either serve as an approximate description of real motions or be used as an initial approximation in solving the problem if it is stated more accurately.

On the other hand, the mechanical model in consideration may be studied independently from its direct applications. This process can be useful for both improvement of the methods and accumulation of our knowledge.

Solving the equations of motion in consideration can be carried out in three different ways. The first one is to obtain the general analytical solution depending on four arbitrary constants. This method provides formulae for direct calculations of a moving particle's coordinates for any instance of time if initial conditions are given. In this case, the solution can be derived with any pre-assigned accuracy. Errors do not depend on a time interval between initial and final positions. The second method is numerical integration of the equations of motion. For given initial conditions, it gives particular solution successively for a series of instances at a limited time interval. But the methods of numerical integration have two essential drawbacks. First, the expansion of time interval results in the increase of errors. Secondly, the accuracy of the solution is always limited, whatever perfect method is used. The third method of solving the equations in consideration is based on treating the motions as perturbed Keplerian motions. The theory of perturbations provides us with an approximate solution on the assumption that the cosmic vacuum's influence is small compared to the mutual gravitational attraction of the bodies. A feature of perturbed Keplerian motion in this problem is that zero initial eccentricities always result in non-zero perturbations in eccentricities. Thus, we see an apparent contradiction to the fact of existence of circular motions in this model. In fact, circular motions are provided here not by non-zero eccentricity, as it is in many other problems, but by the fact that, if the eccentricity is considerable, the mean anomaly stops changing at all so that circular motion is produced by the secular variation of the argument of latitude, i.e. by the apse line motion.

In this paper, we deduce the exact analytical solution of the two-body problem against the cosmic vacuum background and check the solution by comparing it with the results of the numerical integration of the equations of motion. The perturbed Keplerian motion is not examined here.

## 2 REDUCTION OF THE SOLUTION TO QUADRATURES

The equations of motion can be written in the form

$$ \frac{d^2x}{dt^2} = \frac{\partial U}{\partial x}, \quad \frac{d^2y}{dt^2} = \frac{\partial U}{\partial y}, $$
where

$$U = \frac{\mu}{r} + \frac{H^2}{2} r^2$$

is the potential function.

Let us examine some features of these equations and their solutions.

The equations of motion have the following energy integral:

$$\left(\frac{dx}{dt}\right)^2 + \left(\frac{dy}{dt}\right)^2 = \mu \left(\frac{2}{r} + \frac{H^2}{\mu} r^2\right) + 2\alpha_1$$

and the area integral

$$x \frac{dy}{dt} - y \frac{dx}{dt} = 2\alpha_2,$$

where $\alpha_1$ and $\alpha_2$ are arbitrary constants.

There are stationary solutions corresponding to the moving particle's positions on the circle with the radius

$$\bar{r} = \left(\frac{\mu}{H^2}\right)^{\frac{1}{3}}$$

centred in the origin of coordinates. This is the stationary solution circle. It can also be called zero acceleration circle since any motion of the particle on this circle has zero acceleration.

Inside the zero acceleration circle, there are circular motions corresponding to the particular solutions of the equations. If $r = r_c \le \bar{r}$, the linear velocity $V_c$ of the circular motion is defined by the relationship

$$V_c = \sqrt{\frac{\mu}{r_c} - H^2 r_c^2}.$$

On the zero acceleration sphere $V_c = 0$. The circular motion has the period of revolution

$$T_c = \frac{2\pi r_c}{V_c}.$$

As will be shown further, finite motions, i.e. the motions with distances always less than some finite value, are possible only within the zero acceleration circle. So, in this problem, it is quite reasonable to choose the radius of the zero acceleration circle $\bar{r}$ as a unit for the measurements of distances.

Let us reduce the equations of motion to canonical form. For this, we replace the sought functions x, y by distance r and latitude $\lambda$ using the formulae

$$x = r \cos \lambda, \quad y = r \sin \lambda.$$

Treating the new variables as generalized coordinates, we introduce corresponding generalized impulses

$$p_r = \frac{dr}{dt}, \quad p_\lambda = r^2 \frac{d\lambda}{dt}.$$

Now the canonical form of the equations of motion can be written as follows:

$$\frac{dr}{dt} = \frac{\partial F}{\partial p_r}, \quad \frac{dp_r}{dt} = -\frac{\partial F}{\partial r},$$

$$\frac{d\lambda}{dt} = \frac{\partial F}{\partial p_\lambda}, \quad \frac{dp_\lambda}{dt} = -\frac{\partial F}{\partial \lambda},$$

where the Hamiltonian is

$$F = \frac{1}{2} \left( p_r^2 + \frac{p_\lambda^2}{r^2} \right) - \frac{\mu}{r} - \frac{H^2}{2} r^2.$$

Since the Hamiltonian does not explicitly depend on time, the canonical system has the first integral

$$F = \alpha_1,$$

where $\alpha_1$ is an arbitrary constant. Hence, we have

$$\alpha_1 = \frac{1}{2} V^2 - \frac{\mu}{r} - \frac{H^2}{2} r^2, \qquad (2)$$

where V is the modulus of the velocity.

In the particular case of circular motion we have $V = V_c$, $r = r_c$, and the constants $V_c$, $r_c$, and $\alpha_1$ are interrelated in the following way:

$$V_c^2 = \frac{\mu}{r_c} - H^2 r_c^2, \quad \alpha_1 = -\frac{1}{2} \frac{\mu}{r_c} - H^2 r_c^2.$$

The corresponding Hamilton–Jacobi equation is

$$\frac{\partial S}{\partial t} + \frac{1}{2} \left[ \left( \frac{\partial S}{\partial r} \right)^2 + \left( \frac{\partial S}{\partial \lambda} \right)^2 \right] - \frac{\mu}{r} - \frac{H^2}{2} r^2,$$

and its full integral can be expressed in the form

$$S(t, r, \lambda) = -\alpha_1 t + \alpha_2 \lambda + \int_{r_0}^{r} \sqrt{2\alpha_1 r^2 + 2\mu r + H^2 r^4 - \alpha_2^2} \frac{1}{r} dr,$$

where $\alpha_1$, $\alpha_2$ are arbitrary constants, $r_0$ is the minimum value of the distance r in any particular solution. Then the full integral can be defined using the relationships

$$\frac{\partial S}{\partial \alpha_1} = \beta_1, \quad \frac{\partial S}{\partial \alpha_2} = \beta_2, \qquad (3)$$

$$\frac{\partial S}{\partial r} = p_r, \quad \frac{\partial S}{\partial \lambda} = p_\lambda. \qquad (4)$$

where $\beta_1$, $\beta_2$ are two more arbitrary constants.

The first two relationships give us

$$\int_{r_0}^{r} \frac{r \, dr}{\sqrt{2\alpha_1 r^2 + 2\mu r + H^2 r^4 - \alpha_2^2}} = t + \beta_1,$$

$$\alpha_2 \int_{r_0}^{r} \frac{1}{\sqrt{2\alpha_1 r^2 + 2\mu r + H^2 r^4 - \alpha_2^2}} \frac{1}{r} dr = \lambda - \beta_2.$$

Here the dependence of coordinates on time is given by quadratures, so that it is impossible to directly calculate coordinates at given instances using these formulae. The deduction of coordinates as explicit functions of time from these relationships is called inversion of quadratures. The integrals given above are expressed using incomplete elliptic integrals of the third kind. The inversion of such integrals via known functions is impossible.

To carry out inversion of quadratures, we introduce auxiliary variable u using the differential relationship

$$\frac{dr}{\sqrt{2\alpha_1 r^2 + 2\mu r + H^2 r^4 - \alpha_2^2}} = du.$$

Then, we have

$$dt = r \, du, \quad d\lambda = \frac{\alpha_2}{r} \, du.$$

First, we should make inversion of the quadrature

$$\int_{r_0}^{r} \frac{dr}{\sqrt{2\alpha_1 r^2 + 2\mu r + H^2 r^4 - \alpha_2^2}} = u - u_0, \qquad (5)$$

where $u_0$ is an arbitrary constant. In result, we obtain the function $r(u)$.
After that, we obtain the function $u(t)$ via inversion of the integral

$$t - t_0 = \int_{u_0}^{u} r(u) du,$$

and the function $\lambda(u)$ is found by integrating the relationship

$$\lambda - \lambda_0 = \int_{u_0}^{u} \frac{\alpha_2}{r(u)} du,$$

where $\lambda_0$ is an arbitrary constant. If $t = t_0$, we have $u = u_0$, $\lambda = \lambda_0$.

Inversion of the quadrature (5) can be performed using elliptic functions. We rewrite this relationship in the form

$$\int_{r_0}^{r} \frac{dr}{\sqrt{W(r)}} = u' - u'_0,$$

where

$$W(r) = r^4 + 2 \frac{\alpha_1}{H^2} r^2 + 2 \frac{\mu}{H^2} r - \frac{\alpha_2^2}{H^2}, \quad u' = Hu.$$

Obviously, motions are possible only if $W(r) \ge 0$. Hence the roots of the equation $W(r) = 0$ define areas of possible motions. The roots depend on parameters $\mu$ and $H$ as well as on arbitrary constants $\alpha_1$ and $\alpha_2$, that is, on initial conditions. Analysis of the areas of possible motions is made further.

The fact that such differential equations can be integrated is not a new result. Moreover, these equations turn out to be particular case of a more general type of equations examined by Poleshchikov & Zhubr (2008). Not only did these authors demonstrate that the equations admit integration but they also carried out inversion of quadratures. The potential function of the problem solved by Poleshchikov & Zhubr (2008) has the form

$$U = \frac{\mu}{r} + \frac{1}{r} \sum_{k=1}^{N} [A_k(r + b^T x)^k + B_k(r - b^T x)^k],$$

where $x = \{x, y, z\}$ is the position vector of a particle, $b$ is an arbitrary constant unit vector, $A_k, B_k (k = 1, 2, \dots, N)$ are arbitrary constants. This potential function is reduced to that of our problem if we put $b = 0, A_3 = B_3 = H^2/4$ and turn the remaining constants $A_k, B_k (k = 1, 2, 4, \dots, N)$ to zero. Poleshchikov & Zhubr (2008) have proved that the equations admit integration for any value of $N$ and carried out the inversion of quadratures for $N \le 3$, which includes our particular case.

There is a principal possibility to use the results of Poleshchikov & Zhubr (2008) to solve our problem. But their formulae representing coordinates as functions of time have rather general and lengthy form, so that one cannot deduce from them all the features of the motions. Their reduction to the particular case of our problem is rather complicated task. So we carried out inversion of quadratures and obtained the sought formulae directly from our equations of motion. This way seems to be more simple and illustrative.

## 3 EXAMINATION OF THE AREAS OF POSSIBLE MOTIONS

Before solving the equations of motion, we can examine some important features of possible motions in our problem. For this, suffice it to make analysis of the roots of the equation

$$W(r) = 0.$$

The values of the roots are necessary when constructing solution by inversion of quadratures.

It is easy to demonstrate that the roots admit two possibilities. The first one is the case of four real roots, the second one is that

<figure><img src="image_2.png" /><figcaption>Figure 1. Two variants of the function $W(r)$. For the function $W_1(r)$, the equation $W_1(r) = 0$ has four real roots. The function $W_2(r)$ has two real roots and two complex conjugate roots.</figcaption></figure>

when two roots are real and another two are complex conjugate. If all four roots are real, only one of them is negative. Moreover, multiple roots are also possible. Fig. 1 shows two variants of the function $W(r)$. The first one, $W_1(r)$, demonstrates the case when all four roots are real, while the function $W_2(r)$ has only two real roots.

### 3.1 The case of four real roots

Let us first consider the case of four real roots. We designate them as $r_{-1}, r_0, r_1, r_2$ ($r_{-1} \le r_0 \le r_1 \le r_2$).

As was noticed earlier, we exclude from consideration the case of rectilinear motion when $\alpha_2 = 0$ and the case when two bodies collide. Then it is obvious that $r_{-1} < 0$ and $0 < r_0 \le r_1 \le r_2$. The motions are possible only if $r_0 \le r \le r_1$ and if $r \ge r_2$. In the first case we have finite motions, in the second case the particle can move away at infinite distance.

The relationship between the roots of the equation $W(r) = 0$ and constants $\alpha_1, \alpha_2$ can be obtained from the expression

$$r^4 + 2 \frac{\alpha_1}{H^2} r^2 + 2 \frac{\mu}{H^2} r - \frac{\alpha_2^2}{H^2} = (r - r_{-1})(r - r_0)(r - r_1)(r - r_2).$$

Since the polynomial $W(r)$ has no terms of the third degree, we have a simple formula describing the relationship of the roots:

$$r_{-1} + r_0 + r_1 + r_2 = 0.$$

If only finite motions are considered, the roots $r_0$ and $r_1$, which set the limits of possible values of the distance, can be taken as initial arbitrary constants instead of $\alpha_1$ and $\alpha_2$. Then $\alpha_1$ and $\alpha_2$ are defined through $r_0$ and $r_1$ using the formulae

$$\alpha_1 = -\frac{1}{2} H^2 (r_0^2 + r_1^2) - \frac{\mu}{r_0 + r_1}, \quad (6)$$

$$\alpha_2^2 = \frac{2\mu r_0 r_1}{r_0 + r_1} - H^2 r_0^2 r_1^2. \quad (7)$$
In addition, we have

$$ \alpha_1 = -\frac{1}{2}H^2(r_1^2 + r_2^2) - \frac{\mu}{r_1 + r_2}. \qquad (8) $$

The roots $r_{-1}$, $r_2$ are defined through $r_0$, $r_1$ using the relationships:

$$ r_{-1} = -\frac{1}{2}(r_0 + r_1) - \sqrt{\frac{2\mu}{H^2(r_0 + r_1)}} + \frac{1}{4}(r_1 - r_0)^2, $$

$$ r_2 = -\frac{1}{2}(r_0 + r_1) + \sqrt{\frac{2\mu}{H^2(r_0 + r_1)}} + \frac{1}{4}(r_1 - r_0)^2. $$

In finite motions, when the distance $r$ reaches either the minimum value $r_0$ or the maximum value $r_1$, the velocity of the particle is directed perpendicular to the position vector. For our investigation it is necessary to define the velocities at these points of the trajectory with the given values $r_0$ and $r_1$. The velocity $V_0$ at the instances when the distance is at its minimum value $r_0$ can be calculated using the relationship

$$ V_0^2 = 2\mu \frac{r_1}{r_0(r_0 + r_1)} - H^2 r_1^2. $$

Likewise, the velocity $V_1$ corresponding to the maximum distance $r_1$ can be obtained from the relationship

$$ V_1^2 = 2\mu \frac{r_0}{r_1(r_0 + r_1)} - H^2 r_0^2. $$

For further analysis it is convenient to introduce the designations:

$$ \bar{r} = \left(\frac{\mu}{H^2}\right)^{\frac{1}{3}}, \quad r^* = \left(\frac{\mu}{4H^2}\right)^{\frac{1}{3}}. $$

Possible values of $r_0$ are confined to the interval

$$ 0 < r_0 \le r^*. $$

Possible values of the root $r_1$ are defined by the inequalities

$$ r_0 \le r_1 \le r_{\max}(r_0), $$

where $r_{\max}(r_0)$ is obtained from the cubic equation

$$ r_{\max} \left( r_0^2 + 2r_0 r_{\max} + r_{\max}^2 \right) = \bar{r}. $$

The value $r_{\max}(r_0)$ turns out to be the lower limit of the values of the root $r_2$ which is reached if $r_1 = r_2$. In any case, $r^* \le r_2$.

Let us examine the cases of the multiple roots. Let $r_1 = r_0$. Then,

$$ r_2 = \sqrt{\frac{\mu}{H^2 r_1}} - r_1. $$

In this case, from the obvious condition $r_2 \ge r_1$ we have

$$ r_1^3 \le \frac{\mu}{4H^2} = r^{*3}, $$

which means that the multiple root $r_1 = r_0$ can take values only in the range

$$ 0 < r_1 = r_0 < r^*. $$

Now we suppose that $r_1 = r_2$. In this case, the relationship between $r_1$ and $r_0$ is defined by equation

$$ r_1 (r_0^2 + 2r_0 r_1 + r_1^2) = \bar{r}^3. $$

The multiple root $r_1 = r_2$ can take values only in the range $[r^*, \bar{r}]$.

In the case of three equal roots, we have

$$ r_0 = r_1 = r_2 = r^*. $$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Plots of the function $W(r)$ in the case of circular motion when $r_1 = r_0$, and in the case when motion is close to circular when the distance varies in a small interval.</figcaption></figure>

Let us now consider the case of circular motions. Here, we have

$$ \alpha_1 = -\frac{1}{2} \frac{\mu}{r_c} - H^2 r_c^2, $$

where $r_c$ is the radius of the circular orbit. Taking into consideration the equality (6) we conclude that the circular motions take place under the condition

$$ r_0 = r_1 = r_c, $$

where

$$ 0 < r_0 = r_1 = r_c \le r^*. $$

Likewise, from the equality (8) we conclude that the circular motions can also take place if

$$ r_1 = r_2 = r_c, $$

where

$$ r^* \le r_1 = r_2 = r_c < \bar{r}. $$

Depending on the values of the radius $r_c$ belonging either to the interval $(0, r^*)$ or to the interval $(r^*, \bar{r})$, the circular motions have essentially different features. Let us first consider the case where $0 < r_c \le r^*$. Fig. 2 shows the plot of the function $W(r)$ for one particular instance of the circular motion and for the instance which is sufficiently close to this circular motion. When $r_0 \le r \le r_1$, but the difference $r_1 - r_0$ is small enough, the motion does not too much differ from the circular. It takes place between two close circles with the radii $r_0$ and $r_1$.

Let us now consider the circular motions where $r^* < r_c \le \bar{r}$. We show in Fig. 3 the plot of the function $W(r)$ corresponding to the circular motion with $r_1 = r_2 = r_c$ as well as the plot quite close to the former corresponding to the non-circular motion for which $r_1 < r_2$. The plots demonstrate that the area of the motions previously restricted to the circle with the radius $r = r_c$ now has split into two areas. The first of them contains finite motions with quite differing minimum and maximum distances. The second one
<figure><img src="image_4.png" /><figcaption>Figure 3. Plots of the function $W(r)$ in the case of circular motion when $r_2 = r_1$, and in the case close to the latter when the difference $r_2 - r_1$ is small and two types of motions are possible: the first is finite with a great range of variations in the distance, the second is infinite.</figcaption></figure>

is infinite. Thus, we can conclude that, in this case, the slightest variation in arbitrary constants leads to significant transformation of circular motion into non-circular, i.e. circular motions in the area $r^* < r_c \le \bar{r}$ are orbitally unstable.

If all the roots of the equation $W(r) = 0$ are real, infinite motions are possible for $r > r_2$.

## 3.2 The case of two real and two complex conjugate roots

If two roots are complex conjugates, only infinite motions are possible. In these motions, the minimum distance is equal to the real root $r_2$. It is easy to show that at the instance when $r = r_2$ the velocity vector is perpendicular to the position vector. Let $V_2$ be the value of this velocity. Then, to simplify the process of the examination of solutions, we set initial conditions via arbitrary constants $r_2$ and $V_2$.

Later (see below), we shall examine the way in which the parameters $r_2$ and $V_2$ can be obtained on the basis of arbitrary initial conditions, when the particle, at initial moment of time, has arbitrary values of distance and velocity.

If $r_2$ and $V_2$ are given, the constant $\alpha_1$ can be calculated using (2). If we know the values of $\alpha_1$ and of one of the real roots of the equation (let it be the root $r_2$), other roots can be obtained using the cubic equation

$$r_x^3 + r_2 r_x^2 + \left(r_2^2 + \frac{2}{H^2} \alpha_1\right) r_x + r_2^3 + \frac{2}{H^2} \alpha_1 r_2 + \frac{2\mu}{H^2} = 0, \quad (9)$$

where $r_x$ is the sought root of the equation $W(r) = 0$. This cubic equation can have one or three real roots. If it has one real root, it is the root $r_{-1}$. If there are three roots, they are $r_{-1}$, $r_0$, $r_1$.

Let us assume that we know two roots: $r_2$, which is given, and $r_{-1}$, which was obtained as a single real solution of the cubic equation (9). Then, representing the two complex conjugate roots of the equation $W(r) = 0$ in the form

$$a + \sqrt{-1} b, \quad a - \sqrt{-1} b,$$

where $a$ and $b$ are real, we find:

$$a = -\frac{1}{2}(r_{-1} + r_2), \quad b = \sqrt{-\frac{2\mu}{H^2(r_{-1} + r_2)} - \frac{1}{4}(r_{-1} + r_2)^2}.$$

In this case, we have the following expression for the polynomial $W(r)$:

$$W(r) = (r - r_{-1})(r - r_2)[(r - a)^2 + b^2].$$

# 4 INVERSION OF QUADRATURES AND OBTAINING ANALYTICAL SOLUTION

Before we proceed to obtain analytical solution, it is to be recalled that we exclude from consideration the cases of rectilinear motion, collision and circular motions where the roots $r_1$ and $r_2$ are equal. In these particular cases, the formulae to be deduced can turn out to be pointless. The analysis of these types of solutions may become the subject of a special investigation and is not considered here.

The processes of obtaining the analytical solution of the equations of motion and ultimate formulae differ essentially in the following three cases: finite motions, infinite motions with four real roots of the equation $W(r) = 0$, and infinite motions with two real and two complex conjugate roots of the same equation. Let us consider all these cases separately.

## 4.1 Finite motions

Finite motions can take place only if all four roots of the equation $W(r) = 0$ are real. The relationship between the variables $r$ and $u$ is expressed by the quadrature (5). We pass to the variable $u' = Hu$. To simplify our calculations without loosing general character of the solution, we assume that the arbitrary constant $u'_0$ is equal to zero. Then, the function $r(u')$ can be expressed by the quadrature

$$\int_{r_0}^{r} \frac{dr}{\sqrt{(r_2 - r)(r_1 - r)(r - r_0)(r - r_1)}} = u'$$

To inverse this integral, we use formula (254.00) from Gradshteyn & Ryzhik (1963, see p. 256). Hence, we have

$$gF(\varphi, k) = u',$$

where $F(\varphi, k)$ is an elliptic integral of the first kind, $k$ is the modulus of the elliptic integral,

$$\begin{aligned} \sin \varphi &= \sqrt{\frac{(r_1 - r_{-1})(r - r_0)}{(r_1 - r_0)(r - r_{-1})}}, \\ g &= \frac{2}{\sqrt{(r_2 - r_0)(r_1 - r_{-1})}}. \end{aligned} \quad (10)$$

The modulus of the elliptic integral $k$ is defined by the formula

$$k = \sqrt{\frac{(r_1 - r_0)(r_2 - r_{-1})}{(r_2 - r_0)(r_1 - r_{-1})}}.$$

We make the change of variables

$$u'' = \frac{u'}{g}. \quad (11)$$

Then, we have

$$F(\varphi, k) = u''.$$

Now, using Jacobi elliptic function, we get the relationship:

$$\frac{(r_1 - r_{-1})(r - r_0)}{(r_1 - r_0)(r - r_{-1})} = sn^2 u'.'$$
To simplify the calculations, we choose arbitrary constants in such a manner that if $u'' = 0$ we have $r = r_0$. Later (see below), we consider the way in which the parameters of motion are calculated in the case of arbitrary initial conditions.

From the previous formula, we obtain

$$r = \frac{r_0(r_1 - r_{-1}) - r_{-1}(r_1 - r_0)sn^2u''}{(r_1 - r_{-1}) - (r_1 - r_0)sn^2u''}.$$

Now we look for the relationship between the variable $u$ and time $t$. The variable $u$ was introduced by the relationship

$$dt = r(u) du.$$

Taking into consideration the relation between $u$ and $u''$ and substituting here the earlier found function $r(u'')$, we get

$$dt = \frac{g}{H} \frac{r_0(r_1 - r_{-1}) - r_{-1}(r_1 - r_0)sn^2u''}{(r_1 - r_{-1}) - (r_1 - r_0)sn^2u''} du''.$$

For better interpretation of the relations to be obtained, we make the change of variables

$$u'' = \frac{v}{2}$$

and transform the previous relationship in the following way:

$$dt = \frac{g}{2H} \left[ r_{-1} + (r_0 - r_{-1}) \frac{1}{1 - f \sin^2 \frac{v}{2}} \right] dv,$$

where

$$f = \frac{r_1 - r_0}{r_1 - r_{-1}}$$

satisfies the inequality $0 < f < 1$.

Subsequent transformations using the obvious equality

$$\frac{1}{1 - f \sin^2 \frac{v}{2}} = 1 + \frac{f \sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}}$$

give us

$$dt = \frac{g}{2H} \left[ r_{-1} + (r_0 - r_{-1}) \left( 1 + \frac{f \sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} \right) \right] dv,$$

and then

$$dt = \frac{gr_0}{2H} \left( 1 + \gamma \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} \right) dv, \quad (12)$$

where

$$\gamma = f \frac{r_0 - r_{-1}}{r_0} = \frac{r_1 - r_0}{r_1 - r_{-1}} \frac{r_0 - r_{-1}}{r_0}.$$

Let us consider the function

$$\Phi_1(v) = \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}}.$$

It follows from the properties of the Jacobi elliptic functions that $\Phi_1(v)$ is a periodic function with a period $4K$, where $K$ is a complete elliptic integral of the first kind with the modulus $k$. The function $\Phi_1(v)$ takes on only non-negative values. Let $P$ be the mean value of the function $\Phi_1(v)$ at the interval $(0, 4K)$, that is

$$P = \frac{1}{4K} \int_{0}^{4K} \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} dv. \quad (13)$$

It is clear that the mean value of the function $\Phi_1(v) - P$ at any interval with the length $4K$ is equal to zero.

Now we examine the integral

$$\int_{0}^{v} (\Phi_{1}(v) - P) dv$$

and represent its upper limit in the form

$$v = 4Kn + v', \quad \text{where} \quad v' = \mod(v, 4K),$$

and $n$ is an integer. Then, the following relationships are valid:

$$\int_{0}^{v} (\Phi_{1}(v) - P) dv = \int_{0}^{v'} (\Phi_{1}(v) - P) dv,$$

$$\left| \int_{0}^{v'} (\Phi_{1}(v) - P) dv \right| \le \frac{1}{1 - f} 4K.$$

Now we reduce the equation under consideration (12) to the form

$$n_v dt = \left[ 1 + \frac{\gamma}{1 + \gamma P} \left( \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} - P \right) \right] dv,$$

where

$$n_v = \frac{2H}{gr_0(1 + \gamma P)}.$$

Without loss of generality, we assume that the variable $v$ at the instance $t=0$ takes on a value of zero. Then, integrating the equation in view of the previous relationships, we obtain

$$n_v t = v + \frac{\gamma}{1 + \gamma P} \int_{0}^{v'} \left( \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} - P \right) dv.$$

Let us introduce the function

$$\Delta_v(v) = \frac{\gamma}{1 + \gamma P} \int_{0}^{\mod(v, 4K)} \left( \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} - P \right) dv. \quad (14)$$

Hence, we have

$$n_v t = v + \Delta_v(v). \quad (15)$$

Obviously, the function $\Delta_v(v)$ is periodic with the period $4K$.

The integrals appearing in (13) and (14) are expressed via elliptic functions using the relationship

$$\int_{0}^{v} \frac{\sin^2 \frac{v}{2}}{1 - f \sin^2 \frac{v}{2}} dv = \frac{2}{f} \Pi \left( am \left( \frac{v}{2}, k \right), f, k \right) - \frac{1}{f} v,$$

where $\Pi(x, f, k)$ is an incomplete elliptic integral of the third kind according to its definition:

$$\Pi(x, f, k) = \int_{0}^{x} \frac{1}{(1 - f \sin^2 t)\sqrt{1 - k^2 \sin^2 t}} dt,$$

and $am(u, k)$ is the Jacobi elliptic function. Thus, the function $\Delta_v(v)$ can be expressed in the form

$$\Delta_v(v) = \frac{\gamma}{1 + \gamma P} \left[ \frac{2}{f} \Pi \left( am \left( \frac{v'}{2}, k \right), f, k \right) - \frac{1}{f} v' - P v' \right], \quad (16)$$

and the constant $P$ can be calculated using the formula

$$P = \frac{1}{2Kf} [\Pi(\pi, f, k) - 2K].$$

We consider equation (15) as transcendental for $v$. It is clear that

$$|\Delta_v(v)| \le 4K \frac{\gamma}{1 + \gamma P} \frac{1}{1 - f}.$$

Numerical solution of the equation for $v$ can be obtained with any pre-given accuracy by means of successive approximations method using the iterative formula

$$v^{(n+1)} = n_v t - \Delta_v(v^{(n)}), \quad (n = 0, 1, 2, \dots).$$
Let us now return to the examination of the dependence of the
distance *r* on time. The formulae given above permit us to write

$$
r = \frac{r_0 - r_{-1} f s n^{2 \frac{v}{2}}}{1 - f s n^{2 \frac{v}{2}}}.
$$

Taking into account the properties of the Jacobi elliptic func-
tions, it is easy to show that the function *r*(*v*) is periodic with
the period 4*K*. In addition, we have *r* (0) = *r*<sub>0</sub>, *r* (2*K*) = *r*<sub>1</sub>. To
calculate *r*(*t*), it is necessary at first to calculate *v*(*t*), as pointed
above. Substituting this function into the obtained function *r*(*v*),
we get *r*(*t*).

Let us examine now how to calculate the longitude λ at any
instance of time. Earlier, we have obtained

$$
d\lambda = \frac{\alpha_2}{r(u)} du.
$$

In view of the changes of variables that were made and the expres-
sion for *r*(*u*′), the latter equality can be written in the form

$$
d\lambda = \frac{g\alpha_2}{2Hr_0} \left[ 1 + \gamma' \frac{sn^{2\frac{v}{2}}}{1 + f' sn^{2\frac{v}{2}}} \right] dv, \quad (17)
$$

where

$$
f' = - \frac{r_{-1}(r_1 - r_0)}{r_0(r_1 - r_{-1})}, \quad \gamma' = \frac{r_0 - r_{-1}}{r_{-1}} f'.
$$

Now we examine the function

$$
\Phi_2(v) = \frac{sn^{2\frac{v}{2}}}{1 + f' sn^{2\frac{v}{2}}}.
$$

From the properties of the Jacobi elliptic functions, it follows
that Φ₂(v) is a periodic function with the period 4K, where K
is complete elliptic integral of the first kind with the modulus
k. The function Φ₂(v) takes on only non-negative values. Let Q
be the mean value of the function Φ₂(v) at the interval (0, 4K),
that is

$$
Q = \frac{1}{4K} \int_{0}^{4K} \frac{sn^{2} \frac{v}{2}}{1 + f' sn^{2} \frac{v}{2}} dv. \quad (18)
$$

Obviously, the mean value of the function Φ₂(v) − Q at any interval
with the length 4K is equal to zero.

Without loosing the general character of our examination, we
may assume that if v = 0 the longitude λ = 0. Then, integration of
equation (17) gives us

$$
\lambda = \frac{g\alpha_2}{2Hr_0} (1 + \gamma'Q)v + \frac{g\alpha_2}{2Hr_0} \gamma' \int_0^{v'} \left( \frac{sn^{2\frac{v}{2}}}{1+f'sn^{2\frac{v}{2}}} - Q \right) dv, \quad (19)
$$

where, as earlier, $v' = \text{mod}(v, 4K)$.

The integrals appearing in the formula for *Q* (18) and in equa-
tion (19), expressing the relation of λ and *v*, are expressed via
elliptic functions by means of the relationship

$$
\int_0^v \frac{sn^{2\frac{v}{2}}}{1 + f' sn^{2\frac{v}{2}}} dv = \frac{v}{f'} - \frac{2}{f'} \Pi \left( am \left( \frac{v}{2}, k \right), -f', k \right),
$$

where Π(x, -f', k) is incomplete elliptic integral of the third kind.
The constant *Q* can be calculated using the formula

$$
Q = \frac{1}{2Kf'} \left[ \Pi(\pi, -f', k) - 2K \right], \quad (20)
$$

and the relation between $\lambda$ and $v$ takes the form

$$
\lambda = n_{\lambda} v + \Delta_{\lambda}(v), \quad (21)
$$

where

$$
n_{\lambda} = \frac{g\alpha_{2}}{2Hr_{0}}(1 + \gamma'Q),
$$

$$
\Delta_{\lambda}(v) = \frac{g\alpha_2}{2Hr_0} \gamma' \left[ \frac{v'}{f'} - \frac{2}{f'} \Pi \left( am \left( \frac{v'}{2}, k \right), -f', k \right) - Qv' \right] \quad (22)
$$

It is clear that the function $\Delta_{\lambda}(v)$ is periodic with the period 4K.

When calculating the incomplete elliptic integral in (20) and (22),
it is necessary to take into account that *f* can take on values more
than 1, while the integral's parameter can take on any negative val-
ues. In this case, one should use the formulae reducing incomplete
elliptic integrals to the parameter assuming values in the range from
$k^2$ to 1. For this we use formulae (17.7.15), (17.7.16), and (17.7.17)
from Abramovitz & Stigan (1979). To give a complete picture, we
reprint these formulae here using the designations adopted in this
paper:

$$
\Pi(\varphi, -f', k) = \left\{ \left[ (1 - N_0) \left( 1 - \frac{k^2}{N_0} \right) \right]^{1/2} \Pi(\varphi, N_0, k) \right. \\
\qquad \left. + \frac{k^2}{p_0} F(\varphi, k) \right. \\
\qquad \left. + \arctan \left[ \frac{p_0 \sin 2\varphi}{2\Delta_0(\varphi)} \right] \right\} \left[ (1 + f') \left( 1 + \frac{k^2}{f'} \right) \right]^{-1/2},
$$

where

$$
\begin{align*}
N_0 &= \frac{k^2 + f'}{1 + f'}, & p_0 &= \sqrt{\frac{f'(k^2 + f')}{1 + f'}}, \\
\varphi &= am\left(\frac{v'}{2}, k\right), & \Delta_0(\varphi) &= dn\left(\frac{v'}{2}, k\right).
\end{align*}
$$

Calculation of the longitude λ(t) is carried out in the following
way. First, for a given instance of time *t*, *v*(*t*) is calculated in the way
pointed above. Then, substituting this value into (21), we obtain the
longitude λ(t).

4.2 Properties of analytical solution for finite motions

The particle's position in a plane of motion is defined by polar
coordinates *r* (radius) and *λ* (longitude). Let us examine as to how
they depend on time.

If time varies by the value $T_v = \frac{4K}{n_v}$, the variable $v$ varies by 4K.
In its turn, $r$ assumes its previous value. Hence, it is clear that $r$ is a
periodic function of time with the period $T_v$. The longitude $\lambda$ is non-
linear in time. The particle's complete revolution around the central
body can take differing periods of time. An average rate of the
longitude's change in time is defined by an angular velocity $n_v n_\lambda$.
However, at any period of time $T_v$, the longitude $\lambda$ gets the same
increase by $n_\lambda 4K$, where $K$ is, as pointed above, complete elliptic
integral of the first kind with the modulus $k$ which is obtained from
initial conditions by the formulae given above.

From the formulae deduced in the previous subsection, it can
be concluded that the motion under consideration is conditionally
periodic with two frequencies. The first frequency, equal to 2π/Tv,
determines the rate of change of the distance r. The second one,
equal to n<sub>v</sub>n<sub>λ</sub>, determines the rate of the particle's revolutions around
the central body. If the frequencies are commensurable, that is

$$
\frac{2\pi}{4K} n_1 = n_2 n_\lambda,
$$

where $n_1$ and $n_2$ are positive integers, we get periodic motion.
## 4.3 Analytical solution for infinite motions in the case of four real roots

We use the same way of building solution as we have used for the case of finite motions. Let us take the relationship between $r$ and $u'$ in the form

$$ \int_{r_2}^{r} \frac{dr}{\sqrt{(r_2 - r)(r_1 - r)(r - r_0)(r - r_{-1})}} = u'. $$

In the case of infinite motions the variable $r$ takes on values $r > r_2$. To inverse the integral, we use formula (258.00) from Gradshteyn & Ryzhik (1963, see p. 256). Hence, we have

$$ gF(\varphi, k) = u', $$

where $F(\varphi, k)$ is an elliptic integral of the first kind and $k$ is its modulus,

$$ \sin \varphi = \sqrt{\frac{(r_1 - r_{-1})(r - r_2)}{(r_2 - r_{-1})(r - r_1)}}, \\ g = \frac{2}{\sqrt{(r_2 - r_0)(r_1 - r_{-1})}}. \qquad (23) $$

The elliptic integral's modulus $k$ is defined by the relation

$$ k = \sqrt{\frac{(r_1 - r_0)(r_2 - r_{-1})}{(r_2 - r_0)(r_1 - r_{-1})}}. $$

Note that the formulae for $g$ and $k$ are the same as in the case of finite motions.

Let us apply here the change of variables (11). We get

$$ F(\varphi, k) = u'', \qquad (24) $$

Now, using Jacobi elliptic functions, we can write

$$ \frac{(r_1 - r_{-1})(r - r_2)}{(r_2 - r_{-1})(r - r_1)} = sn^2 u''. $$

We choose arbitrary constants in the way that if $u'' = 0$ the variable $r = r_2$. From the previous formula, we find

$$ r = \frac{r_2 - r_1 s \, sn^2 u''}{1 - s \, sn^2 u''}, \qquad (25) $$

where

$$ s = \frac{r_2 - r_{-1}}{r_1 - r_{-1}} $$

and it is obvious that $s > 1$.

Now we must find the relation between the variable $u$ and time $t$. The variable $u$ was introduced by the relationship

$$ dt = r(u) du. $$

Afterwards, we also introduced the variables $u'' = \frac{u'}{g}$ and $u' = uH$. Making the change of variables and taking into account (25), we obtain

$$ dt = \frac{g}{H} \frac{r_2 - r_1 s \, sn^2 u''}{1 - s \, sn^2 u''} du''. $$

Integrating this equation, after some transformations we find

$$ t = \frac{gr_1}{H} u'' + \frac{g(r_2 - r_1)}{H} \int_0^{u''} \frac{du''}{1 - s \, sn^2 u''}. $$

From here we get

$$ n_v t = u'' + v \int_0^{u''} \frac{du''}{1 - s \, sn^2 u''}, $$

where

$$ n_v = \frac{H}{gr_1}, \quad v = \frac{r_2 - r_1}{r_1}, $$

where arbitrary constants are chosen so that if $t = 0$ the variable $u''$ takes on zero value.

The found relationship between $t$ and $u''$ contains an integral that can be expressed via incomplete elliptic integral of the third kind in the following way:

$$ \int_{0}^{u''} \frac{du''}{1 - s \, sn^2 u''} = \Pi(am(u'', k), s, k). $$

Then, we have

$$ n_v t = u'' + v \Pi(am(u'', k), s, k). \qquad (26) $$

Since the latter elliptic integral involves parameter $s$ greater than 1, for its calculation we should use the formula of reduction of incomplete elliptic integral of the third kind to the parameter with the value less than 1. For this, we use formulae (17.7.7) and (17.7.8) from Abramovitz & Stigan (1979). In our notations, we have

$$ \Pi(\varphi, s, k) = -\Pi(\varphi, N, k) + F(\varphi, k) + \frac{1}{2p_1} \ln \frac{\Delta(\varphi) + p_1 \tan \varphi}{\Delta(\varphi) - p_1 \tan \varphi}, $$

where

$$ \varphi = am(u'', k), \quad \Delta(\varphi) = dn(u'', k), $$

$$ N = \frac{k^2}{s}, \quad p_1 = \sqrt{(s-1)\left(1 - \frac{k^2}{s}\right)}. $$

It turns out that $N < 1$. Here $F(\varphi, k)$ is incomplete elliptic integral of the first kind, and $dn(u'', k)$ is Jacobi elliptic function.

Using the latter formulae, we express the relation between the variable $\varphi$ and time $t$ in the form

$$ \ln \frac{\Delta(\varphi) + p_1 \tan \varphi}{\Delta(\varphi) - p_1 \tan \varphi} = E_1(t, \varphi), $$

where

$$ E_1(t, \varphi) = 2p_1 \left[ \frac{n_v t}{\nu} - \left(1 + \frac{1}{\nu}\right) F(\varphi, k) + \Pi(\varphi, N, k) \right]. $$

From here we can deduce the relationship

$$ \varphi = \arctan \left( \frac{\Delta(\varphi)}{p_1} \frac{e^{E_1(t, \varphi)} - 1}{e^{E_1(t, \varphi)} + 1} \right), $$

which can be used to calculate $\varphi$ at a given instance $t$ by the method of successive approximations. However, in this case convergence of successive approximations cannot be ensured. Convergence is provided if we reduce the latter relationship to the form

$$ \varphi = \arctan \left[ \frac{1}{1 + \frac{\Delta(\varphi)}{\nu}} \left( \frac{\Delta(\varphi)}{p_1} \frac{e^{E_1(t, \varphi)} - 1}{e^{E_1(t, \varphi)} + 1} + \frac{\Delta(\varphi)}{\nu} \right) \right]. $$

Taking in the right-hand side $\varphi = 0$ as zero-order approximation, we obtain the first-order approximation $\varphi$. Substituting it into the right-hand side of the latter relationship again, we get next order of approximation and so on. After $\varphi$ is obtained, we find $u''$ from the formula

$$ u'' = F(\varphi, k). $$

Having got relation between $t$ and $u''$, the function $r(t)$ can be obtained by substituting $u''(t)$ into (25).
Now let us find the dependence of the longitude $\lambda$ on time $t$. From the general form of integral, we have

$$d\lambda = \frac{\alpha_2}{r} du.$$

Substituting the expression (25) for $r(u'')$ and changing $u$ to $\frac{g}{H}u''$, we obtain

$$d\lambda = \frac{g\alpha_2}{H} \frac{1 - s \, sn^2 u''}{r_2 - r_1 \, s \, sn^2 u''} du''$$

or, after some simple transformations,

$$d\lambda = \frac{g\alpha_2}{Hr_1} du'' - \frac{g\alpha_2}{Hr_1} \frac{r_2 - r_1}{r_2} \frac{du''}{1 - s' \, sn^2 u''} du''$$

where

$$s' = \frac{r_1}{r_2}.$$

It is easy to prove that $0 < s' \le 1$.

Integrating the latter equation in assumption that $\lambda = 0$ when $v = 0$, we obtain

$$\lambda = n_l u'' - n_l \frac{r_2 - r_1}{r_2} \int_0^{u''} \frac{du''}{1 - s' \, sn^2 u''} du''$$

where

$$n_l = \frac{g\alpha_2}{Hr_1}.$$

Expressing the integral contained herein via incomplete elliptic integral of the third kind, we find

$$\lambda = n_l u'' - n_l \frac{r_2 - r_1}{r_2} \Pi(am(u''), s', k).$$

Substituting the earlier found function $u''(t)$, we obtain the desired dependence of $\lambda$ on time $t$.

## 4.4 Analytical solution for infinite motions in the case of two real and two complex conjugate roots

Let us carry out the inversion of quadratures in the case when the equation

$$W(r) = 0$$

has two real and two complex conjugate roots. We denote the real roots as $r_{-1}$ and $r_2$. It is clear that $r_{-1} < 0$ and $r_2 > 0$. We have already denoted the complex conjugate roots earlier as $a + \sqrt{-1} b$ and $a - \sqrt{-1} b$, where $a$ and $b$ are real numbers. Then, we can write the quadrature relating the variables $r$ and $u'$ in the form

$$\int_{r_2}^{r} \frac{dr}{\sqrt{(r - r_{-1})(r - r_2)[(r - a)^2 + b^2]}} = u'$$

Here, as earlier, arbitrary constants are chosen in the way that if $u' = 0$ the variable $r$ takes on the value $r_2$.

In the case of infinite motions the variable $r$ assumes values satisfying the inequality $r > r_2$. Hence, to inverse the integral, we use formula (260.0) from Byrd & Friedman (1954). We have

$$g_1 F(\varphi, k_1) = u',$$

where

$$\cos \varphi = \frac{A(r - r_{-1}) - B(r - r_2)}{A(r - r_{-1}) + B(r - r_2)},$$

$$g_1 = \frac{1}{\sqrt{AB}},$$

$$k_1 = \sqrt{\frac{(A + B)^2 - (r_2 - r_{-1})^2}{4AB}}. \qquad (27)$$

In these formulae, we use notations:

$$A = \sqrt{(r_2 - a)^2 + b^2}, \quad B = \sqrt{(r_{-1} - a)^2 + b^2}.$$

$F(\varphi, k_1)$ is incomplete elliptic integral of the first kind, so we have

$$F(\varphi, k_1) = \frac{H}{g_1} u = u'', \quad \cos \varphi = cn(u'', k_1). \qquad (28)$$

Now we can write the function $r(u'')$ in the following way:

$$r = \frac{Br_2 - Ar_{-1} + (Ar_{-1} + Br_2)cnu''}{B - A + (A + B)cnu''}.$$

The variable $u$ was introduced by the differential relationship $dt = r \, du$. Changing over to the variable $u''$ and carrying out integration in assumption that $u'' = 0$ if $t = 0$, we find

$$t = \frac{g_1}{H} \int_{0}^{u''} r(u'') du''.$$

We substitute the earlier found function $r(u'')$ into this relationship and, after some transformations, obtain:

$$n'_{v}t = u'' + \eta \int_{0}^{u''} \frac{1}{1 + \beta \, cn \, u''} du'', \qquad (29)$$

where

$$n'_{v} = \frac{H}{g_1} \frac{A + B}{Br_2 + Ar_{-1}}, \quad \beta = \frac{A + B}{B - A},$$

$$\eta = \frac{2AB(r_2 - r_{-1})}{(Br_2 + Ar_{-1})(B - A)}.$$

Here, it turns out that $\beta > 1$.

Calculation of the integral in (29) is carried out using the formula

$$\int_{0}^{u''} \frac{1}{1 + \beta \, cn \, u''} du'' = - \frac{\Pi(\varphi, p, k_1)}{\beta^2 - 1} + \frac{1}{2} \frac{\beta}{\sqrt{\beta^2 - 1}} \ln \frac{\sqrt{\beta^2 - 1} \Delta_1(\varphi) + m_1 \sin \varphi}{\sqrt{\beta^2 - 1} \Delta_1(\varphi) - m_1 \sin \varphi}$$

where

$$p = \frac{\beta^2}{\beta^2 - 1}, \quad m_1 = \sqrt{k_1^2 + \beta^2 - k_1^2 \beta^2},$$

$$u'' = F(\varphi, k_1), \quad \Delta_1(\varphi) = dn(u'', k_1).$$

Since the incomplete elliptic integral of the third kind contained herein involves parameter $p$ greater than 1, we should reduce it to the parameter $N_1 < 1$ using formulae (17.7.7) and (17.7.8) from
Abramovitz & Stigan (1979) in the following way:

$$ \Pi(\varphi, p, k_1) = -\Pi(\varphi, N_1, k_1) + F(\varphi, k_1) + \frac{1}{2p'_1} \ln \frac{\Delta_1(\varphi) + p'_1 \tan \varphi}{\Delta_1(\varphi) - p'_1 \tan \varphi}, \quad (29) $$

where

$$ N_1 = \frac{k_1^2}{p}, \quad p'_1 = \sqrt{(p-1)\left(1 - \frac{k_1^2}{p}\right)}. $$

In addition, it turns out that

$$ p'_{1} = \frac{\sqrt{k_{1}^{2} + \beta^{2} - k_{1}^{2}\beta^{2}}}{\beta\sqrt{\beta^{2} - 1}} = \frac{m_{1}}{\beta\sqrt{\beta^{2} - 1}}. $$

Now the relationship between the variable $\varphi$ and time $t$ can be represented in the form

$$ \frac{(\Delta(\varphi) - p'_1 \tan \varphi)(\Delta_1(\varphi) + p_2 \sin \varphi)}{(\Delta(\varphi) + p'_1 \tan \varphi)(\Delta_1(\varphi) - p_2 \sin \varphi)} = e^{E_2}, \quad (30) $$

where

$$ E_2 = 2p'_1(\beta^2 - 1) \left[ \frac{n'_v t}{\eta} + \left( \frac{1}{\beta^2 - 1} - \frac{1}{\eta} \right) F(\varphi, k_1) - \frac{1}{\beta^2 - 1} \Pi(\varphi, N_1, k_1) \right]. \quad (31) $$

Here, we have introduced the notation

$$ p_2 = \frac{m_1}{\sqrt{\beta^2 - 1}} $$

and used the relationships

$$ u'' = F(\varphi, k_1), \quad \sin \varphi = sn u''. $$

For a given time $t$, we have obtained transcendental equation for $\varphi$. It can be solved numerically by the successive approximations method. However, convergence of iterations can be assured only if the above relationships have specific forms. Thus, when $\varphi$ assumes values satisfying the condition

$$ \varphi \le \arcsin \frac{\Delta_1(\varphi)}{p_2}, $$

to perform iterations, we should use the relationship

$$ \varphi = \arcsin \left[ \frac{\Delta_1(\varphi) C_2 e^{E_2} - 1}{p_2 C_2 e^{E_2} + 1} \right], $$

where $C_2$ is calculated by the formula

$$ C_2 = \frac{\Delta_1(\varphi) + p'_1 \tan \varphi}{\Delta_1(\varphi) - p'_1 \tan \varphi}. $$

In the case when

$$ \varphi > \arcsin \frac{\Delta_1(\varphi)}{p_2}, $$

we should use the relationship

$$ \varphi = \arctan \left[ \frac{\Delta_1(\varphi) \left(1 - C_3 e^{E_2}\right)}{p'_1 \left(1 + C_3 e^{E_2}\right)} \right], $$

where $C_3$ is calculated by the formula

$$ C_3 = \frac{\Delta_1(\varphi) - p_2 \sin \varphi}{\Delta_1(\varphi) + p_2 \sin \varphi}. $$

Since, for an arbitrarily given $t$, we do not know in advance which values of $\varphi$ we will get during iterations, we should at first assume

$\varphi = 0$ and start iterations using the first formula. If, performing iterations, we get negative value for $C_2$, we should set the initial value for $\varphi$ to be equal to $\pi/2$ and resume iterations using the second formula.

Our calculations show that convergence of iterations in the proposed methods for numerically solving the considered transcendental equations can be very slow. For some values of initial parameters, these iterations do not converge at all. However, what appears to be more effective and reliable is dichotomy method (method of division in halves) of solving the transcendental equations. Let us consider as to how this method can be applied in our problem.

The equation for $\varphi$ at a given instance of time should be written in the form

$$ G(\varphi) = e^{E_2} - \frac{(\Delta_1(\varphi) - p'_1 \tan \varphi)(\Delta_1(\varphi) + p_2 \sin \varphi)}{(\Delta_1(\varphi) + p'_1 \tan \varphi)(\Delta_1(\varphi) - p_2 \sin \varphi)} = 0. $$

The sole root of this equation is to be searched for at the interval $(0, \bar{\varphi})$, where $\bar{\varphi}$ is to be calculated using the relationships

$$ \sin \bar{\varphi} = \frac{1 + p'_1^2 + k_1^2 - \sqrt{(1 + p'_1^2 + k_1^2)^2 - 4k_1^2}}{2k_1^2}, $$

$$ \cos \bar{\varphi} = -\sqrt{1 - \sin^2 \bar{\varphi}}. $$

If $t > 0$, we have $G(0) > 0$, and if $\varphi = \bar{\varphi}$, the function $G(\varphi)$ turns into infinity. However, this does not lead to insuperable problem of calculations because the dichotomy method should be used to calculate the value of the function at the middle of the interval. For the values of $\varphi$ less than but sufficiently close to $\bar{\varphi}$, the value of $G(\varphi)$ is considerably negative. The number of necessary iterations in the dichotomy method can be easily evaluated by the formula

$$ N_{\text{iter}} = -\log_2 \varepsilon, $$

where $\varepsilon$ is a preset accuracy with which the equation's root is defined. For instance, if accuracy is $10^{-10}$, we need 34 iterations.

Now we proceed to derive the dependence of the second coordinate, the longitude $\lambda$, on $\varphi$. We start out with the differential relationship

$$ d\lambda = \frac{\alpha_2}{r} du = \frac{\alpha_2 g_1}{rH} du''. $$

We substitute here the earlier obtained relation $r(u'')$ and integrate the above relationship assuming that $\lambda = 0$ if $u'' = 0$. After some simple transformations, we obtain

$$ \lambda = n'_{1} \left( u'' + q \int_{0}^{u''} \frac{u''}{1 + \beta' c n''} du'' \right), $$

where

$$ n'_{1} = \frac{\alpha_2 g}{H} \frac{A + B}{Br_2 + Ar_{-1}}, \quad q = \frac{2AB(r_{-1} - r_2)}{(A + B)(Br_2 - Ar_{-1})}, $$

$$ \beta' = \frac{Br_2 + Ar_{-1}}{Br_2 - Ar_{-1}}. $$

The integral contained herein is expressed via elliptic functions. In result, we get the final expression for the coordinate $\lambda$ in the following form:

$$ \lambda = n'_{1} F(\varphi, k_1) + q \left[ \frac{1}{1 - \beta'^2} \Pi(\varphi, p', k_1) - \frac{\beta'}{2m_2 \sqrt{1 - \beta'^2}} \ln \frac{\Delta_1(\varphi) + p'_2 \sin \varphi}{\Delta_1(\varphi) - p'_2 \sin \varphi} \right], $$
where

$$p' = \frac{\beta'^2}{\beta'^2 - 1}, \quad m_2 = \sqrt{\beta'^2 + k_1^2 - \beta'^2 k_1^2},$$

$$p_2' = \frac{m_2}{\sqrt{1 - \beta'^2}}.$$

## 4.5 Calculation of the particle's coordinates in the case of arbitrary initial conditions

So far we chose arbitrary constants in a way that was convenient to carry out calculations and obtain the solutions. Let us now consider the way in which the parameters of motion are calculated if we have arbitrary initial conditions.

Suppose that for the initial moment of time $t_0 = 0$, the values of the following parameters are given: initial distance $r^{(0)}$, initial radial velocity $\dot{r}^{(0)}$, initial longitude $\lambda^{(0)}$ and initial angular velocity $\dot{\lambda}^{(0)}$.

First, we calculate the initial velocity

$$V^{(0)} = \sqrt{(\dot{r}^{(0)})^2 + (r^{(0)} \dot{\lambda}^{(0)})^2},$$

and, after that, the constants $\alpha_1$ and $\alpha_2$ are calculated using the formulae

$$\alpha_1 = \frac{1}{2}(V^{(0)})^2 - \frac{\mu}{r^{(0)}} - \frac{H^2}{2}(r^{(0)})^2, \quad \alpha_2 = (r^{(0)})^2 \dot{\lambda}^{(0)}.$$

Then, we should find the only negative root $r_{-1}$ of the equation:

$$r^4 + 2 \frac{\alpha_1}{H^2} r^2 + 2 \frac{\mu}{H^2} r - \frac{\alpha_2^2}{H^2} = 0.$$

This can be done using the dichotomy method. After that, we find the solution of the cubic equation:

$$r_x^3 + r_{-1} r_x^2 + (r_{-1}^2 + \frac{2}{H^2} \alpha_1) r_x + r_{-1}^3 + \frac{2}{H^2} \alpha_1 r_{-1} + \frac{2\mu}{H^2} = 0, \quad (32)$$

where $r_x$ is the unknown. If we obtain three real roots $r_0, r_1, r_2$, it is necessary to check if $r^{(0)} \le r_1$.

If the latter condition is satisfied, we have the case of finite motion. Then, we calculate $\varphi$ using the relationship (10) where we put $r = r^{(0)}$. Note that the value of $\varphi$ should be taken in the interval $(0, \pi/2)$ if $\dot{r}^{(0)} > 0$, and in the interval $(\pi/2, \pi)$ if $\dot{r}^{(0)} < 0$. Having calculated necessary constants, we obtain the value of $v'$ from the relationship

$$v' = 2F(\varphi, k),$$

where $F(\varphi, k)$ is incomplete elliptic integral of the first kind. Next, we calculate $\Delta_v(v)$ using (16) and obtain the moment $t_0$ from relationship (15) where we put $v' = v, t = t_0$. Now, to calculate the particle's coordinates for any given moment of time $t$, one can use the above formulae beginning from relationship (15) where $t$ should be replaced by $t + t_0$.

Otherwise, if the condition $r^{(0)} \le r_1$ is not satisfied, the inequality $r^{(0)} \ge r_2$ is true, and we have the case of infinite motion with for real roots $r_{-1}, r_0, r_1, r_2$. To define the coordinates in this case, we first find $\varphi$ from relationship (23), where $r$ is replaced by $r^{(0)}$, and $u''$ using (24). Then $t_0$ is calculated from relationship (26) with $t = t_0$. Hence, for any given moment of time $t$, the coordinates can be calculated using the above formulae beginning from relationship (26) where $t$ should be replaced by $t + t_0$.

If, solving equation (32), only one real root $r_2$ is obtained, we have the case of infinite motion. In this case, the complex roots

should be defined in such a manner as it was specified in Section 3.2. Then, calculations are carried out in the way similar to that in the previous case, using formulae (27), (28), (30) and (31).

## 5 EXAMPLES OF TRAJECTORIES

To perform calculations by the formulae of analytical solution of the equations of motion, a calculating program in c programming language was written. The program works without using additional subroutine libraries and has been compiled by Visual Studio 2008.

To demonstrate the properties of the trajectory characteristic of this problem and the possibility to do calculations by using the deduced formulae, we offer a number of examples. In these examples, the units for measuring distance and time were chosen so that $\mu = 1$ and $H = 1$. Then, it turns out that $\bar{r} = 1$ and $r^* = 0.629960524947$.

It is quite natural in our problem to divide all possible trajectories into typical families. First, we consider the family of the finite trajectories where the motions are such that the particle revolves around the attracting centre. In these motions, the distance reaches periodically its maximum and minimum values, and it is at these instances when velocity vector is perpendicular to the position vector. In the given examples, all trajectories begin in the point where the distance is minimum. The trajectories of this family are shown in Figs 4–6. In Fig. 6, we see the commensurability between the period of the longitude and that of the distance's variations between the minimum and maximum values.

The next three families are segments of the infinite trajectories. In all cases, the motion begins at the point where the distance has its minimum value $r_2$ and initial velocity vector is perpendicular to the position vector. In the first family of the infinite trajectories, initial point has minimum distance $r_2 = 0.5 < r^*$. In our examples, the particle starts with differing values of initial velocities. These examples are shown in Fig. 7. In the second family, shown in Fig. 8, the minimum distance satisfies the inequality $r^* < r_2 = 0.9 < \bar{r}$.

<figure><img src="image_5.png" /><figcaption>Figure 4. Example of a trajectory of finite motion. The value of $\bar{r}$ was taken as a unit of distance. The distance varies in the range 0.1–0.6. The dotted line marks the circle with the radius $r^*$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Example of a trajectory of finite motion. The value of $\bar{r}$ was taken as a unit of distance. The distance varies in the range 0.01–0.8. The dotted line marks the circle with the radius $r^*$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Example of a trajectory of finite motion in the case of the com-
mensurability between the frequency of the longitude and that of the dis-
tance's variation: when the particle makes six revolutions around central
body, it reaches minimum distance five times. The value of $\bar{r}$ was taken as
a unit of distance. The dotted line marks the circle with the radius $r^*$.</figcaption></figure>

In the third family, the particle starts with the minimum distance $r_2 = 1.2 > \bar{r}$. These examples are shown in Fig. 9.

The calculation program and files containing examples of calcu-
lations may be obtained via e-mail. For this, send your request to
the first author at the address emelia@sai.msu.ru.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Examples of trajectories of infinite motion. The value of $\bar{r}$ was taken as a unit of distance. The time unit was chosen so that $\mu = 1$ and $H = 1$. In each instance, the particle is starting off from the x-axis at a distance 0.5 from the origin of coordinates, perpendicular to the x-axis. Initial velocities for each instance are given near the plots. The circle corresponds to the distance $r = 1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Examples of trajectories of infinite motion. The value of $\bar{r}$ was taken as a unit of distance. The time unit was chosen so that $\mu = 1$ and $H = 1$. In each instance, the particle is starting off from the x-axis at a distance 0.9 from the origin of coordinates, perpendicular to the x-axis. Initial velocities for these instances were taken as follows: 0.64, 0.8, 1.36, 2.20. The circle corresponds to the distance $r = 1$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Examples of trajectories of infinite motion. The value of $\bar{r}$ was taken as a unit of distance. The time unit was chosen so that $\mu = 1$ and $H = 1$. In each instance, the particle is starting off from the x-axis at a distance 1.2 from the origin of coordinates, perpendicular to the x-axis. Initial velocities for each instance are given near the plots. The circle corresponds to the distance $r = 1$.</figcaption></figure>

## 6 NUMERICAL INTEGRATION OF THE EQUATIONS OF MOTION AND TESTING THE VALIDITY OF THE ANALYTICAL SOLUTION

Equations (1) can be solved by numerical integration methods. This can be done to get particular solutions with different initial conditions if the time interval is not big and if we do not need high accuracy of the solution. We have carried out numerical integration of the equations to test the obtained analytical solution. Among available methods and corresponding calculating programs, we have chosen the method of Everhart (1974), which is widely and successfully applied in celestial mechanics.

First attempts to integrate equations (1), straightforwardly resulted in obviously wrong solutions in the cases when the particle passes close to the zero acceleration circle. These failed attempts can be explained by the analysis of the right-hand sides of equations (1): the right-hand sides contain the factor $H^2 - \mu/r^3$ which is equal to zero at the zero acceleration circle. At every step of the numerical integration, the right-hand sides are calculated with the given value of $r$. The calculations are always performed with limited accuracy with which numbers are stored in computer memory. Any number is represented in exponential form as mantissa and exponent. Suppose that, for any number, computer memory stores $M$ decimal places of mantissa. Close to the zero acceleration circle, calculation of the factor $H^2 - \mu/r^3$ is done by subtraction of two close positive numbers, and if their first $M_1$ decimal places coincide, the subtraction gives the number for which only $M - M_1$ decimal places are correct. Thus, accuracy of calculations sharply decreases. Since coordinates obtained at every previous step of integration are used at following steps, calculation errors are accumulated. As a result, after the particle passes close to the zero acceleration circle, the solution has significant errors.

To overcome this difficulty, we should express the equations of motion in other variables. We choose as sought-for functions the variables $\xi$, $\lambda$ related to rectangular coordinates in the following way:

$$x = r \cos \lambda, \quad y = r \sin \lambda, \quad \xi = r - \bar{r}.$$

We also use the first integral of the equations of motion:

$$x \frac{dy}{dt} - y \frac{dx}{dt} = r^2 \frac{d\lambda}{dt} = M,$$

where arbitrary constant $M$ can be calculated from initial conditions. Then, the equations of motions expressed in the new variables are as follows:

$$\frac{d\lambda}{dt} = \frac{M}{(\bar{r} + \xi)^2},$$

$$\frac{d^2\xi}{dt^2} = \xi H^2 \left[ 1 + \frac{\bar{r}}{\bar{r} + \xi} + \left( \frac{\bar{r}}{\bar{r} + \xi} \right)^2 \right] + \frac{M^2}{(\bar{r} + \xi)^3}.$$

Numerical integration of these equations does not result in the loss of accuracy when motion takes place close to the zero acceleration circle.

In our examples of calculations initial conditions were defined by the parameters $r_0$, $V_0$ in the case of finite motions and by the parameters $r_2$, $V_2$ in that of infinite motions. Hence, for the variables $\xi$, $\lambda$, initial conditions are calculated by the formulae

$$\xi_0 = r_0 - \bar{r}, \quad \frac{d\xi_0}{dt} = 0, \quad \lambda_0 = 0, \quad M = (\xi_0 + \bar{r})V_0$$

in the case of finite motions. For infinite motions, the formulae are similar.

For every example of calculations by the formulae of analytical solution, comparison has been made with the results of numerical integration. The comparison has shown that each time the results coincide with the accuracy $10^{-9}$. The analysis of the calculations has shown that the accuracy is mainly restricted by the program for calculations of incomplete elliptic integrals of the third kind.

## 7 CONCLUSION

In this paper, we have examined the properties of motions in the two-body problem against the cosmic vacuum background. It is shown that finite motions take place only within a circle with the radius $\bar{r} = (\mu/H^2)^{1/3}$ centred in the attracting body. The minimum distance in finite motions cannot exceed the value $r^* = (\mu/4H^2)^{1/3}$. Circular motions can take place in the area of finite motions. The circular motions with the radii in the range $(r^*, \bar{r})$ are orbitally unstable. This means that slightest variations in the initial conditions turn circular motion either into finite motion with a great amplitude of distance's variation or into infinite motion. There are also infinite motions which can have any value of the minimum distance.

The principal result of our examination is the precise analytical solution of the two-body problem against the cosmic vacuum background. We have deduced formulae for the calculation of a moving point's coordinates at any moment of time for given initial conditions. First, we reduced the solution to quadratures, then the quadratures were inverted using Jacobi elliptic functions. The resulting formulae contain incomplete elliptic integrals of the first and third kind. The chain of the calculation formulae includes a transcendental equation analogous to that of Kepler in the Keplerian motion but containing elliptic integrals. This equation can
be numerically solved with any pre-given accuracy. The calcula-
tion formulae are programmed using c programming language. A
number of examples of calculations have been made, and the corre-
sponding plots are shown in the figures. The calculation programs
and files containing examples of calculations may be sent by the
first author of the paper at request.

We have also composed a program for numerical integration of
the equations of motion. It is shown that numerical integration of
the equations of motion expressed in rectangular coordinates can
result in sharp loss of accuracy when the point passes close to the
circle with the radius $\bar{r}$. However, a change in variables has been
found leading to the equations which, when numerically integrated,
give no loss of accuracy. By comparing the calculation results by
the formulae of the analytical solution with those of numerical
integration, we have proven that the deduced formulae are correct.

In order to save space, when deducing the formulae, we have ex-
cluded from consideration the case of rectilinear motions as well as
some cases of circular motions, because they would require special
attention and analysis.

ACKNOWLEDGMENT

This work was supported by the Russian Foundation for Basic
Research, project no. 13-02-00137.

REFERENCES

Abramovitz M., Stigan I., 1979, Spravochnik po specialnym funkciyam,
Nauka, Moscow

Byrd P. F., Friedman M. D., 1954, Handbook of Elliptic Integrals for Engineers and Physicists, Springer, Berlin

Chernin A. D., 2001, Physics Uspekhi, 44, 1099

de Sitter W., 1917, MNRAS, 78, 3

Dolgachev V. P., Domozhilova L. M., Chernin A. D., 2003, Astron. Rep., 47, 728

Everhart E., 1974, Celest. Mech., 10, 35

Friedmann A. A., 1922, Z. Phys., 10, 377

Friedmann A. A., 1924, Z. Phys., 21, 326

Gradshteyn I. S., Ryzhik I. M., 1963, Tablicy integralov, summ, ryadov i proizvedeniy, Gosudarstvennoe izdatelstvo fiz.-mat. literatury, Moscow

Lyk'yanov L. G., 2009, Astron. Rep., 53, 9

Nandra R. A., Lasenby A. N., Hobson M. P., 2012, MNRAS, 422, 2945

Perlmutter S. et al., 1999, ApJ, 517, 565

Poleshchikov S. M., Zhubr A. V., 2008, Cosm. Res., 46, 202

Riess A. G. et al., 1998, AJ, 116, 1009

This paper has been typeset from a TeX/LA TeX file prepared by the author.