Asteroids and Their Mathematical Methods

Abstract

In this paper, the basic classification of asteroids and the history and current situation of asteroid exploration are introduced. Furthermore, some recent research progress on the orbital dynamics of asteroids, including models of the gravitational potential field, the dynamics near asteroids, hopping motion on the surface, and bifurcations under varying external parameters, is reviewed. In the meanwhile, the future research development such as the configuration and evolution of binary or triple asteroid systems and near-Earth asteroid defense is briefly discussed.

1. Introduction

Exploring, understanding, and trying to conquer one’s own environment is the most fundamental driver of human progress. During this process, human beings have discovered new knowledge, created new technologies, and discovered the truth through practice. Cognition guides practice, and this new knowledge and these new technologies actively develop from perceptual knowledge to rational knowledge, thus actively guiding practice, which in turn guides human beings to transform the world and promotes the further expansion of human territory. From Alexander’s expedition to the East to Zhang Qian’s mission to the Western Regions, Zheng He’s seven voyages to the West, and the space race in the middle of the 20th century, civilizations that have gained advantages in exploration have established scientific and technological advantages due to their exploration.

Although the exploration of travelers has been limited to land and sea due to technical reasons until modern times, the sages have never given up looking up at the stars. Moreover, the observation and understanding of the vast universe has constantly changed the human world view. Johannes Kepler derived Kepler’s laws using the detailed observational data of Tycho Brahe, making Nicholas Copernicus’s heliocentric theory recognized; Isaac Newton proposed the law of universal gravitation to explain the mathematical laws of celestial motion. In the three centuries after Newton and Kepler, the world’s greatest mathematicians studied celestial mechanics very well. In the decades before the first launch event of the Soviet Union in 1957, celestial mechanics did not necessary appear in college courses. However, before human spaceflight missions, celestial mechanics for ancient mathematicians, mechanics, and astronomers was limited to predicting the orbits of naturally existing celestial bodies in the solar system. Only in recent decades has the problem of orbital design for visiting target planets under complex constraints appeared [1].

With the advancement in science and technology, the scope of the world recognized by mankind in the past 200 years has also reached an unprecedented level. William Herschel discovered Uranus with the help of a telescope in 1781. This is the first time that humans discovered a planet with a telescope. Before this, the six major planets have been known to humans since ancient times. Twenty years later, Giuseppe Piazzi discovered the first asteroid, (1) Ceres, on 1 January 1801. It was initially thought to be a new planet and is now within the orbit of Neptune under the current definition. Since Piazzi interrupted the observation on 11 February 1801 due to illness, people only lost the asteroid with the observation data of Ceres’ orbit of about 3° in 41 days. Gauss used an orbit determination method he developed to determine Ceres’ orbit nearly a year after its disappearance using short-term observational data. Thanks to the work of Gauss, Franz Xavier von Zach rediscovered Ceres. Gauss’s calculation of Ceres clearly demonstrated that no assumptions were required and that the orbits of celestial bodies could be determined fairly accurately with just a few days of good observational data, followed by the discovery and orbit determination of (2) Pallas, (3) Juno, and (4) Vesta further verifies the efficiency of this method. This is also the first scientific progress made by humans based on small celestial bodies.

From the beginning of the 19th century to the middle of the 20th century, human beings gradually discovered more small objects in the solar system, but the related research was not given much attention because of their small masses. In the early 1950s, the development of intercontinental rockets allowed humans to use artificial spacecraft for spaceflight, and humans gradually began to plan space missions. Since the Soviet Union successfully launched Luna 1 in 1959, various major powers have made breakthroughs in near-Earth exploration and research on the Moon in the past 60 years of space exploration history. The main belt asteroids, the Jupiter system, the Saturn system, Pluto, the Kuiper Belt, and other solar system celestial bodies have been explored using various methods [2,3,4,5,6,7,8,9,10,11]. Especially since the 1970s, human beings have successively carried out many missions for small celestial bodies and have conducted in-depth and comprehensive research on small celestial bodies, both theoretically and practically. By studying the geological properties of celestial bodies and their space environment and exploring the formation and evolution history of the solar system, the chemical composition and internal structure of small celestial bodies have been preliminarily determined. The shapes of many small celestial bodies have been obtained through optical observations, radar observations, and photographs taken by fly-by missions. Among them, due to the irregular shape of small celestial bodies, the complexity of its nearby dynamic behavior, and its important research significance, the exploration of small celestial bodies has become an important part of deep space exploration, and countries have carried out space missions related to small celestial bodies. It has attracted a large number of scholars to study the dynamics of small celestial bodies [7,12,13,14,15]. It is widely accepted that small celestial bodies relatively completely retain the early information of the formation of the solar system [16,17,18]. Some small celestial bodies may also contain abundant rare metals and other resources needed by human beings, which have potential space mining value. The long-term effects of various perturbation forces in space and possible collisions between near-Earth celestial bodies may cause the small celestial bodies to change their original orbits and fall to the Earth, bringing devastating disasters to life. Some events are well-known—for example, the Tunguska explosion in 1909 and the meteorite impact in Chelyabinsk, Russia on 15 February 2013. Therefore, the observation and defense of near-Earth small objects is also a highly valued research subject [19,20]. In the past 50 years, the research on small celestial bodies in the solar system has already broken through the content of past astronomy, involving the development of nonlinear dynamics, including aerospace engineering, geology, biology, weapons research, and other fields.

For deep space explorations, it is inseparable from the content of close-range fly-by, imaging, landing, sampling, and returning. Compared with pure radar observations, these on-site explorations can provide more intuitive and richer details and evidence for scientific research on small celestial bodies and the defense of near-Earth small objects. The resource exploitation of small celestial bodies is even more inseparable from actual field detection. Some important missions that achieved sample returns from celestial bodies other than the Moon are Stardust, Hayabusa, Hayabusa 2, and OSIRIS-REx. One of the landmarks made by China to become one of the top aerospace powers around 2030 is to realize the sampling return of small celestial bodies.

In-depth research on the shape of small celestial bodies, nearby periodic orbital motions, and surface transition dynamics and quantitative analysis methods for quasi-periodic orbits is the scientific basis for completing the above tasks. In the process of close-range explorations, the orbits of probes near small celestial bodies are significantly different from those near large planets, and the internal mechanisms and laws have not been completely sorted out so far. What is certain is that the regularity of the geometric shape of the small celestial body largely determines the difference between the gravitational field near the small celestial body and the familiar spherical gravitational field near the large planet, which brings about a series of different dynamics and control problems. Therefore, exploring the indicators for scientifically describing the geometric shape of small celestial bodies has a great reference value for the preliminary analysis of the gravitational fields. In addition, considering that the mass of small celestial bodies is generally small, many space perturbations, including the Sun’s gravity, may also have a relatively large influence on the motion of the probe. The richness of the orbital dynamics near small celestial bodies also provides a good soil for discovering and studying nonlinear dynamics in real systems and making them useful. The difficulties of the orbital dynamics and the limitations of the measurement methods make the data parameters obtained before approaching the target small celestial body still have a certain error with the real situation. These factors and the complex mechanical environment near the small celestial body jointly affect the exploration. The orbital design and control of orbiters near small bodies pose great challenges. Even for periodic orbits, considering various gravitational and perturbative effects, they may be perturbed into quasi-periodic orbits or produce chaotic motions. A thorough and accurate analysis of these orbits will aid in orbital design for exploration missions. Due to the strong irregular shape, complex topography, and fast rotation rate of small celestial bodies, the surface escape velocity of small celestial bodies is usually small. To truly realize the sampling return of small celestial bodies, it is necessary to select the soft-landing region and optimize the trajectory, which requires us to deeply study the dynamic laws of surface motions on small celestial bodies.

Although humans have now begun to explore the solar system, the small celestial bodies that account for the vast majority of solar system celestial bodies have only been visited by very few spacecraft. The continued exploration of the vast universe by mankind in the future is inseparable from the in-depth detection and research of small celestial bodies. By studying the shape regularity and nearby periodic orbits of small celestial bodies, the dynamic phenomena in the theory can be verified in the real system in science, enriching the scientific connotation of modern celestial mechanics and nonlinear dynamics [21,22].

2. Basic Classification and Exploration of Small Celestial Bodies

Since the 1980s, many deep space exploration missions related to small celestial bodies have been carried out. This section will introduce the general situation of small celestial bodies, deep space exploration missions, etc.

2.1. Overview of Small Celestial Bodies

