INERTIAL CHARACTERISTICS OF HIGHER ORDERS AND DYNAMICS IN A PROXIMITY OF A SMALL CELESTIAL BODY Текст научной статьи по специальности «Физика»

Russian Journal of Nonlinear Dynamics, 2020, vol. 16, no. 2, pp. 259-273. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd200203

NONLINEAR PHYSICS AND MECHANICS

MSC 2010: 70K20, 70K42, 70F05

Inertial Characteristics of Higher Orders and Dynamics in a Proximity of a Small Celestial Body

As is well known, many small celestial bodies are of a rather complex shape. Therefore, the study of the dynamics of a spacecraft in their vicinity, based on terms up to the second order of smallness in the expansion of the potential of attraction, seems to be insufficient for an adequate description of the observed dynamical effects related, for example, to positioning of the libration points.

In this paper, such effects are demonstrated for spacecraft dynamics in the vicinity of the asteroid (2063) Bacchus. The libration points are computed for various approximations of the gravitational potential. The results of this computation are compared with similar results obtained before for the so-called Sludsky-Werner-Scheeres potential. The dependence of the structure of the regions of possible motions on approximation of the gravitational potential is also studied.

Keywords: (2063) Bacchus, gravitational potential expansion, libration points, region of possible motion, Hill's region, zero-velocity locus

Received December 20, 2019 Accepted May 05, 2020

This work is partially supported by the Program of the President of the Russian Federation for Federal Support of Young Russian Scientists, Candidates of Sciences (project no. MK-1712.2019.1) and RFBR (project no. 18-01-00335).

Alexander A. Burov jtm@narod.ru Vasily I. Nikonov nikon_v@list.ru

Federal Research Center "Computer Science and Control" ul. Vavilova 40, Moscow, 119333 Russia National Research University "Higher School of Economics" ul. Myasnitskaya 20, Moscow, 101000 Russia

A. A. Burov, V. I. Nikonov

1. Introduction

As is well known, the expansion of the gravitational potential in a series of harmonic functions is one of the key research tools in celestial mechanics (see, e.g., [1-6]). As is clear, the use of approximations for the gravitational potential of a higher order allows more subtle dynamical effects to be identified.

The first-order approximation of the gravitational potential allows one to identify "Keplerian orbits" located at a sufficient distance from the small celestial body, in particular, stationary orbits for which the spacecraft is located above the same point on the surface in its "equatorial plane".

Using the second-order approximation of the gravitational potential allows one to verify the complex nonlinear dynamics in the vicinity of a small celestial body, which occurs due to the lack of central symmetry in the body's mass distribution. In particular, it turns out that the stationary orbit is decaying, and instead of it, the so-called libration points can be found. The aforementioned libration points are nothing more than equilibrium positions relative to the coordinate frame uniformly rotating with the small celestial body. However, due to the symmetry of the potential with respect to the equatorial plane, the specified libration points are located either in this plane (as is the case in the example considered below) or symmetrically with respect to it.

Using higher-order approximations in the decomposition of the gravitational potential allows one to identify more subtle effects associated with its complete asymmetry, typical for small celestial bodies. In particular, it turns out that the libration points "leave" the equatorial plane. The paper is devoted to the problem of comparing the above-mentioned dynamical effects with similar effects that can be detected when applying the "exact potential", which in this situation can be attributed to the so-called Sludsky-Werner-Scheeres potential (SWS-potential) [7-9].

As is well known, the so-called regions of possible motion closely relate to the libration points. The study of these regions goes back probably to G.W.Hill [10]. Hill paid particular attention to the boundaries of regions of possible motion, calling them loci of zero-velocity (see [10, p. 25]). The refinement of the structure of zero-velocity surfaces due to the refinement of the gravitational potential approximation is also discussed in the work within the example of the asteroid (2063) Bacchus.

Computing or, more correctly, estimating the coefficients of the gravitational potential decompositions is a difficult task due to the fact that information about both the surface geometry of a small celestial body and the geometry of its mass distribution is rather poor. Moreover, the complexity of computing of the corresponding decomposition increases both with increasing order of decomposition and with the refinement of the surface geometry and mass distribution properties. At the same time, while the general formulas for such coefficients are known (see., e.g., [5, 6, 11-16]), the methods of computing the coefficients using approximations of the surface by polyhedra were developed and discussed in numerous publications (see., e.g., [8, 9, 17-25]).

