Objectives The irregular geometry of small celestial bodies significantly affects gravity potential estimation. However, the commonly used polyhedral model (PM) suffers from high computational complexity, which limits its application to theoretical research and engineering practice.

Methods To address this issue, a gravity potential estimation method based on Monte Carlo (MC) integration is proposed for irregularly shaped small celestial bodies. The proposed method aims to reduce computational complexity while maintaining estimation accuracy.

Results The proposed MC method shows good agreement with the PM in gravity potential estimation, demonstrating its effectiveness. In contrast, although the spherical harmonic method is computationally efficient, it is less suitable for irregularly shaped small celestial bodies because data transformation and the limited capability of spherical harmonic expansion to represent irregular geometries inevitably introduce estimation errors. The gravity potential estimation results for Phobos and Comet 67P further demonstrate the reliability of the proposed method. Compared with the PM, the MC method substantially reduces computational complexity while maintaining comparable estimation accuracy. Application of the proposed method to Psyche reveals a distinct zonal distribution of gravity anomalies between the equatorial and polar regions. The results indicate that oblateness-induced perturbations may significantly affect the orbit determination of orbiting spacecraft. In addition, conducting landing and sampling missions in the polar regions of similar small celestial bodies may help reduce fuel consumption.

Conclusions The proposed MC method provides reliable gravity potential estimation with low computational complexity and is well suited for gravity potential and gravity anomaly estimation of irregularly shaped small celestial bodies. It provides technical support for related theoretical research and future exploration missions.