Small celestial bodies usually refer to celestial bodies that exclude large planets and their satellites in the solar system, mainly including dwarf planets, asteroids, and comets [23,24]. Small celestial bodies orbit the Sun but are much smaller in size and mass than large planets. According to data from the IAU Minor Planet Center, as of 10 July 2021, a total of 1,091,258 small celestial bodies have been discovered in the solar system, of which there are 5 dwarf planets and 4603 comets, and the rest are asteroids [25]. Among these asteroids, 567,132 have been permanently numbered (orbits have been calculated) (more than 90% of them have been newly discovered in the past 20 years; see Figure 1) [26]. A total of 22,568 asteroids have been named [27].

The semi-major axis of near-Earth asteroids is similar to that of Earth. There are currently 26,351 near-Earth asteroids, of which 940 are larger than 1 km in size. Scientists believe that the mass extinction 65 million years ago was caused by an asteroid about 10 km in size hitting the Earth [20]. The defense of near-Earth asteroids is also an important part of the field of deep space exploration. The orbital inclinations of near-Earth asteroids range from 0.02° to 154°, and the orbital eccentricity ranges from 0.062 to 0.999 [28,29]. According to its orbital semi-major axis a, perihelion distance q, and the relationship between the aphelion distance Q and the Earth, it can be divided into the Atira type, Aten type, Apollo type, and Amor type, as shown in

Table 1. Classification of Near-Earth Asteroids.

Type	Semi-Major Axis a	Perihelion Distance q	Aphelion Distance Q
Atira	a < 1 AU	-	Q < 0.983 AU
Aten	a < 1 AU	-	Q > 0.983 AU
Apollo	a > 1 AU	q < 1.017 AU	-
Amor	a > 1 AU	1.017 AU < q < 1.3 AU	-

.

Atira-type and Amor-type asteroids are less dangerous to Earth because their orbits and Earth’s orbits are inward and outward, respectively. Arten-type and Apollo-type asteroids are small celestial bodies with a greater potential danger to the Earth due to their inward and outward swept orbits and Earth’s orbits, respectively.

Near-Mars asteroids are divided into Hungarian-type asteroids and Mars-orbiting asteroids. The semi-major axes of the Hungarian-type asteroids are between 1.78 and 2.00 AU, and they are located inside the Kirkwood gap, which is in a 1:4 resonance with Jupiter. Their orbital periods are about 2.5 years, roughly 3:2 resonant with Mars and 2:9 resonant with Jupiter. Their orbital eccentricities are less than 0.18, and their orbital inclinations are between 16° and 34°. The perihelion distance of the Mars orbit crossing the asteroid orbit is between the distance between the perihelion and the aphelion of Mars, that is, 1.381 AU < q < 1.666 AU. Asteroids with a perihelion distance q < 1.3 AU are classified as near-Earth asteroids. According to this classification standard, a total of 18,043 Mars-orbiting asteroids have been discovered so far.

Asteroids in the middle solar system are divided into main belt asteroids, Jupiter Trojan asteroids, and Hilda asteroids. The main-belt asteroids are located between the orbits of Mars and Jupiter, and the orbital semi-major axis a is between 2.1 and 3.3 AU. The orbital eccentricities of most main-belt asteroids are less than 0.4, and the orbital inclination angles are less than 30°. The main belt is the area with the densest distribution of asteroids, and a total of 1,022,771 asteroids have been observed in the main belt. It is generally believed that the main-belt asteroids are the remnants of the original astrolabe that failed to form large planets due to the perturbation of Jupiter’s huge gravitational force during the evolution of the solar system. According to the current definition, the three largest asteroids in the main belt are small celestial bodies: (4) Vesta, (2) Pallas, and (10) Hygiea. The Jupiter Trojan asteroids are located near the L4 and L5 points of the circular restricted three-body system. The L4 and L5 points are regarded as the stable equilibrium points of this system. The motion periods are basically the same as that of Jupiter, with a phase difference of about 60°. So far, 10,470 Jupiter Trojan-type asteroids have been observed. In addition to the Sun–Jupiter system, there are also four and six Trojan asteroids in the Sun–Mars and Sun–Neptune systems, respectively. So far, only one Trojan asteroid, 2010 TK7, has been discovered in the Sun–Earth system; this was done in 2010. It is located near the L4 point of the Sun–Earth system. A total of 4978 Hilda-type asteroids have been observed. Their semi-major axes are between 3.7 and 4.2 AU. Their orbital eccentricities are less than 0.3, and their orbital inclinations satisfy i < 20°. The Hilda-type asteroids are in a 2:3 resonance with the orbit of Jupiter and approach the L3, L5, and L4 points of the Sun–Jupiter system, in turn, in three orbital periods.

Asteroids in the outer solar system include Centaur and extra-Neptunian asteroids. Centaurs are small celestial bodies whose perihelions are outside the orbit of Jupiter and whose semi-major axes are smaller than Neptune’s semi-major axis by 30 AU. Because the small celestial bodies here have the characteristics of asteroids and comets, they are mostly named after the centaur gods in Greek mythology. For example, (2060) Chiron and (60588) Echeclus have the comet numbers 95P/Chiron and 174P/Echeclus due to coma activity. Extra-Neptunian asteroids refer to the celestial bodies in the solar system whose semi-major axes are greater than 30 AU. Excluding the currently discovered (134340) Pluto, (136108) Haumea, (136472) Makemake, and (136199) Eris and four other dwarf planets, a total of 4,053 extra-Neptunian asteroids have been discovered so far. Most of these celestial bodies contain methane, ammonia, and water, which are volatile. The region outside Neptune between 30 and 50 AU from the Sun is called the Kuiper Belt. Similar to the main-belt asteroids, Kuiper Belt objects are also the original remnants that failed to form large planets. The interesting relations between Kuiper Belt objects and comets can be found in reference [30]. In addition, research on the orbital dynamics of Kuiper Belt objects plays an important role in the process of the human search for the ninth largest planet in the solar system. Based on the orbital eccentricity vector and angular momentum vector of six Kuiper Belt objects, Batygin [31] and Brown [32] inferred that there may be an unknown planet with a semi-major axis a ≈ 700 AU and an eccentricity e ≈ 0.6. This study has prompted scholars from various countries to conduct in-depth research and sky survey observations to find this potential ninth planet in the solar system [33].

According to the spectral characteristics, there are 17 types of asteroids: A, B, C, D, E, F, G, K, L, M, O, P, Q, R, S, T, V. Asteroids are mainly divided into three groups: C, S, and X [34,35,36], and a few other types are not included in these three groups. Group C asteroids contain a large amount of carbon, accounting for about 75% of the total number of asteroids in the solar system; group S asteroids contain a large number of silicates, accounting for about 17% of the total number; group M asteroids contain a large amount of iron and nickel and other metallic elements, which are considered to be the debris of the asteroid’s impacted core and the source of iron meteorites [37]. The specific classification of asteroids according to their spectra can be found in

Table 2. Spectral Classification of Asteroids.

Group	Type	Criterion	Representative
C	B	The general properties are the same as the C type, but the ultraviolet absorption below 0.5 μm is smaller, and the slight blueness is more obvious than the redness in the spectrum. The albedo also tends to be greater than the darker C-type.	(2) Pallas
	C	There is moderate absorption at UV wavelengths of 0.4–0.5 μm, and there are no obvious features but slight reddening at longer wavelengths. There is a mineral feature indicative of hydration known as water absorption around the wavelength of 3 μm.	(10) Hygiea
	F	Similar to B-type asteroids but lacks water absorption features indicative of hydrated minerals around wavelengths around 3 μm and differs from the B-type in the low-wavelength UV portion below 0.4 μm.	(704) Interamnia
	G	Similar to C-type asteroids but has strong absorption characteristics for ultraviolet wavelengths below 0.5 μm. There may also be absorption properties around 0.7 μm, implying the presence of layered silicate minerals such as clay and mica.	(1) Ceres
S	A	Significant olivine features at a 1 μm wavelength and strong reddening at wavelengths below 0.7 μm.	(446) Aeternitas
	K	There is moderate reddening at wavelengths below 0.75 μm and slight bluing at wavelengths above 0.75 μm.	(221) Eos
	L	There is strong reddening at wavelengths below 0.75 μm, and the spectrum is flat at wavelengths above 0.75 μm. Compared with the K type, the redness is more obvious in the visible band, and the spectrum in the infrared band is more gentle.	(83) Beatrix
	Q	There are prominent features of olivine and pyroxene in the 1 μm band, and their spectral changes indicate the possible presence of metallic substances. There is an absorption spectrum at 0.7 μm.	(1862) Apollo
	R	There are distinct olivine and pyroxene features at 1 μm and 2 μm. The spectrum is strongly reddened at wavelengths below 0.7 μm.	(349) Dembowska
	S	There is moderate spectral variation at wavelengths shorter than 0.7 μm and moderate spectral absorption at 1 μm and 2 μm wavelengths. There is also a shallow but broad spectral absorption around 0.63 µm.	(3) Juno