2. Problem statement

The problem of the motion of a spacecraft in the vicinity of a small celestial body B that performs uniform rotations around one of its principal central axes of inertia is considered. As is assumed, the size of the spacecraft is negligible, and its motion does not affect the rotation of the body B. This means that the spacecraft can be considered as a massive point P. The massive point P moves under attraction from the body B. Let Ox\x2x3 be a reference frame fixed in the body B, where the point O coincides with the body's B center of mass. Suppose that the

position of the point P is given by the radius vector = r = (r1, r2, r3)T. In the case when the body B rotates with a constant angular velocity w around the axis Ox3, the motion of the spacecraft P is described by the well-known equations

dU dU dU . .

ri = 2wr2-—, r-2 = —2wri — ——, r3 = - —. (2.1)

dri dr2 dr3

These equations can be written as Lagrange equations with the Lagrange function

L = L2 + L i + Lo, l'2 = \ (r'l + f\ + f'D, L\ = to (r2ri - nf2), L0 = -U. (2.2)

Here the augmented potential

w2

U(r-,u) = Uc + UN, Uc = -—{r\ + rl), UN = UN( r) (2.3)

is formed by the potential of the centrifugal forces Uc and the potential of the attractive forces UN.

Equations (2.1) admit the Jacobi integral

Jo = \ {fl + f'l + rl) +U = h, (2.4)

determining in the configuration space R3 the regions of possible motion

Hh: U(r) < h,

known in celestial mechanics as Hill's regions. Here and below h is a constant of the Jacobi integral. It has units m2fs2.

The boundary of the region Hh is a set of points from R3 such that

dHh: U (r) = h.

It is called a zero-velocity surface. The topology of Hill's regions depends significantly on the so-called libration points, i.e., on equilibrium positions of a massive point P with respect to the uniformly rotating reference frame fixed in the celestial body (see, e.g., [26]). Moreover, the regions Hh significantly depend on the used approximation of the potential of attraction.

In the modern mechanics of small celestial bodies one actively uses the so-called Sludsky-Werner - Scheeres potential [7-9] for description of the field of attraction. This potential is reconstructed from triangulation grids approximating the surfaces of such bodies. The expression for the SWS potential contains a very large number of terms, depending on the number of faces in the triangulation grid and increases rapidly as the size of a simplex decreases. We intend to compare the positions of libration points obtained before using the SWS potential and appropriate positions obtained with use of the potential in approximations up to terms of the third and fourth order of smallness. The small parameter expresses the ratio of the characteristic size of the body to the distance to the studied point in space. This comparison is carried out in Section 3.1 for asteroid (2063) Bacchus. The precalculated inertial parameters are given in the Appendix C. The regions of possible motion calculated for various approximations of the gravitational potential are also compared in Section 3.2.

As is well known, the gyroscopic forces generated by the velocity term L1 in the Lagrange function (2.2) play a special role in the dynamics of mechanical systems. In particular, as is pointed out in [27, 28], when solving the problem of entering the boundary of the trajectory

of a Lagrangian system with a Lagrange function (2.2), the form quadratic in velocity plays a key role

F(ri,r2) = 4(h - Lo)L2 - L\.

In the problem under consideration, this form has the form

( 2(h - U) - u2r22

F = (nr ,r3)

2

w2r1 r2

2

w2r1 r2

2 (h - U) - w2r2

0 0

\ Zri\

r2

V

2 (h - U) ) \h )

This quadratic form is positively defined when the following conditions are fulfilled:

Ko = h - U > 0, Ki = 4 (h - U) - w2 (r2 + r2) > 0, K2 = 2 (h - U) (2 (h - U) - w2 (r2 + r2)) > 0.

(2.5)

(2.6)

These conditions define the Kh domain in the configuration space.

According to [27, 28], if Kh contains a connected compact component of a region of possible motions, then there is a trajectory that connects any internal point of this component to some point of the boundary of the region of possible motion.

As one can see, if the first condition is met, the third condition can be represented as

