Gravitational Interaction of Objects Moving in Intersecting Orbits

(Short title: Gravitational Interaction of Objects)

S. I. Ipatov

Keldysh Institute of Applied Mathematics, Russian Academy of Sciences,
Moscow, Russia

@Abstract - The evolution of two, three, and one hundred gravitating
objects, i.e., material points moving around the Sun in crossing orbits, is
investigated numerically, mainly by the action spheres method.
It is demonstrated that in the
case of three identical objects, maximum eccentricities can exceed those for
two like objects by a factor of several tens. Eccentricities of
heliocentric orbits of objects moving in a sphere of action are examined to
reveal the dependence of their mean variations on the starting data, i.e.,
masses, eccentricities, and semimajor axes. The results of these
examinations and analytical evaluations provide a basis for treating cases,
where, in the course of the evolution of a disk comprising many small objects
and that comprising a lesser number of bigger objects, the increments of average
eccentricities are the same. The obtaind results suggest that, due to their 
gravitational interaction, the possibility of
migration of several bodies from  to Neptune's from the belt beyond the Neptune exists.

Received June 27, 1994

INTRODUCTION

Recently I have investigated the evolution of orbits of two gravitationally
interacting objects that are moving around the Sun in close orbits [19, 20].
Most of the results were obtained by numerical integration of the equations
of motion for the planar three-body problem. The time interval covered was
as much as 25 000 revolutions of the objects around the Sun.

In this paper I present the outcomes of my investigations concerning the
interaction of two, three, or several hundreds of objects that move around
the Sun in intersecting orbits. It is the employment of the action
spheres method, in which we are dealing with two two-body problems, that
made it possible to examine the evolution of orbits for many objects on
large time intervals. The results obtained by the use of this method are
confronted with those derived from numerical integration. Findings gained in
our computations and analytical evaluations of average eccentricity variations
in an evolving disk allow better insight into the evolution of planetesimal
orbits in a protoplanet disk as well as into the evolution of orbits for
bodies in the asteroidal belt and the belt beyond Neptune. In what follows,
we shall look upon the objects mainly as being mass points. @Bodies@, which
congregate under collisions, and the evolution of disks of such bodies were
discussed in the literature reviewed in [17].

TWO GRAVITATIONALLY INTERACTING OBJECTS

The evolution of orbits of two gravitationally interacting objects moving
around the Sun was studied first. The results obtained by numerical
integration of the equations of motion in the framework of planar
three-body problem were described in [19, 20]. Ranges of magnitude of
[epsilon]@0 = (@a@2@0 - a@1@0)/a@1@0,
corresponding to different types of orbital element variations,
were examined for several values of [mu]@1, [mu]@2, and [phi]@0. Here
[mu]@1 and [mu]@2 are the ratios of the object masses to the mass of the Sun,
[phi]@0 is the initial angle between directions to the first and second
object, with the apex at the Sun, and @a@1@0, @a@2@0 are the semimajor axes
of the initial orbits of objects. Maximum values of [epsilon]@0 were
determined: [beta] = max [epsilon]@0 for the motions around triangular points
of libration, [gamma] = max [epsilon]@0 for the motions that allow close
encounters of objects, and [delta] = max [epsilon]@0 for the case of chaotic
variation of orbital elements; [beta] < [gamma]@[delta].

For large time intervals the method of spheres is used extensively in
studies of orbit evolution of gravitating objects. In this approach, it
is assumed that outside of the sphere the objects orbit the Sun in
unperturbed Keplerian orbits. Inside the sphere relative motion of two
objects is treated in the framework of a two-body problem, whereas their
center of mass orbits the Sun in an unperturbed Keplerian orbit. In
general use is the sphere of action, in foreign literature often referred
to as the Tisserand sphere [31]. The radius of this sphere is
@r@[sa] = @R[mu]@[2/5], where [mu] = [mu]@1 + [mu]@2
and @R is the solar distance of objects. The algorithm for
applying the method of spheres used by me is described briefly in [3]
and [15]. A more extensive description of this algorithm was presented in
the Report of Keldysh Institute of Applied Mathematics of the Academy of
Sciences of the USSR, 1985, no. 1211.

The actual motion of objects inside of the sphere of action is poorly
approximated by the method of action spheres, when relative velocities
@v@r of the objects approaching each other are small [27, 40]. Moreover,
even with rather close approximation of the motion inside the sphere, the
velocities of objects leaving the sphere of action may differ
considerably from the actual ones [30]. However, in tentative studies of
the orbit evolution over the course of many encounters one can employ the
method of spheres, even though the values of @v@r are small. I have
compared several tens of plots representing the time variations of
of heliocentric orbital elements for two gravitating objects, which were
obtained by the method of spheres, with similar plots get from numerical
integration of the equations of motion [4, 7]. The examination of these plots
demonstrated that the method of spheres allows the limits and major
trends of heliocentric orbital element variations to be found. In doing so
the sphere of action is to be used in the case @e@M > 2[mu]@[1/3],
where @e@M is the
major eccentricity of two heliocentric orbits of the objects approaching
each other. Otherwise, a sphere with the radius
@r@s = R[delta]@1.5R[mu]@[2/7] or @r@s = R[gamma]@2.4R[mu]@[1/3]
is to be used, as the problem in hand requires.

One should take @r@s = R[gamma], if one needs to avoid the "superfluous"
collisions. When @r@s = R[delta] is used in the case @[epsilon]@,
the encounters involving change of the sign of the quantity
[epsilon] = (@a@2 - a@1)/a@1 are
to be discarded and simulated once more with other pseudorandom numbers [3].
The subscript @b above labels the values of [epsilon] and
@e@ = max(@e@1, e@2) before the encounter;
here @e@1 and @e@2 are the eccentricities of the object orbits.

For objects moving in almost circular orbits one has to take into account
that for @[epsilon]@0@ < [beta] and @[Delta][phi]@0@ > 4[mu]@[1/3]
close encounters are impossible [19, 20].
Gravitational influence of the Jupiter on the evolution of
resonant asteroidal orbits can result in a strong increase of their
eccentricities [14, 16, 42, 43]. The method of spheres gives no way to
account for such effects.