X	E	The albedo is greater than 0.3, the spectrum is flat and reddish, and there are no obvious features.	(44) Nysa
	M	The albedo is between 0.1 and 0.2, there are subtle spectral absorption lines in the bands above 0.75 μm and below 0.55 μm, and the overall spectrum is flat and slightly reddened, lacking obvious features.	(16) Psyche
	P	The albedo is less than 0.1, and the color is redder than that of the S-type asteroid, but it is not reflected in the spectral properties.	Sylvia
Not grouped	D	Very low albedo and featureless, light red electromagnetic spectrum.	(624) Hektor
	O	Strong spectral absorption in the band above 0.75 μm	(3628) Božněmcová
	T	The spectrum is moderately reddened, darker, and has moderate spectral absorption in the band below 0.85 μm.	(114) Kassandra
	V	There is strong spectral absorption in the bands above 0.75 μm and 1 μm and strong reddening in the bands below 0.7 μm.	(4) Vesta

. It should be noted that the diversity of the spectra of asteroids can be affected by space weathering [38,39].

Comets can be divided into the nucleus, coma, and tail. Comet nuclei are composed of loose water ice, rubble piles, solid carbon dioxide, methane, ammonia, etc. [40] Comets usually have long-period, highly eccentric orbits. Therefore, as the comet approaches the Sun, the water ice and volatile matter in the comet’s nucleus will be heated and turned into gas, forming an observable atmosphere called a coma. The coma is affected by the solar wind and solar light pressure to produce a long tail facing away from the Sun, called a comet tail. Comets are perturbed by the gravitational force of Jupiter and other large planets during their operation, which may cause dramatic changes in their orbits or their own shapes or even disintegration: Comet Shoemaker-Levy 9 (Shoemaker-Levy 9) was destroyed by Jupiter in 1994. The gravity ripped apart into 21 pieces and crashed into Jupiter. Hsieh and Jewitt [41] inferred the existence of a population of comets originating in the main asteroid belt based on the optical data. In the past, people often used the presence of volatile gas as a criterion to distinguish asteroids and comets, but with the discovery of active centaurs, especially after the discovery that Ceres also has water vapor [42], the difference between asteroids and comets becomes less clear.

2.2. Exploration of Small Celestial Bodies

With the development of science and technology, the exploration of small celestial bodies has gradually developed from optical observation with telescopes in the 19th century to radar observation and visits through spacecraft. These three methods provide different kinds of information for the research and exploration of small celestial bodies.

Since Piazzi discovered Ceres, the vast majority of small objects have been discovered through optical observations. The introduction of astrophotography and scintillation comparators allowed optical observations to move away from relying on the naked eye to identify asteroids. Using advanced orbiting telescopes and observatories, it has been possible to obtain basic images of large-scale small celestial bodies through optical observation. The double asteroid system (45) Eugenia was discovered in this way [42]. Orbit information can also be used to calculate the size of small celestial bodies by observing the apparent magnitude. Moreover, the rotation period and spatial orientation of the rotation axis can be calculated from the light change information of small celestial bodies. The temperature and spectral information of asteroids can be determined through optical observations in visible light and infrared bands [39]. Different from the optical observation, radar observation is an active observation method, and radar observation can provide the orbit data of small celestial bodies with a higher relative accuracy and information such as small celestial body shape, rotation speed, and albedo. A higher-precision model of small celestial bodies (on the order of 10 m) can also be reconstructed through radar observations. Since the first high-precision shape model of (4769) Castalia was reconstructed in 1994, more and more small celestial body models have been obtained by this method. However, due to the attenuation of radar echoes, ideal radar observations require small celestial bodies to be close enough to Earth, so radar observations are mostly concentrated in near-Earth asteroids.

With the deepening of deep space exploration activities, the United States, the Soviet Union, and Europe have carried out space exploration activities for small celestial bodies since the 1980s. Japan and China have joined the deep space exploration team one after another.

Table 3. Asteroid exploration missions.

Spacecraft	Agency	Start Date	Asteroid	Mission Type
International Comet Explorer	NASA ESA	1982	21P/Giacobini-Zinner	Fly-by
Vega 1/2	IKI	1984	1P/Halley	Fly-by
Pioneer/Comet	JAXA	1985	1P/Halley	Fly-by
Giotto	ESA	1985	1P/Halley 26P/Grigg–Skjellerup	Fly-by
Galileo	NASA	1989	(951) Gaspra (243) Ida	Fly-by
Near-Shoemaker	NASA	1996	(253) Mathilde (433) Eros	Fly-by/Orbiting/Landing
Cassini-Huygens	NASA	1997	(2685) Masursky	Fly-by
Deep Space 1	NASA	1998	(9969) Braille (19P/Borrelly)	Fly-by/Orbiting
Stardust	NASA	1999	(5535) Annefrank 81P/Wild 2 9P/Tempel 1	Fly-by
Comet Nucleus Tourer (Failed)	NASA	2002	2P/Encke 73P/Schwassmann-Wachmann 3 6P/d’Arrest	Fly-by
Hayabusa	JAXA	2003	(25143) Itokawa	Orbiting/Landing/Sample return
Rosetta	ESA	2004	(2867) Steins (21) Lutetia 67P/Churyumov-Gerasimenko	Orbiting/Landing/Sample return
Deep Impact/EPOXI	NASA	2005	9P/Tempel 1 103P/Hartley 2	Impact/ Fly-by
New Horizons	NASA	2006	132524 APL (134340) Pluto 2014 MU69	Fly-by
Dawn	NASA	2007	(4) Vesta (1) Ceres	Orbiting
Chang’e 2	CNSA	2010	(4179) Toutatis	Fly-by
Hayabusa 2	JAXA	2014	(162173) Ryugu	Orbiting/Landing/Sample return
OSIRIS-REx	NASA	2016	(101955) Bennu	Orbiting/Landing/Sample return
Don Quixote (in progress)	ESA	-	2003 SM84	Fly-by/Impact
Double Asteroid Redirection Test (DART)	NASA	2021	(65803) Didymos	Impact
Lucy	NASA	2021	15094 Polymele 21900 Orus	Fly-by
Psyche mission	NASA	2022	(16) Psyche	Orbiting
Tianwen 2 (in progress)	CNSA	2025	(469219) Kamo‘oalewa (2016 HO3) 311P/PanSTARRS	Orbiting/Landing/Sample return Orbiting

lists the small object missions that have occurred and may be carried out in the future. Through these exploration missions, human beings have gained further understanding of the geological characteristics of small celestial bodies, the space environment, and the formation and evolution of the solar system.

The early exploration activities of small celestial bodies were mainly affected by the return of Halley’s Comet (1P/Halley) in 1986, and the fly-by of the comet was the mainstay. The first exploration of small celestial bodies by humans was the International Cometary Explorer (ICE), jointly conducted by ESA and NASA in 1982. The predecessor of ICE was the first International Sun–Earth Explorer-3 (ISEE-3) located at the Sun–Earth L1 point. It was renamed as the International Comet Explorer to conduct comet exploration activities. After a low-altitude fly-by of the Moon on 22 December 1983 for gravity assistance, the International Comet Explorer passed through the tail of Comet 21P/Giacobini-Zinner at a distance of 7800 km from the nucleus in 1985, while the geomagnetic field downstream of the long tail blown by the solar wind was also detected during the Earth–Moon gravitational assistance and passed through the tail of Halley’s Comet in 1986 [43,44,45]. From 1984 to 1985, when Halley’s Comet returned, the Soviet Union launched Vega-1 and Vega-2 successively. In the process of exploring Venus, the two spacecraft were placed at distances of 10,000 km and 3000 km to conduct fly-bys of Halley’s Comet. Japan also launched two probes, Sakigake and Suisei, to conduct fly-bys of Halley’s Comet at distances of 7,000,000 km and 150,000 km. In 1985, ESA launched the Giotto probe to observe Halley’s Comet. Giotto flew by Halley’s Comet at a distance of 596 km in March 1986 and was the first probe to observe the comet at close range [46,47,48]. Giotto flew by Comet 26P/Grigg-Skjellerup at a distance of 200 km in 1992 after gravitational assistance in 1990. Vega 1, Vega 2, Pioneer, Comet, and Giotto are known as the “Halley Fleet” for their continuous exploration of Halley’s Comet. In 1989, NASA launched the Galileo probe on its way to Jupiter. It flew by (951) Gaspra and (243) Ida in 1991 and 1993, respectively, and discovered the moons (Dactyl) of Ida. This is the first time that humans have explored asteroids and double asteroid systems [49,50,51,52].