h -U-^- (r{ + r'l) > 0.

(2.7)

Moreover, if condition (2.7) is met, then conditions (2.5) and (2.6) are also met. Thus, the sufficient condition [27, 28] for reaching the boundary from any inner point of the region of possible motion is not met. The mutual disposition of the regions Hh n {r3 = 0} and Kh n {r3 = 0} is studied in Section 3.3 for the asteroid (2063) Bacchus within the second approximation of the gravitational potential.

2

3. Libration points, regions of possible motions and approximation of the potential

3.1. Sensitivity of the libration points to approximation of potential

In [29] libration points were calculated for the asteroid (2063) Bacchus, based on the assumption of its homogeneity. The Sludsky-Werner-Scheeres approximation of the exact potential [7-9] by the potential was used within these calculations. Quite unexpected, at first glance, was the fact that these points are not located in the "equatorial plane" of the asteroid, i.e., in the plane passing through the center of mass and perpendicular to the axis of rotation. In this paper, libration points are calculated for potential expansions up to terms of the second, third, and fourth order of smallness. The results of calculations are concentrated in Tables 2-5. They provide data for libration points E\, E2, E3 and E4, where E\, E3 are saddle points of potential (2.3), while E2 and E4 are maximum points of potential (2.3). For each of the blocks in these tables, the first, second, and third rows contain the results of calculations carried out using approximations for the potential up to the second, third, and fourth order terms, respectively. The last columns of the tables contain the constant h (introduced above) of the Jacobi integral at the equilibrium point. The fourth row contains data borrowed from [29] and

calculated there using the Sludsky - Werner-Scheeres potential approximation. Naturally, in the second approximation, all libration points turn out to be in the equatorial plane. However, the deviation of libration points from this plane is detected already in the third approximation. The accuracy of calculations can be estimated using the relations

ök =

\гз{к) - r3| (r3(fc)+r3)/2

• 100%,

where k is the order of the potential approximation. position of libration points is shown in Table 1.

Such an accuracy in determining the

Table 1. Libration points: accuracy in determining the position

Ei E2 E3 E4

к = 3: ö3 13.1% 27.8% 28.9% 0.5% к = 4: ö4 5.2% 22.6% 12.2% 6.4%

Table 2. Positions of the equilibrium point Ei for the asteroid (2063) Bacchus

model 2nd order 3rd order 4th order [29]

ri

1.139272396 1.141764567 1.128832898 1.14738,

Г2 0

0.01514411367 0.01741607561 0.0227972

Г3 0

-0.0009823587802 -0.0009071046484 -0.000861348

h

-2.54224903720888•10~2 -2.54551778269219•10~2 -2.53277880925007•10~2 -2.55101029859997-10~2

Table 3. Positions of the Equilibrium point E2 for asteroid (2063) Bacchus

model

2nd order 3rd order 4th order [29]

ri 0

0.02388165946 0.04052974431 0.0314276

r2

1.071115157 1.069949559 1.071108412 1.07239

гз

h

0 -2.43732397327281 • 10~2

0.0009412669899 -2.43639645631026•10~2

0.0008923210581 -2.43781662306405•10~2

0.000711379 —2.43861943518558•10~2

Table 4. Positions of the Equilibrium point E3 for asteroid (2063) Bacchus

model 2nd order 3rd order 4th order [29]

ri

-1.139272396 -1.136529881 -1.123079849 1.14129

r2 0

0.01833173178 0.01696204884 0.00806235

гз 0

-0.001057875823 -0.001252553566 0.00141486

h

-2.54224903720888•10~2 -2.53899333512537•10~2 -2.52593412347194•10~2 —2.54353644094217•10~2

Table 5. Positions of the Equilibrium point E4 for asteroid (2063) Bacchus

model 2nd order 3rd order 4th order [29]

ri 0

0.02599401445 0.02758791033 0.0203102

r2

-1.071115157 -1.071545893 -1.073266552 1.07409

гз 0

0.0008543770594 0.0009063739343 0.000849894

h

-2.43732397327281•10~2 -2.43812423974624•10~2 -2.43961350747501•10~2 -2.44020000723179-10~2

Remark 1. The study of libration points of uniformly rotating celestial bodies dates back probably to publications [30-36], in which the existence and stability of libration points for a uniformly rotating ellipsoid were studied. According to recent studies (cf. [29]), the number of libration points near small celestial bodies can take quite unexpected values. The study of this phenomenon with use of truncated potentials described above is a subject of particular interest.

3.2. Sensitivity of regions of possible motions to approximation of the potential

Consider the boundaries of the regions of possible motion. These surfaces are determined as levels of the augmented potential. Their cross-sections formed by the intersection with the plane r3 = 0, located far from the asteroid's surface and constructed for truncations of the gravitational potential, are visually indistinguishable from appropriate curves in the case of the Sludsky-Werner-Scheeres gravitational potential. Similar sections for second-, third- and and fourth-order approximations of the gravitational potential are presented in Figs. 2-6, respectively.

As one can see in Fig. 2, constructed for the second-order approximation of the gravitational potential, within this approximation the points of libration are located in the equatorial plane of the asteroid, and their positions are symmetric with respect to appropriate coordinate planes.

However, even for the third-order approximation of the gravitational potential this symmetry disappears. In particular, the points of libration no longer belong to the coordinate axes. In order to show this, we draw sections of boundaries of the regions of possible motion which are formed by the intersection with the planes passing through the above-mentioned points of libration — see Figs. 3-4. The magnified regions, close to appropriate points of libration, are drawn in extreme left and right parts of Figs. 3-4 in order to demonstrate the aforementioned deviations of the points of libration from the coordinate axes.

Similar pictures, drawn for the fourth-order approximation of the gravitational potential in Figs. 5-6, demonstrate more precisely appearing asymmetries.

3.3. The reachability of the zero-velocity surfaces for the asteroid (2063) Bacchus

In this section we assume that the motion of the test point P in the vicinity of the asteroid is described within the second-order approximation of the gravitational potential. In this case, r3 = 0 is a symmetry plane of the problem. Then, within this section, we assume that the point P remains in the plane r3 = 0 throughout the motion. This means that at the initial time the conditions r3 = 0, r3 = 0 are true, and for our goals the regions Hh and Kh will be represented by their sections Hh n {r3 = 0} and Kh n {r3 = 0}, respectively.

The value h = -2.52416582487888 • 10_2 of the energy integral is less than its maximum value corresponding to the libration points E2 and E4, but exceeds the value of the energy integral corresponding to the saddle points E1 and E3. This level corresponds to a zero-velocity curve consisting of two closed components rh. This curve is the boundary of the area Hh n {r3 = 0}: rh = d (Hh n {r3 = 0}) (Fig. 7 (left)). The curve Ah = d (Kh n {r3 = 0}), which is the boundary of the region Kh n {r3 = 0}, is also drawn in this figure. In this case, the region of possible motion is a plane without two disks.

In Fig. 7 (center) these regions and their boundaries are drawn for the value h1,3 = = -2.54224903720888 • 10"2 corresponding to the saddle points E1 and E3. In this case, the region of possible motions is a plane without two disks glued together by a pair of points.

Fig. 1. From (a) to (d): sections of zero-velocity surfaces formed by the intersection with the plane r3 = 0 constructed for the second, third, and fourth approximations of the potential, respectively, as well as for the Sludsky-Werner-Scheeres potential. The color code represents the augmented potential per unit mass (cf. Fig. 1i from [29]).

Fig. 2. Sections of zero-velocity surfaces formed by the intersection with the planes r2 = 0 (left) and ri =0 (right), constructed for the second-order approximation of the gravitational potential.

Fig. 3. Sections of zero-velocity surfaces formed by the intersection with the planes r2 = r2(E3) (a) and r2 = r2(Ei) (b), as well as magnified vicinities of E3 (c) and Ei (d) are drawn for the third-order approximation of the gravitational potential.

Fig. 4. Sections of zero-velocity surfaces formed by the intersection with the planes ri = ri (E4) (a) and ri = ri(E2) (b), as well as magnified vicinities of E4 (c) and E2 (d) are drawn for the third-order approximation of the gravitational potential.

0.0015

T3 0.0010

0.0005

0

-0.0005 -0.0010 -0.0015

Fig. 5. Sections of zero-velocity surfaces formed by the intersection with the planes r2 = r2 (E3) (a) and r2 = r2 (E1) (b), as well as magnified vicinities of E3 (c) and E1 (d) are drawn for the fourth-order approximation of the gravitational potential.

0.0015 0.0010 0.0005 0

-0.0005 -0.0010 -0.0015

Fig. 6. Sections of zero-velocity surfaces formed by the intersection with the planes r1 = r1 (E4) (a) and r1 = r1 (E2) (b), as well as magnified vicinities of E4 (c) and E2 (d) are drawn for the fourth-order approximation of the gravitational potential.

- r3 v /

- / \ T2

1.070/ 1.071 \o72

Fig. 7. Zero-velocity curve and the curve Ah are colored in red and black, respectively. Region Kh n {r3 = 0} is white. Regions where no motion is possible are dark-gray. (The colors are valid only for online version of the article.)

Finally, in Fig. 7 (right) these regions and their boundaries are drawn for the value h = = —2.56033224953888 • 10"2 < hi 3. In this case, the region of possible motion is a plane without a ring.

4. Concluding remarks

As is well known, there exist different approaches to approximating the attraction potentials for real celestial bodies. The classical approach provides the use of the truncated potential containing terms up to the second order of smallness. According to V. V. Beletsky and A. V. Rodnikov (cf. [37-39]), the attraction field for various celestial bodies can be described by the potential of pairs of homogeneous balls, making up the so-called "gravitating dumbbell". Within this representation of celestial bodies, numerous qualitative analytical results related to dynamics in their vicinity were obtained. The question is whether comparable qualitative results can be obtained considering the truncated potential with terms up to the third and fourth order of smallness.

References

[1] Duboshin, G.N., Theory of Attraction, Moscow: Fizmatlit, 1961 (Russian).

[2] Brour, D. and Clemence, G., Methods of Celestial Mechanics, New York: Acad. Press, 1961.

[3] Beletskii, V. V., Motion of an Artificial Satellite about Its Center of Mass, Jerusalem: Israel Program for Scientific Translations, 1966.

[4] Demin, V.G., Motion of an Artificial Satellite in an Eccentric Gravitation Field, Springfield, Va.: NASA, 1970.

[5] Dubochin, G.N., Celestial Mechanics: Main Problems and Methods, 3rd ed., Moscow: Nauka, 1975 (Russian).

[6] Doubochine, G.N., Sur le developpement de la fonction des forces dans le probleme de deux corps finis, Celestial Mech., 1976, vol. 14, no. 2, pp. 239-281.

[7] Sludsky, 0., On Deviation of Plumb-Lines, Master's Thesis, Moscow University, 1863 (Russian).

[8] Werner, R. A., The Gravitational Potential of a Homogeneous Polyhedron or Don't Cut Corners, Celest. Mech. Dyn. Astron., 1994, vol. 59, no. 3, pp. 253-278.

[9] Werner, R. A. and Scheeres, D. J., Exterior Gravitation of a Polyhedron Derived and Compared with Harmonic and Mascon Gravitation Representations of Asteroid 4769 Castalia, Celest. Mech. Dyn. Astron., 1996, vol. 65, no. 3, pp. 313-344.

[10] Hill, G. W., Researches in the Lunar Theory, Amer. J. Math., 1878, vol. 1, no. 1, pp. 5-26.

[11] Pines, S., Uniform Representation of the Gravitational Potential and Its Derivatives, AIAA J., 1973, vol. 11, no. 11, pp. 1508-1511.

[12] Schutz, B.E., The Mutual Potential and Gravitational Torques of Two Bodies to Fourth Order, Celestial Mech., 1981, vol.24, no. 2, pp. 173-181.

[13] Sincarsin, G. B. and Hughes, P. C., Gravitational Orbit-Attitude Coupling for Very Large Spacecraft, Celestial Mech., 1983, vol.31, no. 2, pp. 143-161.

[14] Soler, T., A New Matrix Development of the Potential and Attraction at Exterior Points As a Function of the Inertia Tensors, Celestial Mech., 1984, vol. 32, no. 3, pp. 257-296.

[15] Ashenberg, J., Proposed Method for Modeling the Gravitational Interaction between Finite Bodies, J. Guid. Control Dynam, 2005, vol.28, no.4, pp. 768-774.

[16] Ashenberg, J., Mutual Gravitational Potential and Torque of Solid Bodies via Inertia Integrals, Celest. Mech. Dyn. Astron., 2007, vol. 99, no. 2, pp. 149-159.

[17] Fantino, E. and Casotto, S., Methods of Harmonic Synthesis for Global Geopotential Models and Their First-, Second- and Third-Order Gradients, J. Geodesy, 2009, vol. 83, no. 7, pp. 595-619.

[18] Hirabayashi, M. and Scheeres, D.J., Recursive Computation of Mutual Potential between Two Polyhedra, Celest. Mech. Dyn. Astron., 2013, vol.117, no. 3, pp. 245-262.

[19] Wang, Y. and Xu, S., Gravity Gradient Torque of Spacecraft Orbiting Asteroids, Aircr. Eng. Aerosp. Tec, 2013, vol. 85, no. 1, pp. 72-81.

[20] Wang, Y. and Xu, S., Equilibrium Attitude and Stability of a Spacecraft on a Stationary Orbit Around an Asteroid, Acta Astronaut., 2013, vol. 84, pp. 99-108.

[21] Hou, X., Scheeres, D. J., and Xin, X., Mutual Potential between Two Rigid Bodies with Arbitrary Shapes and Mass Distributions, Celest. Mech. Dyn. Astron., 2017, vol. 127, no. 3, pp. 369-395.

[22] Saraswati, A. T., Cattin, R., Mazzotti, S., and Cadio, C., New Analytical Solution and Associated Software for Computing Full-Tensor Gravitational Field due to Irregularly Shaped Bodies, J. Geodesy, 2019, vol. 93, no. 12, pp. 2481-2497.

[23] Dobrovolskis, A. R., Inertia of Any Polyhedron, Icarus, 1996, vol. 124, no. 2, pp. 698-704.

[24] Mirtich, B., Fast and Accurate Computation of Polyhedral Mass Properties, J. Graph. Tools, 1996, vol. 1, no. 2, pp. 31-50.

[25] Werner, R. A., Spherical Harmonic Coefficients for the Potential of a Constant-Density Polyhedron, Comput. Geosci, 1997, vol. 23, no. 10, pp. 1071-1077.

[26] Kozlov, V. V., Calculus of Variations in the Large and Classical Mechanics, Russian Math. Surveys, 1985, vol.40, no. 2, pp. 37-71; see also: Uspekhi Mat. Nauk, 1985, vol.40, no. 2(242), pp. 33-60, 237.

[27] Kozlov, V. and Polekhin, I., On the Covering of a Hill's Region by Solutions in the Restricted Three-Body Problem, Celest. Mech. Dyn. Astron., 2017, vol.127, no. 3, pp. 331-341.

[28] Kozlov, V. and Polekhin, I., On the Covering of a Hill's Region by Solutions in Systems with Gyroscopic Forces, Nonlinear Anal., 2017, vol. 148, pp. 138-146.

[29] Wang, X., Jiang, Y., and Gong, Sh., Analysis of the Potential Field and Equilibrium Points of Irregular-Shaped Minor Celestial Bodies, Astrophys. Space Sci., 2014, vol. 353, no. 1, pp. 105-121.

[30] Zhuravlev, S. G., Stability of the Libration Points of a Rotating Triaxial Ellipsoid, Celestial Mech., 1972, vol. 6, no. 3, pp. 255-267.

[31] Zhuravlev, S. G., About the Stability of the Libration Points of a Rotating Triaxial Ellipsoid in a Degenerate Case, Celestial Mech., 1973, vol. 8, no. 1, pp. 75-84.

[32] Kosenko, 1.1., On Libration Points near a Gravitating and Rotating Triaxial Ellipsoid, J. Appl. Math. Mech., 1981, vol.45, no. 1, pp. 18-23; see also: Prikl. Mat. Mekh, 1981, vol.45, no. 1, pp. 26-33.

[33] Kosenko, I. I., Libration Points in the Triaxial Gravitating Ellipsoid Problem: Geometry of the Stability Region, Kosmicheskie Issledovaniya, 1981, vol.19, no. 2, pp. 200-209 (Russian).

[34] Kosenko, I. I., Non-Linear Analysis of the Stability of the Libration Points of a Triaxial Ellipsoid, J. Appl. Math. Mech., 1985, vol.49, no. 1, pp. 17-24; see also: Prikl. Mat. Mekh, 1985, vol.49, no. 1, pp. 16-24.

[35] Kosenko, I. I., On a Power Series Expansion of the Graviational Potential of an Inhomogeneous Ellipsoid, J. Appl. Math. Mech., 1986, vol. 50, no. 2, pp. 142-146; see also: Prikl. Mat. Mekh, 1986, vol. 50, no. 2, pp. 194-199.

[36] Kosenko, I. I., On the Stability of Points of Libration of an Inhomogeneous Triaxial Ellipsoid, J. Appl. Math. Mech., 1987, vol.51, no. 1, pp. 1-5; see also: Prikl. Mat. Mekh., 1987, vol.51, no. 1, pp. 3-8.

[37] Beletsky, V. V., Generalized Restricted Circular Three-Body Problem as a Model for Dynamics of Binary Asteroids, Cosmic Research, 2007, vol.45, no. 5, pp. 408-416; see also: Kosmicheskie Issledovaniya, 2007, vol.45, no. 5, pp.435-442.

[38] Beletskii, V. V. and Rodnikov, A. V., Stability of Triangle Libration Points in Generalized Restricted Circular Three-Body Problem, Cosmic Research, 2008, vol.46, no. 1, pp. 40-48; see also: Kosmicheskie Issledovaniya, 2008, vol.46, no. 1, pp. 42-50.

[39] Rodnikov, A. V., Triangular Libration Points of the Generalized Restricted Circular Problem of Three Bodies for Conjugate Complex Masses of Attracting Centers, Russian J. Nonlinear Dyn., 2014, vol.10, no. 2, pp. 213-222.

[40] Sretensky, L.N., Theory of the Newton Potential, Moscow: Gostekhizdat, 1946 (Russian).

[41] Burov, A.A. and Nikonov, V. I., Computation of Attraction Potential of Asteroid (433) Eros with Accuracy up to the Terms of the Fourth Order, Dokl. Akad. Nauk. Fizika, Tekhn. Nauki, 2020, vol.492, no. 1, pp. 56-60 (Russian).

[42] Nikonov, V. I., Gravitational Fields of Small Celestial Bodies, Moscow: Belyi Veter, 2020 (Russian).

[43] Small Body Radar Shape Models v2.0., Neese, C. (Ed.), NASA Planetary Data System, 2004.

A. Approximations of the Newtonian potential

According to the general theory (see, e. g., [1, 40]), at the point P, external to the body B, the potential of attraction UN(r) and the gravitational field strength g(r) are given by the relations

b ' '

Here p = p(x) is the body's density and G is the gravitational constant.

The expansion of the potential in a series of harmonic polynomials (see, e. g., [5]), has the

G

UN(r) = U0(r) + ... + Un(r) + ..., = __<?„( r), (A2)

Here and below Cn(r), n = 1,2,... are polynomials in the components of the vector r. Determination of a general form for these coefficients for the first five terms in (A.2) and their computation for the small celestial bodies listed above is a subject for further investigation.

As is well known [2] (see also [41, 42]), the functions Ck(r) read

Co = Jooo, Ci = - ( JiooH + Jowl'2 + Jooir-з),

r

1 3

c2 = -- (J200 + J020 + J002) + ~~9 (Jiwnr2 + Jioinr-i + J0iir2r3) +

2 r2

+ (J200rf + J020^2 + Лэ02?з),

3

C3 = -— [( J-зоо + Л20 + J\02) ri + (J210 + Лэзо + J012) r2 + ( J201 + -Лш + Лэоз) Г3] +

+ [Jioorf + Joso^i + </003^'3 + 6Jllirir2r3 +

+ 3 (J210rfr2 + J20Ir^r3 + J120Г1 r| + J102Г1Г2 + J02irir3 + J012Г2Г|)] , 3

C4 = - [J400 + <Л)40 + <Л)04 + 2(J22O + ^202 + <Лш)] — 8

15 2 2 2

— [(^400 + J220 + ^202)^'l + (<Л)40 + J220 + Jo22)l'2 + (<Л)04 + ^202 + 3 +

+ 2 (( J310 + J112 + J130)r1 r2 + ( J301 + J121 + J103)r1 Гз + (J211 + J031 + J013)r2ГЗ)] + 35

+ [>hoor\ + Jo4or42 + Joo4?3 + Щ'Ьпг'1г2г3 + Ji2irir2r3 + Jm^i^rf) + + 4(J130 r1r3 + J 103r1 r3 + J301 Г3ГЗ + J310r3r2 + J031 Г|ГЗ + J013r3r2) + + 6( J220r^ + J202r^ + J022rlr3)] . The sequence of computation of the coefficients Jklk2k3 is presented below.