I have used the method of action spheres to examine the evolution of
orbits of two gravitationally interacting mass points that circle the Sun
in the same plane. For a number of starting data the results were
presented in [3]. In the series of computations being considered
starting values were chosen as follows:
10@[-9]@[mu]@1@10@[-5], [mu]@1@[mu]@2, [epsilon]@0 = 1.01,
0.05@e@0 = e@1@0 = e@2@0@0.15, [pi]@1@0 = 1 rad, and [pi]@2@0 = 4 rad,
where e@k@0 and [pi]@k@0 are starting values for the eccentricity @e
and the perihelion longitude [pi] of the @kth object's orbit, respectively.

For the case [mu]@1 = [mu]@2, as long as sufficiently large time interval
of evolution is considered, maximum eccentricities @e@max of object orbits
depend mainly on their initial eccentricities rather then on the masses of
objects, with @e@[max]@2e@0 and [delta][pi] = [pi]@2 - [pi]@1@180@[degree],
whereas the aphelion of one
orbit is close to the perihelion of the other. Thus, for instance, with
@e@0 = 0.05 and [mu]@1 = [mu]@2 = 10@[-9] the eccentricities increased
up to @e@[max]@2e@0
in the course of 700 mutual approaches of the objects up to the sphere of
action radius @r@[sa]; they occured in 10@7 revolutions of objects
around the Sun.
The plots of time variations of semimajor axes @a@1 and @a@2 were "symmetric"
with respect to 2@a@1@0a@2@0/(a@1@0 + a@2@0),
so that the objects interchanged periodically the
values of their semimajor axes. Such an "exchange" of orbits is
attributable to the fact that in some encounters the velocity @v@r of the
second object relative to the center of mass is small, and the direction of
the vector
@v@r at exit is almost opposite to that at entry into the action sphere,
whereas the absolute value of @v@r changes scarcely. The "symmetry" of
variations of @a@1 and @a@2 follows from the energy integral. In the case of
[mu]@1 = [mu]@2 the time variation plots were almost identical for both @e@1
and @e@2.

In the case of [mu]@1@[mu]@2 the variation of orbital elememts of the second
object was greater then that for [mu]@1 = [mu]@2. Thus, for example, with
@e@0 = 0.15 and [mu]@1 = 10@[-7]@[mu]@2 the eccentricity of the second
object's orbit
increased up to 0.5 in the course of 10@3 approaches up to @r@[sa], which
occured in 10@6 revolutions around the Sun, whereas [delta][pi] took on
different values in the course of evolution.

Numerical integration of three-body problem show that the semimajor axis of
an asteroidal orbit, which crosses the orbit of Ceres with the mean angle
between the planes of orbits [delta]@i = 10@[degree], could be altered by
[Delta]a@0.01 a.u. during the
lifetime of the Solar System [3]. In the planar model the rate of change of
@a is three orders of magnitude higher. It is believed that the mass of
Ceres may exceed the overall mass of other asteroids. The above results
allow the possibility that several percent of the asteroids from the main
belt entered the Kirkwood gaps during the lifetime of the Solar System.
Some of these asteroids could leave the asteroidal belt under the influence
of planets [14, 16, 42, 43].

INTERACTION OF THREE OR MORE OBJECTS

For three gravitationally interacting objects moving in close orbits
maximum eccentricities can be several tens times greater than those for two
objects. By way of example assume that two gravitationally interacting
objects have masses close to that of Pluto and circular initial orbits.
Maximum eccentricities of their orbits generally appear to be less than
0.01. For three such objects with initially circular orbits I investigated
the evolution of their orbits in the course of 4@10@5 revolutions, i.e.,
during @10@8 years at a@50 a.u. For this purpose, the equations of motion
of planar four-body problem were considered. Numerical treatment of these
equations was accomplished using an extrapolation algorithm with variable
step size, namely, the routine BULSTO due to Bulirsch and Stoer [25].
Computed orbital elements @a, @e, and [pi] versus @N@r, the number of
revolutions of the first object around the Sun, are shown in Fig. 1, where
[phi]@1, [phi]@2, and [phi]@3 denote starting values of angles between a
fixed axis @Ox and directions to the objects. At @N@r = 2@10@5 maximum
eccentricity exceeded 0.04, and perihelion distance of one of the objects
decreased from 50 to 45 a.u. (see Fig. 1a). Because of close encounters
of the objects, the specific appearance of the plots in Figs. 1a - 1c
depends on the relative tolerance [epsilon]@r used in a particular run.
Nevertheless, the trends and the limits of variation of orbital elements
were approximately the same for [epsilon]@r@[10@[-9], 3@10@[-7]].
The accuracy of obeying the conservation laws in such computations
has been examined in my earlier paper [11].

For large time intervals I have used the method of action spheres to
study the evolution of orbits of three objects with masses close to that of
Pluto in a series of computations with different starting data [4]. An
example is shown in Fig. 2: In the course of 10@6, 3@10@6, and 10@7
revolutions the maximum eccentricity @e@[max]
grows up to 0.1, 0.2, and 0.4, and, if semimajor axes
of the initial orbits are close to 50 a.u., the minimum perihelion distance
@r@[pi]@[min]
equals 45, 30, and 15 a.u., respectively. In a number of other samples in
this series of computations the growth of @e@[max] was not as large as
above. If gravitational interactions of objects outside their action
spheres were taken into account, the growth rate of @e@[max] might be
increased. The positions of the objects on their orbits at the instance
they enter the sphere of action were chosen independently of their starting
positions. In the algorithm employed they were specified using a generator
of pseudorandom numbers [3, 15].

The semimajor axis of one of the objects in Fig. 2 did not alter over a
period of time. During this period the aphelion of the second object was
close to the perihelion of the third object. Thus, for this pair of
objects, the sum of angles with their apices at the Sun, inside which the
separation of objects is less than @r@[sa], is considerably greater than
that for
two other pairs of objects [4]. Therefore, for the second and third objects
the probability of approaching each other up to the distance @r@[sa] was
substantionally greater than the probabilities of their approaching the
first object. In addition, in the case that one of the orbits is "steady"
and the other is "varying", relative velocities of the objects approaching
each other are considerably greater, than in the case of two "varying"
orbits.

Still more rapid increase of @e@[max] and decrease of @r@[pi]@[min] are
obtained
in studying the evolution of four objects: two of them with masses close to
the mass of Pluto, and two other with masses several times smaller [4].
Semimajor axes of neighbouring orbits, which were initially almost
circular, differ by 0.1%. Under assumption that the directions to perihelia
are subject to accidental perturbations during the periods between
encounters, the increase of @e@[max] is accelerated [4]. Such perturbations
can be caused, for instance, by the secular perturbations of planets.