In the 1990s, the exploration of small celestial bodies took various forms such as orbiting, impacting, landing, and sampling return. In 1996, NASA launched the Rendezvous–Shoemaker probe. Shoemaker was originally only planned for an orbital mission, but Dr. Bob Farquhar calculated that Shoemaker could successfully land in a saddle-shaped area on the southern surface of the small body (433) Eros without a soft-landing device. Shoemaker landed undamaged and continued to work for another 16 days, making it the first probe to soft-land an asteroid [53,54,55,56,57]. In 1997, NASA launched the Cassini-Huygens probe to probe Saturn. Cassini flew by (2685) Masursky on its way to Saturn in 2000, confirming that its diameter is between 15 and 20 km [58]. In 1998, the Deep Space 1 probe launched by NASA first flew by (9969) Braille in July 1999 and then rendezvoused with Comet Borrelly (19P/Borrelly) in September 2001 for observations [59,60,61,62]. In 1999, NASA launched the Stardust probe, which flew by (5535) Annefrank in 2002-11 and flew by 81P/Wild 2 in January 2004. The dust was sampled and returned. In February 2011, the probe visited 9P/Tempel 1 [63,64,65,66,67,68,69,70,71].

After entering the 21st century, in 2002, NASA launched the Comet Nucleus Tourer (CONTOUR), which plans to fly by Comet 2P/Encke, 73P/Schwassmann-Wachmann 3, and 6P/d’Arrest. Due to the failed launch, the mission has become the only small celestial object mission that has completely failed so far. In 2003, the Japan Aerospace Exploration Agency (JAXA) launched the Hayabusa probe to visit (25143) Itokawa. It landed on the asteroid in November 2005, collected some asteroid samples, and returned to the Earth. It is the first deep space exploration mission to sample and return asteroids [72,73,74,75,76]. In 2004, ESA launched the Rosetta comet probe to visit 67P/Churyumov-Gerasimenko, and on the way, it flew by (2867) Steins at a distance of 800 km in 2008 and flew over (21) Lutetia at a distance of 3160 km in 2010 [77]. On 12 November 2014, the lander Philae carried by the Rosetta probe landed in the pre-selected J region on the comet, becoming the first probe to land on the surface of a comet’s nucleus [78,79,80,81]. In 2005, NASA launched the Deep Impact program to study the composition of the comet nucleus of Comet Tempel 1. In July of the same year, the impactor was released to complete the impact mission of Comet Tempel 1 and probe in deep space. For the first time in history, the material ejected from the surface of a comet had been measured [82,83,84]. The deep impact was then extended to the EPOXI mission, which conducted a fly-by of Comet 103P/Hartley 2 in November 2010 with the aid of Earth’s gravity [85,86]. In 2006, NASA launched the New Horizon probe, which flew by 132,524 APL at a distance of 100,000 km in 2006 and approached the dwarf planet Pluto and its five moons in January 2015, becoming the first probe in history to visit a dwarf planet. It flew by Pluto at a distance of 12,500 km in July of the same year, after which NASA designated the fly-by of the Kuiper Belt small object 2014 MU69 as an extended mission of New Horizons [87]. In 2007, NASA launched the Dawn probe to visit (4) Vesta and the dwarf planet Ceres. Dawn arrived at the small celestial body (4) Vesta in July 2011. After the exploration, it flew to Ceres and arrived at Ceres in March 2015. This is the first probe in human history to orbit asteroids and dwarf planets in the main belt [88,89,90,91]. In 2010, the China National Space Administration (CNSA) launched Chang’e-2 probe to orbit the Moon. After that, the probe went to (4179) Toutatis in April 2012 and completed the mission of flying over (4179) Toutatis at a distance of 3.2 km in December 2012, which led to the obtention of a clear image of the surface of (4179) Toutatis for the first time [92,93,94]. Through this extended mission, China explored an asteroid for the first time, becoming the fourth country in the world to explore asteroids after the United States, Europe, and Japan. In 2014, after the Hayabusa mission, JAXA launched the Hayabusa-2 probe to detect (162173) Ryugu and used the blasting method to collect its deep samples and return. On 6 December 2020, the “Hyabusa 2” landed in the desert area of southern Australia and obtained 5.4 g of the sample, which aroused enthusiastic attention from all walks of life. In 2016, NASA launched the Pluto probe (OSIRIS-REx) to carry out a sampling return mission to (101955) Bennu. In December 2018, the probe reached the asteroid Bennu. After more than a year of short-range detection, a “touch and go” sampling was carried out to confirm that the sample was collected, and it is planned to arrive on Earth in September 2023.

In order to effectively deal with the potential collision threat of NEOs to Earth, scientists have been studying various means of NEO defense. For asteroid defense, the basic technical way is to use nuclear bombs for interception, use spacecraft for kinetic energy impact, or use laser ablation and other schemes. Other methods such as ion traction, gravitational drag, and mass drive are still in the argumentation stage. ESA launched a preliminary study of the Don Quijote project in 2006 and plans to test asteroid defense technology targeting 2003 SM84 or Destroyer (99942) Apophis in the future. NASA started the Double Asteroid Redirection Test (DART) in November 2021 to test the asteroid impact and defense against the binary asteroid system (65803) Didymos. According to the plan, the DART spacecraft will approach the target asteroid between the end of September and the beginning of October 2022 and finally hit the asteroid head-on at a speed of about 6.6 km/s. Such an impact would change the speed of the Moon in its orbit around the primary by one percent, and the period may be varied by a few minutes. China has been monitoring near-Earth asteroids since the 1950s. At the 7th International Academy of Astronautics Planetary Defense Conference, Yanhua Wu, Deputy Director of the Chinese Space Administration, mentioned that China will establish an organizational system and process to deal with the risk of asteroid impact.

In addition to planetary defense, research on the evolution of the solar system is also a focus related to small celestial bodies. NASA launched the Lucy probe in October 2021 and planned to fly by an asteroid in the inner main belt and five other Jupiter Trojan asteroids. NASA also plans to launch the Psyche spacecraft in 2022 to probe a series of questions related to the evolution of Psyche through its orbit.

The white paper “2016 China’s Aerospace” proposes to deepen the demonstration and key technological breakthroughs in the main tasks of Mars sampling and return, asteroid exploration, Jupiter system and planetary transit detection, etc. In April 2019, the “Asteroid Exploration Mission Payload and Carrying Project Opportunities Announcement” issued by the National Space Administration confirmed that China’s asteroid exploration mission will achieve a near-Earth asteroid sampling return and a main-belt comet orbit through one launch. In 2021, the first Mars exploration mission of China (Tianwen 1) has achieved great success. China plans to launch a small celestial body detector around 2025 and spend 10 years visiting two small celestial bodies: (2016) HO3 near-Earth asteroid and comet (133) P.

3. Research on Small Celestial Body Gravitational Field Environments and Orbital Mechanics

Deep space exploration missions related to small celestial bodies have stimulated a large number of studies on related issues, including analysis of the dynamic environment of small celestial body detectors, orbit design and control, the formation and evolution of celestial bodies, nonlinear dynamic characteristics, and other aspects. This section will introduce the research on the gravitational field model of small celestial bodies, the orbital dynamics, the surface transition dynamics, and the research on the gravitational field environment of small celestial bodies under the change of physical parameters, etc.

3.1. Research on the Gravitational Field Model of Small Celestial Bodies

The study of the periodic and quasi-periodic orbital dynamics around small irregular celestial bodies depends on the proper description of the dynamic model of the small celestial bodies. The work of Hamilton et al. [95] shows that the influence of the planetary gravitational perturbation is very small compared with the Sun’s gravitational perturbation, and its influence can be ignored. Therefore, the dynamic equation of the particle moving near the small celestial body can be expressed as

$$\stackrel{··}{\mathbf{r}}={ \mathbf{a}+\mathbf{a}}_S$$

(1)

In the formula, $r$ represents the position of the particle in the small celestial body coordinate system, $a$ represents the acceleration obtained by the particle due to the gravitational force of the small celestial body, and ${\mathbf{a}}_S$ represents the acceleration obtained by the particle perturbed by the gravitational force of the Sun. Since the particle motion is considered, the influence of solar pressure perturbation can be ignored. For the distance range applicable to the dynamic model of small celestial bodies, generally consider the radius of the sphere affected by the gravitational force of the small celestial body relative to the gravitational force of the Sun:

$$\frac{R_1}{D}={\left(\frac{M_A}{M_S}\right)}^{\frac{2}{5}}$$

(2)

In the formula, $R_1$ is the radius of the influence sphere, $D$ is the revolving distance of the small celestial body around the Sun, $M_A$ is the mass of the small celestial body, and $M_S$ is the mass of the Sun. Within the range of the small celestial body’s influence sphere radius $R_1$, the gravitational force of the small celestial body can be regarded as the main force affecting the motion of the particle, and the Sun’s gravitational force can be regarded as the perturbing force. For a more rigorous estimate, the radius at the point of gravitational neutralization for small bodies and the Sun can be considered:

