Highlights • A new spherical harmonic model of coma gas density is proposed. • An improved orbital model about a comet with coma gas, SRP, and non-spherical gravity is established. • A global grid method is used to search stable orbits around a comet. • The heliotropism and spatial distribution of stable orbits are analyzed.

Abstract As a frontier field of deep-space explorations, the study of comets is of great value to planetary scientific research, resource exploration, and planetary defense. To date, only the ESA's Rosetta spacecraft successfully orbited the comet 67P. The orbital dynamical environment around comets is more complex than that of other small bodies, making the orbital stability around comets an important issue. The parameters of Rosetta's target, comet 67P/Churyumov-Gerasimenko, are used in this study of the orbital dynamics around a comet. The most unique perturbation around a comet is the aerodynamic force caused by the coma gas. However, the current numerical models are too complicated for analytical studies, and the current spherical harmonic model is too simple to reflect the actual distribution of the coma gas field accurately. To overcome the disadvantages of current models, a method is proposed to reduce the numerical model to a simplified analytical model. Then, the simplified model is fitted with a new spherical harmonic model, which is suitable for both numerical simulations and analytical studies. By considering perturbative effects of coma gas, solar radiation pressure (SRP), and non-spherical gravity of the comet nucleus, the orbital dynamical model is then established. Next, with a global search method, stable orbits around 67P/Churyumov-Gerasimenko are searched, and their characteristics are analyzed. It is found that the perturbative effect of the coma gas reduces the stability of the orbits. As the distance between the comet and the Sun increases, the stable orbital region will gradually shift from low inclinations to higher ones, and all of the stable orbits are heliotropic. According to the results, some candidate exploration orbits around the comet are finally identified.