In some cases @e@[max] for three interacting objects exhibited scarcely any
growth over a protracted period (see Fig. 1a at @N@r > 2@10@5). The more is
the number of objects involved, the less is the duration of such
"stabilization" of orbits.

The obtained outcomes suggest one of the possible mechanisms of migration
of bodies from the belt beyond Neptune to the Neptune's orbit. The
gravitational influence of the giant planets is regarded as another
mechanism of migration [28].

My computations were performed for the planar model. A consideration of
orbit inclinations can diminish the rate of eccentricity growth.
Nevertheless, the more is the number of bodies, the greater is the increase
of @e@[max] (see below). Additionally, by the bodies of the beyond-Neptune
belt the eccentricities of their orbits could be enhanced in the course of
the formation of giant planets under the effect of planetisemals coming
into the belt from the feeding zones of these planets.

According to Eneev [23], Fernandez [29], and Whipple [41], initial
solar distances of the bodies in the belt beyond Neptune fall predominantly
in the range between 40 and 60 a.u., and the overall mass of these bodies
@M@[Sigma] is greater or equal to 10 Earth masses @m@. On the other hand,
Hamid et al. [36] reasoned that @M@[Sigma]@1.3m@ and @a@50 a.u. If the
starting circular orbit has @a@40 a.u., and in the course of evolution the
aphelion distance @r@[alpha] remains constant, then the object in its
perihelion reaches the Neptune's orbit at @a@35 a.u. and @e = 1/7.

It is Eneev's opinion that Pluto was initially resided in the
beyond-Neptune belt, but gravitational interactions with other objects with
masses close to the Pluto mass changed the orbit of Pluto, so that its
close encounters with Neptune became possible [23].

TYPICAL VARIATIONS OF ORBITAL ELEMENTS IN A SINGLE ENCOUNTER

Using the method of action spheres I have simulated multiply repeated
individual encounters of two objects orbiting the Sun in the same plane.
Each run comprised up to several hundred thousands of such encounters, and
for a variety of starting data the following quantities were computed
[5, 6, 8, 10]:
@de = @, @da = @, @de@+ = @, @de@*, and @N@e@ = .
Here @e@1@, @e@2@ and @a@1@, @a@2@
are the alterations of eccentricities and semimajor axes,
respectively, of the both objects in the @kth approach each other up to
the sphere of action radius @r@[sa]; @N@e is the total number of
encounters; @N@e@c is the number of collisions; @de@* is the value of @de
for bodies of terrestrial density with @a@1 = 1 a.u. Encounters ending in
collision of bodies were ignored in calculating @de@*. The results for many
runs, or variants, of such computations are exemplified in Tables 1 and 2.
In these tables values of @N@e@ are given for the case
@a@* = (a/a@)([rho]/[rho]@)@[1/3] = 1, where [rho]@
is the terrestrial density, @a@ is the semimajor axis of the Earth orbit,
and [rho] is the bodies' density.

The directions pointing toward the perihelia of the object's orbits were
adopted with pseudorandom number generators. In some variants of
calculations the value of [epsilon]@0 was the same for all @k. In the
variants, which are marked "[Sigma]" in both tables, the quantity
1 + [epsilon]@0 took on different values between
@a@[min]@ = (1 - [mu]@[2/5])(1 - e@1@0)/(1 + e@2@0)
and @a@[max]@ = (1 + [mu]@[2/5])(1 + e@1@0)/(1 - e@2@0)
allowing the objects to come together up to the distance @r@[sa]. In the
variants labeled "[Sigma]" the values of [epsilon]@0 were calculated either
from the formula
1 + [epsilon]@0 = a@[min]@ + (a@[max]@ - a@[min]@)[eta]
or from 1 + [epsilon]@0 = @,
where [eta] took on a psevdorandom value equidistributed in the segment
[0, 1]. Both of these formulas yield about the same results of computations
[10]. Masses of objects and eccentricities of their orbits were taken equal
for all encounters within every variant.

For some number of small values @e@0/h@1@8, where @e@0 = max(@e@1@0,
@e@2@0) and @h@1 = ([mu]@1/3)@[1/3], the values of @e@2@k were computed by
numerical integration of the equations of motion in [32, 33]. The computed
values were in the range between 5 and some fraction of
unity, depending on starting orientations and semimajor axes of orbits.
According to [32], @de@2 for @e@0 = 0 and [mu]@1@[mu]@2, whereas @de@0.4
for @e@0/h@1 = 2. Values of
@de@r = de[mu]@1/([mu]@1 + [mu]@2)h@1 obtained by me for
@e@0/h@1@10 in variants labeled "[Sigma]" fall within the range 0.7 - 1.4.
Thus, they are of the same order of magnitude as those obtained by
numerical integration. For orbits with nonzero eccentricities my values of
@de@r are slightly greater than those from [32]. This can be attributed to
the fact that in [32] all conceivable encounters were used to determine
@de, and not just the approaches to the distance @r@[sa], the share of
which may be small.

In the case @e@0/h@1 > 20 the values of @de@r decrease with the increase of
@e@0/h@1. For the range 200@e@0/h@1@3000 holds the relationship
@de@e@0@[-[lambda]] with [lambda]@0.7 - 0.85. For equal values of @e@0/h@1
with the constraint @e@1@0/e@2@0 = const the ratios @de/[mu]@1@[1/4] are
approximately equal. In the variants that differ from each other only by
values of @e@1@0/e@2@0 < 1 the value of @de can vary within a factor of two.
If [epsilon]@0 is varied also, they can vary within a factor of 10 (cf. the
12th and 14th runs in Table 2).

If @e@0/h@1 exceeds 30, mean variations of eccentricities in encounters are
the more, the less is the minimum separation of the objects. Therefore, the
values of @de for mass points are greater than the corresponding values of
@de@* obtained with @a@* = 1.

In the variants labeled "[Sigma]" for the case @e@1@0 = e@2@0 = e@0
I obtained @da/de@1.5 - 2 with @e@0/h@1@10 @da/de@3 - 4 with @e@0/h@1@3
[5]. Under the variation of [epsilon]@0 this ratio may vary several-fold.
The ultimate values of @da are attained, when the aphelion of one orbit is
close to the perihelion of the other.