$$\frac{R_2}{D}={\left(\frac{M_A}{M_S}\right)}^{\frac{1}{2}}$$

(3)

In the formula, $R_2$ is the radius of the gravitational neutralization point. It is easy to see that R₁ > R₂. Yu [96] gave the gravitational radii of 23 small celestial bodies in his doctoral dissertation.

Table 4. Gravitational radius of asteroids.

Asteroids	MA/MS	D/AU	R1/km	R2/km
(216) Kleopatra	2.33 × 10−12	[2.09, 3.49]	[6.97 × 103, 1.16 × 104]	[4.78 × 102, 7.99 × 102]
(243) Ida	2.11 × 10−14	[2.74, 2.98]	[1.39 × 103, 1.51 × 103]	[5.97 × 10, 6.50 × 10]
(433) Eros	3.36 × 10−15	[1.13, 1.78]	[2.75 × 102, 4.34 × 102]	[9.83 × 100, 1.55 × 10]
(1620) Geographos	1.30 × 10−17	[0.83, 1.66]	[2.19 × 10, 4.38 × 10]	[4.49 × 10−1, 8.98 × 10−1]
(6489) Golevka	1.06 × 10−19	[0.99, 4.02]	[3.81 × 100, 1.55 × 10]	[4.83 × 10−2, 1.96 × 10−1]

lists the radius range of the small celestial body’s gravitational action, calculated according to Formulas (2) and (3). It is not difficult to see that the radius $R_1$ of the gravitational influence sphere is about two orders of magnitude larger than the radius $R_2$ at the gravitational neutralization point.

Small celestial bodies are irregular in shape, rotate faster than large planets, and have smaller masses. They are very different from the gravitational fields around large celestial bodies, showing the characteristics of asymmetry and irregularity. Therefore, in order to study the dynamics of small irregular celestial bodies, it is necessary to approximate their gravitational field with an appropriate model. Common approximate models of gravitational field are: the simple geometry model; the spherical harmonic and ellipsoidal harmonic function model; the particle group model; and the polyhedron model.

The simple geometry model has the characteristics of a simple structure, few shape parameters, and convenient calculation. It is easy to obtain analytical results and qualitative conclusions about shape parameters. The simple geometries commonly used for simulation in current research include the homogeneous thin straight rod model [97], the homogeneous ring model [98], the triaxial ellipsoid model [99], and the dipole model [100]. The early simple geometric models can only reflect the basic characteristics of small celestial bodies and cannot accurately simulate the surrounding gravitational field environment. With the continuous development of the study of simulating irregular gravitational fields with simple geometric models, the dipole model proposed by Zeng et al. [100] can better reflect the gravitational field near the equilibrium point of irregular small celestial bodies. Accuracy is also taken into account.

The main idea of the spherical harmonic model is to use infinite series to approximate the gravitational potential function of celestial bodies. It was first applied in the dynamics of near-Earth satellite orbits and then introduced into the study of gravitational field modeling of small celestial bodies to describe small celestial bodies. Using the spherical harmonic function method, the gravitational potential can be expanded as

$$U(\mathbf{r})=\frac{G{M}_A}{r}\{1+{\displaystyle \sum _{l=1}^{\infty }{\displaystyle \sum _{m=0}^l{(\frac{r_e}{r})}^l}}{P}_{lm}(\sin (\varphi ))[C_{lm}\cos (m\lambda )+S_{lm}\sin (m\lambda )]\}$$

(4)

where G is the gravitational constant G = 6.674 28 × 10⁻¹¹ m³·kg⁻¹·s⁻²; $r$ is the position vector of the particle; $r,\varphi ,\lambda$ are the three components of the vector in spherical coordinates; M_A is the mass of the small celestial body; P_lm is the associative Legendre polynomial; r_e is the radius of the Brillouin sphere, reflecting the range of convergence of the series, that is, the applicable range of Formula (4); and C_lm and S_lm are spherical harmonic coefficients, reflecting the shape irregularity and inhomogeneity of internal mass distribution [101]. The advantage of this method is that the gravitational potential can be given analytically. It is convenient for obtaining theoretical solutions through analytical methods. In addition, once the spherical harmonic function coefficient is obtained, it can be directly substituted in the subsequent numerical calculation, which is convenient to use, especially for inversion calculation through flight data [102]. During the orbiting of Ceres by the Dawn probe, Takahashi et al. [103,104] used the spherical harmonic model to estimate the precise gravitational field of Ceres and iteratively iterated the known spherical harmonics to give the direction of its principal axis. The main limitations are that the model cannot be applied to the region located within the Brillouin sphere. The reason for this is that the series does not converge [105]. The truncation error in calculations may lead to large errors in the obtained gravitational field model in some cases [106].

Considering the convergence region of spherical harmonics, Hobson [107] uses Lamé polynomials to approximate the ellipsoidal harmonics model of the gravitational potential function. Pick [108] established the theory of ellipsoid harmonics on this basis. Using the ellipsoid harmonic method, the gravitational potential of a unit mass particle can be expanded as

$$U(r)=G{M}_A{\displaystyle \sum _{l=0}^{\infty }{\displaystyle \sum _{m=1}^{2l+1}{\alpha }_{lm}}}\frac{F_{lm}({\lambda }_1)}{F_{lm}({\lambda }_e)}{F}_{lm}({\lambda }_2)F_{lm}({\lambda }_3)$$

(5)

In the formula, λ₁, λ₂, λ₃ are the ellipsoid coordinate components of the vector $r$. λ_e is the parameter of the Brillouin ellipsoid, which reflects the range of convergence of the series, that is, the range of use of Formula (5). F_lm is the Lamé equation canonical solution, and α_lm is the ellipsoid harmonic coefficient [105]. A conversion method between spherical harmonics and ellipsoidal harmonics was proposed by Dechambre et al. [109] to simplify the solution process. The ellipsoid harmonic function model expands the convergence region of small celestial bodies while still retaining the characteristics of the spherical harmonic function model for easy calculation [110].

The convergence rate of the two models near the boundary of the convergence domain decreases rapidly as the distance from the particle to the small body decreases. In addition, the spherical harmonic function and ellipsoidal harmonic function models lack the information to judge whether the particle is located inside or outside the small irregular celestial body. Therefore, it cannot meet the calculation requirements of the global gravitational field in the application of studying the dynamics.

Particle swarm models are often used to model dynamic environments in asteroid evolution and near-Earth asteroid orbit collision avoidance problems. This model is a very intuitive method that discretizes the space where the small celestial body is located into a series of particles, calculating the gravitational force or gravitational potential of these particles separately and summing them up to obtain the overall gravitational force or gravitational potential of the small celestial body. Assuming that the small celestial body is divided into N voxels, the position coordinate of the i-th voxel is $r_i$, and the mass is M_i. Then, the gravitational potential can be expressed as

$$U(r)={\displaystyle \sum _{i=1}^N\frac{G{M}_i}{|r-r_i|}}$$

(6)

The advantage of the method is that its algorithm is simple and easy to implement, and it can ensure the convergence. By increasing the number and distribution of divided voxels by appropriate rules, the accuracy of the gravitational field of small celestial bodies can be improved. Real-world problems such as distribution have good scalability [111,112]. However, this method also has some defects: the number of voxels increases rapidly with the accuracy requirements, which leads to a sharp increase in the amount of calculation and a great decrease in the calculation speed; it cannot provide a direct and effective judgment for whether the orbit of the particle intersects with the small irregular celestial body.

The polyhedron model is a numerical modeling method for the gravitational field of small irregular celestial bodies. Since the nineteenth century, in order to describe a rugged terrain in geology, scholars have studied the gravitational field of a simple polyhedron. MacMillan and Waldvogel [113,114] successively gave the analytical form of the gravitational potential energy of the cuboid and the general homogeneous polyhedron. The disadvantage is that the amount of calculation is large. In the 1990s, Werner [115,116] used a polyhedron in which every surface is a triangle to approximate the shape of a small irregular celestial body. Moreover, he used Gauss’s theorem and Green’s formula to simplify the triple integral. Werner also studied the orbital behavior around the regular tetrahedron, which was compared to the orbitals affected by the J₃ and J₃₃ terms in the spherical harmonic model. Then, Werner et al. [117] sorted out the previous work, taking (4769) Castalia as an example, and introduced the modeling method of the polyhedral gravitational field in detail. Mirtich [118] also applied Gauss’s theorem and Green’s formula to replace the integral by summation and calculated the center of mass, the moment of inertia, the product of inertia, and other physical quantities of a homogeneous polyhedron.

