<p>We derive spherical harmonic coefficients for the external gravitational potential generated by a tesseroid with radially variable density. The Newtonian volume integrals required to compute the coefficients are decomposed into 2D surface integrals and 1D radial integrals. Analytical expressions are derived for the radial integrals corresponding to polynomial, exponential, and parabolic density functions, while an adaptive Gauss–Legendre quadrature is used to numerically evaluate the radial integrals for general smooth density distributions. Simplified expressions for the surface integrals are further obtained of special geometries, including spherical zonal bands, spherical sectors, and spherical shells. The derived formulas extend the tesseroid-based spectral technique for gravimetric forward modeling beyond the common assumption of constant density within each tesseroid. Benchmark tests indicate that our approach is reliable and offers a favorable trade-off of accuracy and efficiency. No pronounced polar-singularity and polar-region problems are observed, whereas a near-zone problem appears as evaluation points approach the source. This behavior is primarily attributed to truncation errors and can be mitigated by increasing the truncated harmonic degree, leading to good agreement with a spatial-domain double exponential quadrature method. Finally, we compute the gravitational field of global marine sediments with 3D variable density and provide the resulting potential coefficients model to support related geodetic and geophysical applications.</p>