Analytical estimates show that @N@e@*@e for
1 + 2[Theta]@2[Theta], and @N@e@a@* for 1 + 2[Theta]@1 [12].
Here [Theta] = 0.5(v@p/v@r)@2, @v@r is the
relative velocity of the objects separated by the distance @r@[sa],
@v@p = @, @G is the gravitational constant, @M@ is the mass of the Sun, and
@r@[Sigma] is the sum of the objects' radii. In the case @a@*@1 it was
derived from numerical computations that for @e@0/h@1 < 20 the relation
@N@e@ is accurate within 5%, and for @e@0/h@1 > 5000 the relation @N@e@a@*
is accurate within 3%. The accuracy of the relation @N@e@e
for @a@*@1 and @e@0/h@1 < 30 is within 23%.

Values of @de obtained in studies of the evolution of flat rings comprising
one hundred identical material points orbiting the Sun are tabulated in
Table 3. The method of spheres of action was used, and in each variant of
these computations material points came 10@4 times up to the distance
@r@[sa] of each other. Initial eccentricities @e@0 of their orbits were
equal, and directions pointing towards perihelia of initial orbits were
pseudorandom numbers uniformly distributed in the range between
-180@[degree] and 180@[degree]. Initial semimajor axes were computed from
the formula @a@0@n =@ with @n = 1,..., 100, @a@[min]@ = 1 - @e@0, and
@a@[max]@ = 1 + @e@0.

These computations revealed that after 10@4 encounters the difference
between average eccentricity @e@[av] and @e@0 did not exceed 0.025, if
[mu]@1@3@10@[-10], where [mu]@1 is the ratio of the mass of one of the
material points to that of the Sun. Therefore, @e@[av]@e@0 in most variants
presented in Table 3.

The values of @de tabulated in Table 3 are by a factor of 1.5 - 2.1 greater
than those in Table 1, if @e@0/h@1 > 40. This is due to the fact that in a
disk consisting of material points the probability of an encounter with
large eccentricity variation is greater than that with small one. For equal
values of @e@0 and [mu]@1 the values of @N@e@ obtained in simulation of a
disk's evolution were one half of those tabulated in Tables 1 and 2.

In the case 50@e@0/h@1@1000 two relations were fulfilled simultaneously:
@de@[mu]@1@[1/4] for equal values of @e@0/h@1, and @de@[mu]@1@[1/2]
for equal values of @e@0. At lower values of @e@0/h@1 the dependence of @de
on [mu]@1 is more smooth. For example, with [mu]@1 = 3@10@[-6], the values
of @de for @e@0 = 0.3 were only by a factor 1.5 smaller than those for
@e@0 = 0.02 (cf. the last column in the Table 3). For @e@0@0.1 and
[mu]@1@3@10@[-12], i.e., for @e@0/h@1@1000, an approximative relation was
derived: @de@[mu]@1@[xi]/e@0@[lambda], with [xi]@0.5 and [lambda]@1.
At lower values of @e@0/h@1 the dependence of @de on @e@0 is not as strong.
The outcomes of studying @de and @N@e@ as functions of starting data were
used in my estimates of the evolution of disks comprising many objects [13].

In most variants of computations average eccentricities grew in the course
of evolution. It was found from several variants that for @e@[av]/h@1@10
@de@+ averages 0.06, but because of poor statistics the value of @de@+
diverged considerably in different variants [10].

When @e@[av] grows from 2@h@1 to 10@h@1, the value of @de@+ is greater
than 0.1. Values of @de@+ obtaind in studies of multiply repeated
encounters of objects are greater than those found in studies of disks.
This is due to the fact that the values of @de@+ are small for those
orientations of orbits, where the probability of approach up to @r@[sa] is
large. For instance, if the aphelion of one orbit is close to the
perihelion of the other one, then @de@+ < 0 [5]. When @e@[av] equals
approximately 0.3 - 0.4, intensive ejection of objects into hyperbolic
orbits occurs. Hence, as a rule, values of @e@[av] greater than 0.4 - 0.5
are unattainable [8].

EVOLUTION OF DISKS OF MATERIAL POINTS

The evolution of a flat disk involving @N uniform objects, i.e. material
points orbiting the Sun, was investigated in [12]. From the obtained
results one can derive the variation of the average eccentricity @e@[av]
during a time interval [Delta]@t:

@.(1)

Here @dN@c is the number of objects' approaches up to @r@[sa] during
[Delta]@t, [Delta]@t@ = [Delta]@t/t@0, @t@0 is the period of revolution
of a mass point in the middle of the disk around the gravitating center,
and @[kappa]@N is the mean number of mass points crossing the orbit of one
of them. A typical time interval between successive encounters of objects
in a spatial disk is approximately by a factor of

@(2)

greater than that for a planar disk, with the constraint @i@[av] > r@s@*
[12]. Here @i@[av] is the average inclination of orbits in radians, [eta]@i
= 0.5 + ln(2@i@[av]/r@s@*), @r@s@* = r@[sa]/R, and @R is the solar
distance [12]. Therefore, in the case of a spatial model
@de@[av]@dede@+@[kappa]N(r@s@*)@2[eta]@i[Delta]t@/40i@[av].
It sould be noted that in the spatial case @de can be several times (up to
5 times) smaller than in the planar model.

The above estimates of the variation rate of @e@[av] were derived for a
probabilistic model. In this model the pair of objects approaching each
other up to @r@[sa] was chosen in proportion to the probability @p@[ij] of
their encounter. Actually, an encounter occurs in that pair of objects,
which have the minimal time [tau]@[ij]@1/p@[ij] till the encouter. In the
latter case the evolution time of a disk consisting of several hundreds of
objects can be one order of magnitude less than in the probabilistic model
[18].

As stated above, for large values of @e@[av] we obtained
@de@[mu]@1@[1/2]/e@[av]@[lambda], where [lambda]@0.7. In this case, in view
of @[kappa]@e@[av], one finds from (1) that
@de@[av]2n[mu]@[13/10]e@[av]@[1-[lambda]]@N[mu]@[4/3]e@[av]@[1-[lambda]]@
M@[Sigma][mu]@[1/3]e@[av]@[1-[lambda]],
where @M@[Sigma] is the overall mass of @N objects. From these estimates it
follows, in particular, that the increase of @e@[av] for @N = 10 and [mu]@1
= [mu]@* is as large as that for @N = 10@3 and [mu]@1 =
[mu]@*@10@[-3/2]@[mu]@*/30.