The gravitational potential, gravitational force, and gravitational gradient tensors of the polyhedron at any point outside the homogeneous polyhedron can be expressed as [117,118]:

$$U(r)\frac{1}{2}G\sigma {\displaystyle \sum _{e\in E}{L}_e}(r_e·E_e·r_e)-\frac{1}{2}G\sigma {\displaystyle \sum _{f\in F}{\theta }_f}(r_f·F_f·r_f)$$

(7)

$$\nabla U(r)=G\sigma {\displaystyle \sum _{e\in E}{L}_e}(E_e·r_e)-G\sigma {\displaystyle \sum _{f\in F}{\theta }_f}(F_f·r_f)$$

(8)

$${\nabla }^{2}U(\mathrm{r})=G\sigma {\displaystyle \sum _{e\in E}{L}_e}{E}_e-G\sigma {\displaystyle \sum _{f\in F}{\theta }_f}{F}_f$$

(9)

In the formula, σ is the density of the homogeneous polyhedron P, E is the set of all edges on the surface f, F is the set of faces of the polyhedron P, $r_e$ represents the vector from $r$ to any point on edge e, $r_f$ represents the vector from $r$ to any point on the side f, L_e, E_e, F_f are the quantities related to the edge and the side, and θ_f is the solid angle formed by the points at the side f and r. Its specific calculation formula is

$${\theta }_f=2\arctan \frac{r_1·(r_2\times r_3)}{r_1{r}_2{r}_3+r_3{r}_1·r_2+r_1{r}_2·r_3+r_2{r}_1·r_3}$$

(10)

where $r_1$ $r_2$ $r_3$ are the vector radii from the point at $r$ to the three vertices on the side of the triangle. Denote

$$\Omega ={\displaystyle \sum _{f\in F}{\theta }_f}$$

(11)

when the point is inside the polyhedron P, Ω = 4π; when the point is outside the polyhedron P, Ω = 0. From this, the positional relationship between the point and the polyhedron can be determined.

The polyhedron model has no truncation error; the error only comes from the shape error and numerical calculation error between the model and the real celestial body. The calculation accuracy is high, and the polyhedron is not necessarily a convex polyhedron. Moreover, it can be well simulated near the surface and even inside the asteroid, and it can meet the requirements of global calculations. In the analysis of orbital dynamics, the polyhedron method also easily judges whether the particle is outside the asteroid. The main disadvantage of the polyhedron model is the large amount of calculation. Every time the gravity calculation needs to be performed on all edges and vertices, the calculation speed will be greatly reduced when the number of edges and vertices increases. (101955) Bennu’s geometric shape based on the polyhedron model can be found in Figure 2.

The above methods have their own advantages and disadvantages and need to be appropriately selected according to the characteristics of specific problems.

3.2. Research on Orbital Dynamics near Small Celestial Bodies

Based on the various dynamical models described in the previous section, the dynamic research on the vicinity of small irregular celestial bodies mainly includes the manifold structure and the local motion at the equilibrium point and its vicinity, large-scale periodic orbits and their bifurcations and resonances, quasi-periodic orbits, chaos, dynamic configuration, and the evolution of binary asteroids or multi-star systems. Dynamic equilibrium points, periodic orbits, and quasi-periodic orbits are important ways to study the phase space structure of complex dynamical systems.

The research related to equilibrium points started with the dynamic problem near small celestial bodies and has the most relevant research so far. Early research focused on the existence, quantity, and stability of equilibrium points near special geometries. Zhuravlev [120] first studied the stability of equilibrium points near the three-axis spheroid and calculated the stable and unstable regions.

Scheeres et al. [121,122] also used the spherical harmonic gravitational field model to calculate the position of the equilibrium point of (4769) Castalia and analyze its stability. Elipe et al. [123] found four equilibrium points in the gravitational field of a finite straight segment and analyzed their stability. Scheeres et al. [124] calculated the positions of the four equilibrium points of (25143) Itokawa. Mondelo et al. [125] calculated the positions of four equilibrium points of (4) Vesta and analyzed the stability. Liu et al. [126] analyzed the manifold structure near the equilibrium points in the gravitational field of a rotating homogeneous cube and the heteroclinic orbits between different equilibrium points. Yu et al. [127] calculated the coordinates of four equilibrium points in the gravitational field of (216) Kleopatra, linearized the eigenvalues of the matrix, and analyzed the stability based on this. Scheeres [128] calculated the equilibrium points of (1580) Betulia and 67P/Churyumov-Gerasimenko.

Jiang et al. [129] gave a linearized equation for motions near the equilibrium points, deduced a sufficient and necessary condition for the stability of the equilibrium point, and studied the characteristic root distribution, stability, and topological types of equilibrium points. According to the sub-manifold structure, the non-degenerate equilibrium points are divided into eight categories, which is a major advancement for scientists in correctly understanding the relevant characteristics of the equilibrium points of small celestial bodies. Wang et al. [130] used the polyhedron model to calculate the positions of the equilibrium points of 23 small celestial bodies and analyzed their stability. In particular, they found eight equilibrium points near the asteroid Bennu. This shows that the number and distribution of equilibrium points near small bodies are diverse and cannot be completely divided into two types determined by a simple geometric model.

Regarding the change in equilibrium points with the normalization parameters of the density and the rotational speed of small celestial bodies, Jiang et al. [131] found that equilibrium points always appear or annihilate in pairs, and the number of non-degenerate relative equilibrium points is odd. Wang et al. [132] took Bennu as an example to summarize the bifurcations of equilibrium points.

Since the 1990s, scientists have also conducted many studies on the periodic orbits in the gravitational field of small rotating irregular celestial bodies. The main focus is on the search for periodic orbits, orbit classification and stability analysis, and the dynamic bifurcation behavior due to changes in parameters such as the rotation speed of small celestial bodies and the energy integral of particle motion. Scheeres et al. [121,133] calculated the periodic orbits in the equatorial plane of (4) Vesta and (433) Eros using a three-axis ellipsoid model. Scheeres et al. [134] also used the second-order quadratic gravitational field model to calculate the frozen orbits and periodic orbits near Tutatis and analyzed the effects of the C₂₀ and C₃₀ terms on the frozen orbits. Antreasian et al. [135] and Scheeres et al. [136] successively used the second-order quadratic gravitational field model and the average method to analyze the motions near (433) Eros and found a family of retrograde periodic orbits, which were used for the Shoemaker mission. Shang et al. [137] investigated various periodic orbits near non-principal-axis rotation asteroids.

Scheeres et al. [138,139,140,141,142] also studied the influence of the C₂₀ and C₂₂ terms on the energy and angular momentum of the particle motion and numerically calculated the stable and unstable orbital regions in the parameter space. They further studied the orbital dynamics considering the non-spherical gravity of small celestial bodies, solar radiation pressure, and solar gravity and used the averaging method to find the frozen orbits near small celestial bodies. The property is closely related to the area-to-mass ratio of the spacecraft and the distance from the small celestial body to the Sun. Because the search for periodic orbits is very complicated, it is generally necessary to use symmetry for analysis and research, but the gravitational field of small irregular celestial bodies does not have this feature [128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143]. Yu et al. [144] proposed a global search method for 3D periodic orbit families in the vicinity of irregular small bodies by using the polyhedral model and the hierarchical grid method. The 29 periodic orbits are given by taking asteroid (216) Kleopatra as an example. The periodic orbits are classified into seven types according to the orbits of the four-dimensional symplectic manifold by calculating the eigenvalues of the monodromy matrix of periodic orbits [145]. Topological types, the bifurcation phenomenon, and the stability of periodic orbits with the continuation of energy are studied. Jiang et al. pointed out that periodic orbits move in a six-dimensional symplectic manifold and that its manifold structure is different from that of the four-dimensional case and re-classified the periodic orbits near small irregular celestial bodies into 13 topological types [132]. Applying this theory, Yu et al. [146] found periodic orbit families belonging to different topological types near (243) Ida and the bifurcation behavior in the continuation process in the periodical orbit search and continuation near (243) Ida. Jiang’s theory provides a powerful tool for follow-up research to better understand the type and stability of periodic orbits near irregular small celestial bodies from the topological structure. Non-equatorial equilibrium points near an asteroid with gravitational orbit-attitude coupling perturbation were analyzed in reference [147]. Li et al. [148] calculated the geophysical environments and periodic orbits near 2016 HO3 by using different shape models.

