Gravity forward modeling as a basic tool has been wildly used to topography correction and 3D density inversion. The source region is usually discretized into tesseroid (i.e., spherical prism) to consider the influence of the curvature of planets in global or large-scale problems. Traditional gravity forward modeling methods in spherical coordinates, including the Taylor expansion and Gaussian-Legendre quadrature, are all based on spatial domain, which are mostly low computational efficiency. This study proposes a high-efficient forward modeling method of gravitational fields in spherical harmonic domain, in which the gravity anomalies and gradient tensors can be expressed as spherical harmonic synthesis form of spherical harmonic coefficients of 3D density distribution. A homogeneous spherical shell model has been used to test its effectiveness compared with traditional spatial domain methods. It demonstrates that the computational efficiency of the proposed spherical harmonic domain method has improved by four orders of magnitude and with a similar level of computational accuracy compared with the optimized 3D GLQ method. Finally, the proposed forward method is applied to topography correction of the Moon. The results show that the gravity response of the topography by our method is close to that of the optimized 3D GLQ method, and is also consistent with previous results.