In a planar model with @N = 3 and @a@50 a.u. the time taken to increase
@e@[max] from 0 to 0.2 is one order of magnitude less than the lifetime of
the Solar System. With the use of formula (2) one finds that for
@i@[av]@10@[degree] and for bodies with masses equal to that of Pluto the
mean time interval between approaches to the distance @r@[sa] is about 300
times greater than in a planar model. Consequently, for @i@[av]@10@[degree]
the increase of @e@[max] up to 0.2 in the lifetime of Solar System occurs
at @N@100. It seems likely that bodies migrate to Neptune mainly from the
inner part of the beyond-Neptune belt, and they have rather small
eccentricities and inclinations. Because of mutual gravitational influence
some bodies from the outer part of the belt beyond Neptune can migrate into
its inner part without acquiring large eccentricities.

Although the formula (1) is derived for a disk of bodies, let us compare
the estimates of [Delta]@e@[av] obtained from this formula with the results
of computations shown in Fig. 1. Taking into consideration the data
presented in Table 3 and setting @de@[mu]@1@[1/2], we get @de@0.0005 for
@e@0.01. Let as assume that [Delta]@e@+ = 0.1, @[kappa] = 1, and @N = 3;
then, with [mu]@1 = 7.5@10@[-9] and [Delta]@t = 10@5, we obtain from (1)
that @de@[av]@0.01, in accordance with the data displayed in Fig. 1.

If the eccentricity and inclination of a protoplanet gravitationally
interacting with many small planetesimals are greater than the average
eccentricity and average inclination, respectively, they decrease in the
course of time, as it was demonstrated by computations [3, 15, 17, 34]. In
the course of evolution of a disk consisting of bodies with different
masses the increase of the orbital eccentricities and inclinations is
smaller for the greater bodies. Therefore, the probability of migration of
the largest objects from the beyond-Neptune belt into the zone of the giant
planets is less than that of smaller objects. The orbit inclinations of the
largest objects are likely to be comparatively small, and the rate of the
orbit elements variation for two gravitationally interacting objects is the
greater, the smaller is the angle between the planes of their orbits. Thus,
the objects with small inclinations should migrate from the beyond-Neptune
belt first. The object 1992 QB1 with @a = 44.4 a.u., @e = 0.11, and @i =
2.2@[degree], discovered in this belt [28], has a small inclination, and
its diameter @d is about 240 km [26]. Other known objects that came from
the belt beyond Neptune, e.g., 944 Hidalgo (@a = 5.86 a.u., @e = 0.66, @i =
42@[degree], @d@40 - 60 km) and 2060 Chiron (@a = 13.7 a.u., @e = 0.38, @i =
6.9@[degree], @d@310 - 400 km) [39], changed strongly their orbits under
the influence of giant planets.

The evolution of disks comprised of one hundred identical material points
with circular initial orbits was investigated in [3].
In the course of evolution such disks expanded,
orbital eccentricities of the majority of objects became close to the average
eccentricity, and the objects with small values of @e remained minor in
number.

Also the evolution of a disk with overall mass greater than that of
asteroidal belt by a factor 10@5 and mass distribution of large objects
close to that of large asteroids in the main belt was investigated in [3].
The average eccentricity, being initially equal to 0.15, shows a rise in
the course of evolution, whereas the eccentricities of the largest
objects decrease. The overall expansion of the disk is accompanied by
focussing of massive material points with respect to @a.
The eccentricities of orbits
at the outer edge of disk are @0.4 - 0.6 and exceed appreciably
@e@[av]@0.15 - 0.2, just as in the case of
disks comprised of equal masses.
 
The obtained estimates indicate that, if the average eccentricity of
asteroids in the main belt is to increase in the lifetime of Solar
System from zero to 0.15, i.e. to the value of
contemporary average eccentricity, then
the mass of this belt should be considerably greater than the present-day
mass of the asteroidal belt [3]. This fact, as well as large value
(37.4@[degree]) of the orbit inclination of Pallada, which is one of the
largest asteroids, argues in behalf of formation of the asteroidal belt
under the action of bodies from the giant-planets zone [17, 22].

RELATIVE MOTION OF ENCOUNTERING OBJECTS

If the method of spheres of action is used, the relative motion of objects
inside this sphere is close to actual one provided that the eccentricity
[zeta] of the orbit, along which one of the objects moves about the other
inside of the action sphere, is greater than 4 [27]. Moreover, [zeta]
grows along with @e, and the simulated motion comes more close to the
actual one. Using the method of action spheres I investigated the
evolution of orbits of three objects with masses equal to that of Pluto. It
was demonstrated that the percentage @N@[[zeta]4] of encounters with
[zeta] > 4 amounts about to 90% under the increase of @e@[av] from 0.01 to
0.05. Furthermore, if the everage eccentricity is greater than or equal to
0.05, i.e., if @e@[av]/h@1 is greater than 35, then in the variants
considered @N@[[zeta]4]@99% and, in addition, [zeta] > 100 in more than 80%
of encounters. In this case the motion of an object inside of the sphere of
action for the most part is close to the actual one.

For a number of initial data values of @N@[[zeta]4] in the case of
evolution of orbits of two objects are presented in [4, 5]. The mean values
of @N@[[zeta]4] are roughly 99, 80, and @4%, if @e@[mu] is equal to 50, 10,
and @3, respectively; here @e@[mu] = e@M/[mu]@[1/3] and @e@M is the maximum
eccentricity of two objects approaching each other. At fixed values of
@e@[mu] the values of @N@[[zeta]4] depend on the starting orientations of
orbits and on [epsilon]@0. For example, with @e@1 = e@2 = 0.05 and [mu]@1 =
[mu]@2 = 10@[-5] I found that @N@[[zeta]4] = 9% for [epsilon]@0 = 0.01 and
@N@[[zeta]4] = 0 for [epsilon]@0 = 0.1 [5]. This distinction is due to the
fact that with [epsilon]@0 = 0.1 encounters occur, when the aphelion of the
second object's orbit and the perihelion of the first one are in close
proximity; thus, relative velocities are small. As pointed out above, for
@e@[mu] > 2 the limits and the major patterns of the heliocentric orbital
elements variation are much the same with the method of action spheres and
with the use of numerical integration of the equations of motion. In this
case the values of @N@[[zeta]4] can approach zero.