With the development of orbital dynamics and calculation methods for spacecraft near asteroids, people’s research on small celestial bodies has gradually expanded from single asteroids to binary asteroid systems and even systems comprised of multiple small bodies. Liang et al. [149] used the Poincare section and Jacobi constant to find the homoclinic and heteroclinic orbits connecting the equilibria near the contacting binary asteroids, which provided an important reference for designing low-energy transfer orbits between equilibria. Hou and Xin [150] constructed an explicit first-order solution to the rotations and the orbital motion in the planar two-body problem. Shi et al. [151] studied the equilibrium points and periodic orbits of the binary asteroid system (66391) 1999KW4 by the homotopy method based on the restricted full three-body problem. Wang and Fu [152] constructed a semi-analytical model for spacecraft near the primary of a binary asteroid system based on a perturbed two-body problem. For nonlinear dynamics of multiple body systems, the discrete element model is usually applied to simplify related calculations [153]. Jiang et al. [154] analyzed the dynamical configurations of the five triple-asteroid systems 45 Eugenia, 87 Sylvia, 93 Minerva, 216 Kleopatra, and 136617 1994CC and the six-body system 134340 Pluto. Figure 3 shows the dynamical configuration of binary asteroid system (66391) 1999 KW4. Valvano et al. [155] discussed the stability regions near the triple asteroid system 2001 SN 263 considering the effect of irregular shapes and the solar radiation pressure.

Due to the complex orbital shape near the small celestial body, the theoretically calculated periodic orbit may not be closed. Wang [156] expanded the definition of the resonant orbit from the orbital period satisfying the conventional relationship to the angular velocity of the particle moving around the small celestial body. In addition, some theoretical periodic orbits are affected by various perturbations, and the orbits will not be closed. Scheeres et al. [134] studied the quasi-periodic frozen orbit near (4179) Toutatis under the second-order quadratic gravitational field model. Chanut et al. used the polyhedron model to study the long-term motion of particles near the small celestial bodies (433) Eros and (216) Kleopatra [157,158].

Chaos phenomena are ubiquitous in the natural world and in engineering problems such as chemistry, mechanical systems, financial economics, and nanoscience [159,160,161,162,163,164]. The generation of chaotic phenomena in dynamics is usually closely related to bifurcation and resonance phenomena [165]. Studies on simple geometric models of gravitational fields have uncovered the chaotic behavior of particles. Elipe et al. [123] found the bifurcation caused by 1:1 resonance in the gravitational field of a finite straight segment and chaos due to parameter changes. Lindner et al. [166] discovered the chaotic phenomenon of particles moving around a line segment. Makarov et al. [167] found the chaotic phenomenon of the rotation of triaxial asteroids and minor planets. Since most of the small celestial bodies have strong irregular shapes, the resonance mechanism in the double asteroid or multi-asteroid system is more complicated, and the influence of the mutual coupling of orbit and attitude must be considered. Based on the sphere-ellipsoid model, Nadoushan et al. [168] found that aspheric factors and orbital eccentricity can significantly affect the size of the resonance region. After a certain critical value, the resonance regions intersect and lead to the appearance of chaos.

The YORP effect is a thermal-radiation torque acting on small asteroids and plays a crucial role in their physical and dynamical evolution. A detailed introduction can be found in Section 3.4. Cuk Nadoushan et al. [169] considered the solar gravitational perturbation, the orbital attitude dynamics, and the long-term evolution of the binary asteroid system under the YORP (Yarkovsky-O’Keefe-Radzievskii-Paddack) effect. They found the chaotic behavior of the attitude under the condition of small orbital eccentricity, even with the same rotation period. The widespread existence of chaotic phenomena makes it very important to distinguish between ordered and chaotic motions. Lyapunov Characteristic Exponents (LCE) [170] provide the distinction criterion from theoretical and numerical viewpoints and provide the chaotic intensity quantitative characterization. However, in practice, the numerical computation required to discover chaotic phenomena is time-consuming, especially for some chaotic motions that are very close to ordered motions. Froeschlé et al. [171,172] and Fouchard et al. [173] successively developed the Fast Lyapunov Indicator (FLI) and the Orthogonal Fast Lyapunov indicator (FLI) in response to the shortcomings of LCE. This provides an effective indicator for effectively distinguishing ordered and chaotic motions. Ni et al. [174] proposed to quantitatively analyze the indicators of quasi-periodic orbits from the perspective of frequency domain analysis and specifically analyzed the indicators of orbits in the gravitational field that neither escape nor collide with small celestial bodies. Among them, a complex orbit in the gravitational field of (6489) Golevka can be seen in Figure 4.

3.3. Research on Surface Motion Dynamics and the Capillary Phenomenon of Small Celestial Bodies

After exploring the special orbits near small bodies, people naturally pay attention to the selection of the soft-landing region and trajectory optimization, which are closely related to the transition dynamics on the surface of small celestial bodies. Generally speaking, the matter on the surface of irregular small celestial bodies may undergo transitions or even launch and escape behaviors. These motions are specifically classified into surface equilibria, motion confined on the surface, surface transition, and bouncing. For example, ESA’s Rosetta probe’s lander Philae clearly observed two bounces when it touched down on 67P/Churyumov–Gerasimenko, which means that Philae was observed to make three landings. It is noteworthy that the interval between the first landing and the final resting on the surface of the comet nucleus is 2 h. ESA did not fully consider the contact mechanics, collision, and bounce of the soft-landing process on the surface of irregular bodies. This led to the fact that the final soft-landing position of the lander was hundreds of meters away from the preset position.

A large number of previous studies have studied the physical and chemical properties of the surface particles of small celestial bodies, including the electrokinetic and rotational emission of dust particles in cometary cores [175] and the mineralogy and mineralogy of asteroid dust particles [176]. Moving particles and dust may be caused by a variety of reasons, including the YORP effect [177], the windmill effect on the surface of small celestial bodies [175], the collision and gravitational reconstruction of asteroid families and small moons [178], and the gravel disintegration of rubble-pile asteroids [179,180]. The disintegration of asteroid P/2013 R3 produced more than 10 distinct small bodies, numerous small grains, and a comet-like dust tail [180]. Most asteroids and comet nuclei have irregular shapes, and particles can move on the surface of small irregular bodies. To understand the dynamical behaviors of particles on the surfaces of irregular celestial bodies, we need to study the surface mechanical environment. In addition, the movement of the lander or the rover on the surface of the irregular small celestial body also needs to be studied. If the lander lands on the surface of an irregular small body, the collision and bounce of the lander on the irregular surface also exist [181].

Beginners choose simple shapes, including ellipsoids [182,183] and cubes [184], to help understand the surface motion. Guibout and Scheeres [182] discussed the existence and stability of surface equilibrium points for spheroids of revolution. Bellerose and Scheeres [183] used ellipsoids to simulate the shape and gravity of asteroids, studying hopping on flat surfaces. Belleros et al. [181] considered the motion and control of the surface exploration robot under the gravitational force of the uniformly rotating ellipsoid. Liu et al. [184] calculated the positions and eigenvalues of the surface equilibrium of a uniformly rotating cube. The non-degenerate equilibria in the gravitational field of general irregular celestial bodies can be divided into eight different types [129]. The motion confined on the surface of the irregular small celestial body is different from the motion in the gravitational field [185]. The former needs to consider the irregular gravitational force and the contact force, while the latter only considers the irregular gravitational force. Generally speaking, the transition or landing process of particles or landers on the surface of irregular celestial bodies includes orbital motion, collision, jumping motion, surface motion, and surface equilibria.

Considering a non-smooth surface with a constant coefficient of friction, the equilibria remain but are not as stable as the equilibria on a smooth surface [184]. Considering the precise gravity and irregular shape, Yu and Baoyin [186,187] numerically calculated the motion and particle migration of the rover on the asteroid surface and found that the most stable direction is the rotational pole direction, which can limit the rover’s movement after landing. However, the friction phenomenon on the surface of irregular small celestial bodies has a stick-slip effect [188], and the surface transition particles may be charged [189]. The orbital motion, collision and jumping motion, surface motion, and surface equilibrium of particles released over three different regions (flat surface, concave region, convex region) relative to asteroid 6489 Golevka were investigated in [185]. The results showed that when the particles were released over a flat surface and concave region, the surface equilibrium can be reached in a short time.

In the research of water ice material in small celestial bodies, previous studies in several pieces of literature found that the water on Earth comes from asteroids [190,191]. Kanno et al. [192] analyzed the wavelengths of the infrared spectrum to confirm the presence of water and ice on class D asteroids. Asteroids may release meteoroids, which fall to Earth as meteor showers or meteors, bringing material to Earth [193,194,195]. Treiman et al. [193] studied the Eucrite from the asteroid (4) Vesta and found quartz in the meteorite, arguing that quartz originates from the deposition of liquid water and that water may originate from (4) Vesta. Campins et al. [196] reported the presence of water and ice on the surface of (24) Themis and the widespread distribution of this water and ice. Comets also tend to have water on them. Sunshine et al. [197] detected the presence of water and ice on Comet 9P/Tempel, indicating that the surface deposits of the comet nucleus are loose aggregates. Taylor [198] reported in detail the presence of water in Eucrite meteorites originating from the asteroid (4) Vesta. Zolensky et al. [199] discussed temperatures of alteration, water:rock ratios, and the oxygen isotopic composition of water by analyzing the record of low-temperature alteration in asteroids. Trigo-Rodríguez et al. [200] presented the action of water in asteroids by studying carbonaceous chondrite meteorites.

