<p>Accurate modeling of small-body gravity fields is critical for proximity operations, yet traditional exterior approaches such as spherical harmonics, polyhedral models, and mascon representations often face limitations in accuracy, convergence, or computational efficiency. This work introduces an interior cylindrical Bessel gravity field expansion, derived from Laplace’s equation in cylindrical coordinates, to provide a localized and computationally efficient representation of the gravitational potential. The expansion employs Fourier–Bessel functions, ensuring convergence within a predefined cylindrical region. Expansion coefficients are obtained via least-squares fitting to synthetic datasets, and the formulation is validated against a state-of-the-art polyhedral gravity model. Results demonstrate high-fidelity reconstruction, achieving sub-percent errors in potential and acceleration within the cylindrical region, including near the small body’s surface. Further validation through trajectory propagation and uncertainty prediction via linear covariance confirms the model’s suitability for spacecraft orbit determination during small-body operations. Compared with traditional methods, the cylindrical harmonic expansion achieves superior local convergence and accuracy with fewer coefficients. This framework enables precise gravity modeling in targeted regions such as landing or touch-and-go trajectories and offers promising extensions for improving orbit determination accuracy, analyzing particle dynamics, and conducting gravity popper experiments to support gravity science investigations and interior density inference.</p>