Figure 3 illustrates the relative motion of two objects in Cartesian
coordimate system. Both objects move around the Sun along initially
circular orbits. The origin is placed at the first object. The close
encounter displayed in Fig. 3 occurs at @t@75t@0, where @t@0 is the period
of revolution of the first object around the Sun. The spacing between fat
points marked on the curve corresponds to the time interval 0.1@t@0. After
the encounter shown in this figure the eccentricities of heliocentric
orbits of both objects increased up to 0.001, but in the course of the
encounter they peaked at 0.06 during a short time.

Although the curve shown in Fig. 3 is complex-shaped, it can be seen that
within the sphere of the radius @R[mu]@1@[1/4] it can be approximated by a
portion of a conic section. During the motion inside this sphere the
objects enter the sphere of action twice. Note that the directions of the
relative velocity vector at the first entrance and at the second egress are
closer to those in the collision of elastic spheres instead of the
direction of object motion along a conic within the sphere of action. With
@a@1@0 = 50 a.u. the minimum separation between the objects in this
encounter is @r@[min] = 2@10@[-4] a.u. = 3@10@4 km, or one order of
magnitude greater than the diameter of Pluto.

Various cases of relative motion of objects approaching each other were
considered in studies of planet rotations [21, 37] and collision
probabilities [35, 38]. I used the method of action spheres in studying the
formation of regular components of planet rotations caused by the incidence
of small bodies [5, 6, 9]. It was found that if the eccentricities of
heliocentric orbits of encountering bodies are small, then positive spins of
forming bodies predominate. The fraction of planetesimals that bring
positive axial momentum to a planet diminishes with the increase of their
eccentricities. Similar results were obtained by numerical integration of
the three-body problem [37]. In contrast to numerical integration, the
method of action spheres predicts a predominance of positive axial
momenta of colliding solid bodies that were moving in the Earth zone along
initially circular orbits. Basically, this distinction can be explaned as
follows. In the method of action spheres it is supposed that the orbits stay
osculating circular up to the instant when the separation @dR of bodies
becomes equal to @r@[sa]. To the contrary, in numerical integration such
assuption holds only if @dR is many times @r@[sa], since at @dR = r@[sa]
the orbits are already weakly eccentric.

Relying on the outcomes of my computations with the method of action
spheres, I have drawn some conclusions concerning the eccentricities of
orbits of small bodies [6, 9]. If the rotation of planets is formed under
the incidence of small bodies alone, then the eccentricities of orbits are
bound to be moderate for the majority of small bodies incident upon a
planet for a
number of planets, with exception of Mercury, Venus, Mars, and Pluto. On
the other hand, the eccentricities of orbits for the majority of small
bodies incident on giant planets were, most likely, not small, because of
gravitational perturbations caused by these planets. These deductions argue
for the conjecture of A.V. Vityazev and G.V. Pechernikova [2]. They
supposed that massive solid bodies made a major contribution to the forward
rotation of some planets, e.g., of Neptune.

T.M. Eneev and N.N. Kozlov developped a model that considers the formation
of giant rarefied protoplanets as pooling dilute gas and dust clusters,
which were formerly moving around the Sun in near-circular orbits [24].
Axial momenta of such giant protoplanets are considerably greater than
those of actual planets. For example, suppose that the rarefied protoplanet
and the clusters were homogeneous spheres coinciding with their Hill
spheres, i.e., their radii are @r@s = R([mu]/3)@[1/3]. Then the axial
momentum of the rarefied proto-Earth should be two orders of magnitude
greater than that of the present-day Earth [6]. But in the course of the
contraction of a giant rarefied protoplanet its axial momentum diminishes
substantionally due to the tidal interaction with the Sun, as it was
recognized by
V.V. Beletskii and A.V. Grushevskii [1]; thus, the results obtained in [6]
conform with the idea of planet formation by contraction of these giant
protoplanets.

CONCLUSION

If starting data are alike, graphical representations of results obtained
with the method of spheres and by numerical integration of the equations of
motion indicate that the main trends of the evolution and the limits of the
orbital elements variation for two gravitating objects, which move in
crossing orbits, are much the same. In the method of spheres an action
sphere is to be used, if the maximum eccentricity of heliocentric orbits of
the encountering objects @e@m is greater than 2[mu]@[1/3], and a sphere of
the radius@r@s = R[gamma]@2.4R[mu]@[1/3] for smaller eccentricities; here
@R is the solar distance of objects, and [mu] = [mu]@1 + [mu]@2 is the
ratio between the sum of object masses and the mass of the Sun.

Furthermore, the evolution of orbits for three identical material points
(MP) with masses equal to that of Pluto was investigated. These MPs were
supposed initially moving in close circular orbits. By numerical
integration of the planar four-body problem it was found that the maximum
eccentricity e@[max] grows up to 0.04 in the course of 2@10@5 revolutions
about the Sun. With the method of action spheres in one of the computed
variants for the planar model it was found that e@[max] increased up to 0.1
and 0.4 during 10@6 and 10@7 revolutions, respectively. When two MPs of
the same masses were considered, e@[max] was less than 0.02. In some
computed variants for three MPs e@[max] exhibited scarcely any growth over
a long period of time. The augmentation of the number of MPs taken into
consideration inhibits such stabilization of orbits. These results suggest
that certain bodies in the beyond-Neptune belt can migrate to the Neptune's
orbit by virtue of the gravitational interaction.

For the MPs moving inside the sphere of action the mean variations of
eccentricities and semimajor axes of their heliocentric orbits were
considered; particularly, the dependence of these quantities on starting
data, i.e., on masses, eccentricities, and semimajor axes. The outcomes of
these considerations, along with the analytical formulas derived, were
applied to evaluate the evolution of disks consisting of many objects,
e.g., of planetesimals, asteroids, or bodies from the beyond-Neptune belt.
It was found, for example, that, if the eccentricities are not small, then
for a disk consisting of 10 MPs and for another consisting of 1000 MPs, but
with masses 30 times less than those in the first one, the growth of the
average eccentricities was about the same. Contemporary eccentricities and
inclinations of asteroidal orbits could not be acquired owing to their
gravitational interaction.

Several cases of the objects' motion inside the action spere were
considered. If the average eccentricity @e@[av] of the objects' orbits
exceeds 35@h@1, where @h@1 = ([mu]/3)@[1/3], then inside the sphere of
action one of the objects moves about the other along an orbit, which
eccentricity exceeds 4 in 99% of encounters. In this case, using the
method of spheres, one obtaines the relative motion of objects inside of
the action sphere close to the actual one.