B. Step-by-step computation of coefficients Jklk2k3

Under the assumption of homogeneity of the body, the coefficients Jklk2k3 are calculated in the following order (see, e.g., [41, 42]).

The following order of computation is suggested:

1) computation of the body's mass: mB = J0 = J000;

2) determination of the center of mass Z of the body

OOZ = z = (z1, Z2, Z3)T = m-1(JW0, J010, J001 )T = m-1J1;

3) transfer from the coordinate frame OX1X2X3 to the coaxial coordinate frame: ZX1X2X3, with OXk y ZXk, k = 1, 2, 3;

4) computation of components J200, J020, J002, J110, J011, J101 of the second-order Euler-Poinsot tensor with respect to ZX1X2X3:

( J200 J110 J101 ^ J110 J020 J011 y J101 J011 J002 J

5) computation of the central tensor of inertia I = Tr(J2)E — J2;

6) computation of central principal moments I1, I2, I3 and unit eigenvectors ei, e2, e3 of the tensor I. These eigenvectors specify the principal central axes of inertia of the body;

7) transfer from the coordinate frame ZX1X2X3 to the coordinate frame Z^^ with the axes directed along the central principal axes of inertia;

8) computation of components of Euler-Poinsot tensors of the third and the fourth ranks: Jklk2k3, k1 + k2 + k3 = 3 and k1 + k2 + k3 = 4 with respect to the Z(1(2(3 coordinate frame.

C. Inertial characteristics of the Asteroid (2063) Bacchus

Determination of the inertial characteristics of some small celestial bodies is the subject of this section. Following the plan outlined in Section 3, the volume and position of the center of mass for celestial bodies under consideration are computed first. Then the matrix of the second rank Euler-Poisson tensor is computed to obtain the principal central moments of inertia for the body. The eigenvectors of the matrix define the principal axes of inertia for the body. Eventually, the components of the Euler-Poisson tensors of the third and fourth ranks J3 and J4 are computed in the reference frame composed by the principal axes of inertia.

To compute the inertial characteristics of the asteroid (2063) Bacchus, let us use a triangulation mesh from [43] defined by 2048 vertices and 4092 faces given in some coordinate system OX\X2X3.

The volume and the center of mass of the asteroid are

Vb = Jc/p = 0.1355 (km3), OOZ = J\/mB = (-0.02521593529, -0.001218836242, -0.001367422284)T (km).

The components of the tensor J2 divided by the mass of the asteroid presented in the axes ZX\X2X3 have the form (km2)

J2/mB =

0.0646393784920605 -0.0008902028705339 0.0002409649492319 -0.0008902028705339 0.0121264077275966 0.0002756872990353 0.0002409649492319 0.0002756872990353 0.0116427117309893

The eigenvectors of the matrix J2 determine the principal axes of inertia of the body. These vectors read

e1 = (0.999846899291092, -0.016920863811134, 0.004456719182762)T , e2 = (0.013453751646430, 0.906259358156108, 0.422507955334708)T , e3 = (-0.01118814303739, -0.42238330947420, 0.90634825830527)T .

In this case, the principal central moments of inertia Ie divided by the mass of the asteroid, denoted by ie, read (km2)

i1 = 0.023753, i2 = 0.076167, i3 = 0.076897.

The tensor components J3/mB and J4/mB, given in the axes Z(1(2(3, are written out in Tables 6 and 7, respectively.

Table 6. The components of the tensor J3/mB (km3)

J3oo = 0.002339975, J120 = -0.000148912, J102 = -0.000318486,

J012 = -0.00001405, J030 = -0.000007241, J02i = 0.000084352,

J20i = -0.000263668, J210 = 0.000429630, J003 = -0.000042338,

J111 =0.000121352.

Table 7. The components of the tensor J4/mB (km4)

J400 = 0.008333582, J103 = -0.000001884, J022 = 0.000096268,

J040 = 0.000328445, J301 = 0.000044096, J202 = 0.000617647,

J004 = 0.000288454, J130 = -0.000015269, J211 = 0.000010083,

J310 = 0.000142097, J013 = 0.000001549, J121 = 0.000001216,

J031 = -0.000005927, J220 = 0.000619338, J112 = -0.000005588