The research on capillary action on asteroid surfaces is closely related to the equilibrium points and surface motion of particles. The height of the water in the capillary tube on the surface of the asteroid depends on the irregular shape and the gravitational potential of the asteroid. Different surfaces produce different heights of water in the capillary tube, so the friction coefficients of different surfaces are different, resulting in the stability of the surface balance sex being different. It was found that the gravitational field and spin velocity of asteroids have a significant effect on the height of the liquid in the capillaries on the asteroid surface [201]. This research can be applied in the following four aspects: ① This research helps to further study the distribution of water and ice on the surface of asteroids, and the different distributions of water and ice are related to the different distributions of liquid lengths in the capillary [196]. ② On the scale of millions of years, water can corrode the surfaces of asteroids and the structure inside, and the surface material and shape of a large number of asteroids have changed through surface erosion, especially for rubble-pile asteroids. In other words, long-term erosion could cause asteroids to break up and disintegrate. Aqueous alteration processes in asteroids were investigated by Rotelli et al. [202] to better understand the increasing complexity towards prebiotic chemistry. ③ The different heights of the liquid in the capillary can affect the electrokinetic and rotational jets of gas and dust particles on the surface of small celestial bodies [189]. Under the action of sunlight pressure, the jet may form a mini fountain on the surface of small celestial bodies. The change in the length distribution of the liquid in the capillary causes the height of the fountain and the radius of the fountain envelope to change [189]. ④ The probe can carry some liquid to the surface of the asteroid, so the study of the height of the liquid in the capillary may also be applied to future asteroid missions.

3.4. Dynamic Characteristics under Varying Parameters

The rotational speed, density, shape, internal structure, and other physical parameters of small celestial bodies vary widely, resulting in very different dynamic behaviors in their gravitational fields. Taking the rotational speed as an example, there is only one equilibrium point in the gravity field of 1998KY26 when it is rapidly rotating. For (52760) 1998 ML14, comet 1P/Halley, and 9P/Tempel 1, the external equilibrium points are farther away from the small celestial body. Asteroid (216) Kleopatra has seven relative equilibrium points, which is different from the five relative equilibrium points of most asteroids [131]. The YORP effect of small celestial bodies is closely related to thermosphysical parameters. It is a phenomenon that the rotation axis and rotation speed of small celestial bodies vary slowly due to the photon moment generated by sunlight [203]. Studies have shown that the YORP effect can slowly change the rotation speed of small celestial bodies and even cause the rotational disintegration of small celestial bodies. For example, the YORP effect can make the rotation rate of the small celestial body 54509 (2000PH5) accelerate by $(2.0\pm 0.2)\times {10}^{-4}$ (°)/d² [204], and the rotation rate of the small celestial body (1620) Geographos can be accelerated by $1.15\times {10}^{-8}$ rad/d² [205]. Numerical experiments show that the YORP effect can lead to the disintegration of rubble-like celestial bodies and the formation of small moons. The YORP effect may also indirectly affect the distribution and topological characteristics of the relative equilibrium points in the gravitational field. The shapes of asteroids may also be deformed as landslides and mass shedding occur. Similar to Comet Shoemaker-Levy 9, small rubble-like bodies can change their topography dramatically as they approach a planet, and some even disintegrate. Holsapple et al. [206,207] established the mechanical mechanism of tidal deformation caused by the influence of nearby larger celestial bodies. The variations in the shape, spin, and state during the slowly increasing angular momentum of rubble-pile, self-gravitating, homogeneous ellipsoidal bodies were investigated in reference [208,209]. Zhang et al. [210] found that three typical tidal response outcomes may appear on rubble piles, namely, deformation, scattering, and destruction. During the long-term evolution of small celestial bodies, their density, rotational speed, shape, and internal structure may change. The disintegration and gravitational aggregation of asteroids will also affect their internal structure, average density, and shape. These factors lead to changes in the parameters of small celestial bodies. If the parameters of the primary in binary asteroids with a large-scale ratio vary, it is bound to have an impact on the movement of the small moon. In addition, changes in parameters will also affect the movement of dust, particles, and gravel in the gravitational field of small celestial bodies.

Tanbakouei et al. [211] studied the mechanical properties of particles on the surface of asteroid 25143 Itokawa using the nanoindentation technique. They found that these particles of asteroid regolith can be more compacted than the minerals forming the particular LL chondrite associated with potentially hazardous asteroids. For the DART mission, the impact with the secondary of the Didymos system will cause a momentum transfer from the spacecraft to the binary asteroid. It is expected to change the orbit period of this system and force it to librate in the original orbit [212]. Furthermore, the primary may be reshaped due to landslides or internal deformation during this process. A detailed analysis can be seen in reference [213]. In 2024 October, the ESA Hera mission will be launched to obtain a detailed characterization of the physical properties of binary asteroids and of the crater caused by the DART mission [214]. For ringed asteroids, a change in the parameters may also lead to a variation in the parameters of the ring and even cause the formation or disappearance of the ring.

The influence of parameter changes on the dynamic behaviors in the gravitational field of small celestial bodies is very complex. Many scholars usually assume uniform mass distribution when modeling small celestial bodies. However, there may be mascons inside small celestial bodies, which make the mass distribution uneven. The density and distribution will inevitably have an impact on the gravitational field environment. Based on the polyhedron model, Chanut et al. [215] established a small celestial body model considering mascons. Aljbaae et al. [216] studied the gravitational field environment and the position of equilibrium points of the (21) Lutetia asteroid in the three cases: no mascons, three-layer mascons, or four-layer mascons. Chanut et al. [217] studied the orbital stability in the equatorial plane of the asteroid (101955) Bennu under the conditions of uniform mass and non-uniform mass distribution and found that, for Bennu, it is more appropriate to divide the mascon structure into 10 layers. Jiang [218] considered the position and topology transitions of the equilibrium points in the gravitational field of (2867) Steins with single-layer similarity, single-layer spherical, multi-layer similarity, and multi-layer spherical mascons. Jiang et al. [219] studied the bifurcation types corresponding to the collision and annihilation of equilibrium points during increasing the spin rate of (216) Kleopatra. It was found that the number of non-degenerate equilibrium points is changed from seven to five, then to three, and, finally, to 1. Moreover, they found a conserved quantity about the equilibrium points and deduced that the number of non-degenerate equilibrium points in the gravitational field of small celestial bodies can only be an odd number. Considering the shape effect on the environment of the gravitational field, the homotopy analysis method was used to generate continuous shape variation from small celestial bodies to simple symmetric geometric bodies (such as spheres, ellipsoids, and cubes), and the bifurcation phenomena of equilibrium points in the gravitational field are studied in [220].

4. Summary and Future Development

This paper introduces the general situation of small celestial bodies, summarizes the history and status of international missions, and sorts out the research progress of the orbital mechanics of spacecraft, such as modeling the gravitational field, orbital dynamics, surface motion dynamics, and dynamic properties under varying parameters. First, building an accurate gravitational model of small bodies lays the foundation for research on orbital dynamics. Taking into account the resource constraints of the onboard computer, how to balance computational efficiency and computational accuracy is still a key issue in modeling the gravitational field of asteroids. Second, considering the needs and constraints of real missions, the design of low-energy transfer orbits and corresponding control schemes deserves more attention. It is an important way to make full use of the orbital dynamics near asteroid systems, such as invariant manifolds and heteroclinic orbits. Moreover, studying the evolution of multiple asteroid systems is of great significance for understanding the origin of the solar system. Finally, the impact monitoring of near-Earth objects and estimating the impact probability with the Earth is the first step towards planetary defense. Furthermore, investigations on the means of deflecting a potentially hazardous object and the evaluation of the defense effectiveness also need the jointed efforts of international communities.

China achieved the first fly-by exploration of small celestial bodies in 2012 and plans to implement the mission of orbiting, landing on, and obtaining sampling returns from small celestial bodies in 2025. At present, the exploration mission of small bodies has entered the engineering development stage. The core technology needed to realize the orbiting of and landing on small celestial bodies lies in the orbital dynamics and control of spacecraft. The strong irregular shape of small celestial bodies causes their gravitational fields to be very different. This brings great challenges to the design and control of the orbit. Therefore, it is necessary to conduct a more comprehensive and in-depth study on related orbital mechanics.