ACKNOWLEDGEMENTS

I would like to thank V.N. Zharkov, I.N. Ziglina, and T.M. Eneev for a
number of useful remarks.

This work was supported by the Russian Foundation for Fundamental Research,
project no. 93-02-17035.

REFERENCES

1. Beletskii, V.V. and Grushevskii, A.V., Structure of Rotational Motions
of Celestial Bodies and the Formation of Solar System,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1991, no. 103, p. 32.

2. Vityazev, A.V., Pechernikova, G.V., Solution of the Problem of Planetary
Rotation in the Framework of the Statistical Theory  of Accretion,
@Astronomicheskii Zhurnal@, 1981, vol. 58, no. 4, pp. 869 - 878.

3. Ipatov, S.I., The Evolution of a Flat Ring of Gravitating Mass Points,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1978, no. 2, p. 60.

4. Ipatov, S.I., On an Approximate Method for Studying Mutual Gravitational
Influence of Bodies in a Protoplanet Cloud: A contribution to the Problem
of the Pluto's Orbit Evolution,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1980, no. 43, p. 33.

5. Ipatov, S.I., Numerical Investigations of Rotational Momenta of
Accumulating Bodies,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1981, no. 101, p. 28.

6. Ipatov, S.I., Some Questions of the Formation of Axial Rotation of
Planets,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1981, no. 102, p. 28.

7. Ipatov, S.I., Computer Simulation of the Evolution of Flat Rings of
Gravitating Particles Orbiting the Sun,
@Astronomicheskii Zhurnal@, 1981, vol. 58, no. 5, pp. 1085 - 1094.

8. Ipatov, S.I., Evolution Estimates for Spatial Rings of Gravitating
Bodies Coaggregating under Collisions,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1982, no. 186, p. 28.

9. Ipatov, S.I., Axial Rotations of Accreting Planets,
@O.Yu. Shmidt i Sovetskaya Geofizika 80-kh godov@
(O.Yu. Schmidt and Soviet Geophisics in the 1980s),
Moscow: Nauka, 1984, pp. 239 - 243.

10. Ipatov, S.I., Evolution of Bodies' Eccentricities at the Initial Stage
of Solid Bodies' Accumulation into Planets,
@Preprint of Keldysh Inst. of Appl. Math., USSR Acad. Sci.@, Moscow,
1985, no. 4, p. 28.

11. Ipatov, S.I., Numerical Investigations of Possible Orbits' Evolution of
Pluto and Bodies in beyond-Neptune Belt,
@Kinematika i Fizika Nebesnykh Tel@, 1988, vol. 4, no. 6, pp. 73 - 78.

12. Ipatov, S.I., Time Scales of Planetesimal Disks Evolution,
@Astronomicheskii Zhurnal@, 1988, vol. 65, no. 5, pp. 1075 - 1085.

13. Ipatov, S.I., Evolution of the Orbital Eccentricities of Planetesimals
during Formation of the Giant Planets,
@Astronomicheskii Vestnik@, 1989, vol. 23, no. 3, pp. 197 - 206.

14. Ipatov, S.I., Variations of Eccentricities of the Asteroidal Type
Orbits in the Vicinity of the 2:5 Resonance,
@Pis'ma v Astronomicheskii Zhurnal@, 1989, vol. 15, no. 8, pp. 750 - 760.

15. Ipatov, S.I., Evolution of Orbits of Growing Giant-Planet Embrios
Initially Moving along Strongly Eccentric Orbits,
@Pis'ma v Astronomicheskii Zhurnal@, 1991, vol. 17, no. 3, pp. 269 - 281.

16. Ipatov, S.I., Numerical Model of the Evolution of Asteroid Orbits at
the 2:5 Resonance,
@Astronomicheskii Vestnik@, 1992, vol. 26, no. 6, pp. 26 - 53.

17. Ipatov, S.I., Migration of Bodies in the Accretion of Planets,
@Astronomicheskii Vestnik@, 1993, vol. 27, no. 1, pp. 83 - 101.

18. Ipatov, S.I., Techniques of Choosing Pairs of Contacting Objects in
Studies of the Evolution of Discrete Systems with Binary Interactions,
@Matematicheskoe modelirovanie@, 1993, vol. 5, no. 1, pp. 35 - 59.

19. Ipatov, S.I., Gravitational Interaction of Planetesimals Moving in
Close Orbits,
@Preprint of Keldysh Inst. of Appl. Math., Russ. Acad. Sci.@, Moscow,
1994, no. 23, p. 60.

20. Ipatov, S.I., Gravitational Interaction of Two Planetesimals Moving in
Close Orbits,
@Astronomicheskii Vestnik@, 1994, vol. 28, no. 6 (in press).

21. Kiladze, R.I., @Sovremennoe Vrashchenie Planet kak Resultat Razvitiya
Okoloplanetnykh Roev Melkikh Chastits@ (Contemporary Rotation of Planets as
a Result of Evolution of Small Particle Swarms nearby Planets), Tbilisi:
Metsniereba, 1986.

22. Safronov, V.S., @Evolyutsiya Doplanetnogo Oblaka i Obrazovanie Zemli i
Planet@ (Evolution of Pre-Planet Cloud and Formation of the Earth and
Planets), Moscow: Nauka, 1969.

23. Eneev, T.M., On Possible Structure of External (beyond Neptune) Domains
of the Solar System,
@Pis'ma v Astronomicheskii Zhurnal@, 1980, vol. 6, no. 5, pp. 295 - 300.

24. Eneev, T.M. and Kozlov, N.N., A Model of Planetary Systems Formation in
the Course of Accretion, I: Computational Experiments,
@Astronomicheskii Vestnik@, 1981, vol. 15, no. 2, pp. 80 - 94.

25. Bulirsch, R. and Stoer, J., Numerical Treatment of Ordinary
Differential Equations by Extrapolation Methods,
@Numerische Mathematik@, 1966, vol. 8, no. 1, pp. 1 - 13.

26. Cowen, R., Distant Object Hints at the Kuiper Belt,
@Science News@, 1992, vol. 142, no. 10, p. 196.

27. Cox, L.P., Lewis, J.S., and Lecar, M., A Model for Close Encounters in
the Planetary Problem,
@Icarus@, 1978, vol. 74, no. 2, pp. 415 - 427.

28. Duncan, M. and Quinn, T., The Longterm Dynamical Evolution of the Solar
System,
@Annu. Rev. Astron. Astrophys.@, 1993, vol. 31, pp. 265 - 296.

29. Fernandez, J.A., On the Existence of a Comet Belt beyond Neptune,
@Monthly Notices Roy. Astron. Soc.@, 1980, vol. 92, no. 2, pp. 481 - 491.

30. Greenberg, R., Carusi, A., and Valsecchi, G.B., Outcomes of Planetary
Close Encounters: A Systematic Comparison of Methodologies,
@Icarus@, 1988, vol. 75, no. 1, pp. 1 - 29.

31. Greenberg, R., Bottke, W.F., Carusi, A., and Valsecchi, G.B.,
Planetary Accretion Rates: Analytical Derivation,
@Icarus@, 1991, vol. 94, no. 1, pp. 98 - 111.

32. Greenzweig, Y. and Lissauer, J.J., Accretion Rates of Protoplanets,
@Icarus@, 1990, vol. 87, no. 1, pp. 40 - 77.

33. Ida, S., Stirring and Dynamical Friction Rates of Planetesimals in the
Solar Gravitational Field,
@Icarus@, 1990, vol. 88, no. 1, pp. 129 - 145.

34. Ida, S. and Makino, J., N-Body Simulation of Gravitational Interaction
between Planetesimals and a Protoplanet. II: Dynamical Friction,
@Icarus@, 1992, vol. 98, no. 1, pp. 28 - 37.

35. Ida, S. and Nakazawa, K., Collisional Probability of Planetesimals
Revolving in the Solar Gravitational Field. III,
@Astron. Astrophys.@, 1989, vol. 224, no. 1/2, pp. 303 - 315.

36. Hamid, S.E., Marsden, B.G., and Whipple, F.L., Influence of a Comet
Belt beyond Neptune on the Motion of Periodic Comets,
@Astron. J.@, 1968, vol. 73, no. 8, pp. 727 - 729.

37. Lissauer, J.J. and Kary, M., The Origin of Systematic Component of
Planetary Rotation. I: Planet on a Circular Orbit,
@Icarus@, 1991, vol. 94, no. 1, pp. 126 - 159.

38. Nishida, S., Collisional Processes of Planetesimals with a Protoplanet
under the Gravity of the Proto-Sun,
@Prog. Theor. Phys.@, 1983, vol. 70, no. 1, pp. 93 - 105.
    
39. Olsson-Steel, D., Planetary Close Encounters: Probability Distributions
of Resultant Orbital Elements and Application to Hidalgo and Chiron,
@Icarus@, 1987, vol. 69, no. 1, pp. 51 - 69.

40. Wetherill, G.W. and Cox, L.P., The Range of Validity of the Two-Body
Approximation in Models of Terrestrial Planet Accumulation. II:
Gravitational Cross-Sections and Runaway Accretion,
@Icarus@, 1985, vol. 63, pp. 290 - 303.

41. Whipple, F.L., Evidence for a Comet Belt beyond Neptune,
@Proc. National Acad. Science@, 1964, vol. 51, no. 5, pp. 711 - 718.

42. Wisdom, J., Chaotic Behavior and the Origin of the 3/1 Kirkwood Gap,
@Icarus@, 1983, vol. 56, no. 1, pp. 51 - 74.

43. Yoshikawa, M., Motion of Asteroids at the Kirkwood Gaps. II. On the
5:2, 7:3, and 2:1 Resonances with Jupiter,
@Icarus@, 1991, vol. 92, no. 1, pp. 94 - 117.

TABLES

Table 1. Average variations of orbital eccentricities of two objects with
[mu]@1 = [mu]@2 and @e@1@0 = @e@2@0 = @e@0 obtained in simulations of their
multiple approaches each other up to the sphere of action radius @r@[sa]

Columns: 1. Run

Table 2. Average variations of eccentricities of the second object's
orbit obtained in simulations of two objects with [mu]@1@[mu]@2 multiply
approaching each other up to @r@[sa]

Columns: 1. Run

Table 3. Average variations @de of orbital eccentricities of objects
approaching each other up to @r@[sa] obtained in simulating the evolution
of planar disks comprising one hundred identical material points

FIGURE CAPTIONS

Fig. 1. The orbit elements @a, @e, and [pi] of three gravitationally
interacting objects versus the number @N@r of revolutions of the first
object around the Sun. Results obtained by numerical integration of the
planar four-body problem with different tolerances [epsilon]@r:
(a) [epsilon]@r = 10@[-8], (b) [epsilon]@r = 3@10@[-7], and
(c) [epsilon]@r = 10@[-9]. Starting data: [mu]@1 = [mu]@2 = [mu]@3 =
7.5@10@[-9], @a@1@0 = 50 a.u., @a@2@0 = 50.3 a.u., @a@3@0 = 49.7 a.u.,
@e@1@0 = @e@2@0 = @e@3@0 = 0, [phi]@1@0 = 2 rad, [phi]@2@0 = 4 rad,
and [phi]@3@0 = 6 rad.

In Figs. 1a - 1c legends at the vertical axes (four times repeated):
1. [pi], rad
2. @e
3. @a, a.u.

Fig. 2. The orbital elements evolution for three gravitationally
interacting objects. Outcomes of the method of action spheres.
The abscissa is the number of encounters @N@[en] with mutual approach up to
the radius of the corresponding sphere of action. Starting data:
[mu]@1 = [mu]@2 = [mu]@3 = 10@[-8], @a@1@0 = 50 a.u., @a@2@0 = 50.1 a.u.,
@a@3@0 = 49.9 a.u., @e@1@0 = @e@2@0 = @e@3@0 = 0.005, [pi]@1@0 = 1 rad,
[pi]@2@0 = 3 rad, and [pi]@3@0 = 2 rad.

In Fig. 2 legends at the vertical axes:
1. [pi], rad
2. @e
3. @a, a.u.

Fig. 3. The motion of the second object with respect to the first one in
Cartesian coordinates @x, @y.
Starting data: [mu]@1 = [mu]@2 = 7.5@10@[-9], [epsilon]@0 = 0.004,
@e@1@0 = @e@2@0 = 0, [phi]@0 = 160@[degree].
Time intervals between the points marked on the orbit are 10 times smaller
than the period of the first object's revolution about the Sun.
Fig. 3b shows the portion of Fig. 3a inside the dashed rectangle.