Accelerating gravity calculations using nested surface triangulation grids

For decades, one of the deceptively simple problems in celestial mechanics has been this: given an irregularly shaped body of uniform density, how do you compute its gravitational potential quickly and accurately? For a perfect sphere or an ideal ellipsoid, the answer is a textbook exercise. For a potato-shaped asteroid like 433 Eros, or a dumbbell-shaped cometary nucleus, or a torus of orbiting debris, there is no closed-form solution, and scientists have had to rely on clever numerical workarounds. The most elegant of these, the polyhedral method, treats the body’s surface as a mesh of flat facets and reduces the daunting volume integral of Newtonian gravity to a sum of elementary analytical contributions from each facet. A new study by Jean-Marc Huré of the University of Bordeaux’s Laboratory of Astrophysics, published in Celestial Mechanics and Dynamical Astronomy, now revisits this classic technique and reveals how much faster and more accurate it can become when combined with nested grids and a classical acceleration trick borrowed from numerical analysis.

The polyhedral method traces its lineage back to the mid-1970s, with landmark work by Waldvogel on the gravity of a homogeneous cube and by Paul on three-dimensional geophysical interpretation, and it was later brought into planetary science by Werner in a widely cited 1994 paper whose title—””Don’t Cut Corners””—remains a warning to anyone careless with mesh geometry. The idea is beautiful in its economy: the gravitational potential of a uniform body can be rewritten, via the divergence theorem, as a surface integral over its boundary. When that boundary is tiled with triangular facets, the contribution of each triangle to the potential at any point can be expressed in closed analytical form, involving geometrically meaningful quantities such as the solid angle subtended by the triangle and the perpendicular distance from the point to the triangle’s plane. Because each contribution is exact for the flat facet, the only source of error is the discrepancy between the faceted surface and the true smooth surface of the body. Huré’s new work asks the question that follows naturally: if we refine that mesh, how fast does the numerical potential converge to the true one, and can we do better than brute-force refinement?

To answer it, Huré constructed a battery of six test shapes of genuine astrophysical relevance: a sphere, a spheroid, a triaxial ellipsoid, a lemon-shaped figure, a dumbbell, and a torus. Each of these has known or independently computable reference values for the gravitational potential, which makes them ideal benchmarks. The surfaces were projected onto nested grids—that is, families of triangular meshes in which each finer grid contains the coarser one, with the number of nodes doubling along each direction at every level of refinement. This nesting is not merely a convenience; it is the structural prerequisite for the convergence-acceleration technique that forms the heart of the paper. The analysis was carried out both inside the bodies, at their surfaces, and in the exterior field, covering the full range of situations encountered in practice.

The central numerical finding is clean and satisfying: the polyhedral method, applied to triangular surface meshes, converges to the reference potential at second order in the grid step size. In other words, halving the typical facet size reduces the error in the computed potential by roughly a factor of four, uniformly inside, on, and outside the body. This second-order behavior is the theoretical signature that makes the next step possible, because it tells you exactly how the leading error term depends on the mesh spacing. Knowing the order of convergence is like knowing the slope of a road before you drive it—it lets you extrapolate with confidence rather than guesswork.

That next step is the Repeated Richardson Extrapolation, a venerable technique that goes back to Lewis Fry Richardson in the early twentieth century and is a staple of every serious numerical methods course. The idea is simple to state. Suppose you have computed some quantity at two different resolutions, one twice as fine as the other, and you know the order of the leading error term. You can then combine the two results to eliminate that leading error entirely, obtaining an estimate as accurate as one from a far finer grid. Richardson’s insight was that this procedure can be iterated: by computing the quantity at a sequence of successively finer nested grids and forming a triangular table of corrected values, each pass through the table cancels one more term in the error expansion. When the underlying convergence is second order, the effective error after each extrapolation level drops dramatically, scaling with ever-higher powers of the mesh spacing. In his paper, Huré adapts the standard formulas to the specific structure of his spherical grids, where the triangular tiling near the poles requires a modified set of combinatorial factors in the extrapolation table.

The practical payoff is striking, but only under the right conditions. For smooth bodies—the sphere, the spheroid, the ellipsoid, the lemon, the dumbbell, and the torus of his test suite—the extrapolated values of the potential converge far faster than raw grid values at the same resolution, meaning that a desired accuracy can be reached with a much coarser mesh and therefore much less computing time. This matters enormously in two areas of theoretical astrophysics. The first is the theory of figures: the classical problem, going back to Newton, Maclaurin, Jacobi, and Chandrasekhar, of determining the equilibrium shapes of rotating self-gravitating fluid bodies, whether single-layer or multi-layer. Solving for the figure of a rotating fluid requires evaluating the gravitational potential hundreds or thousands of times as the shape is iterated toward equilibrium, so any speedup in each evaluation compounds across the entire calculation. The second area is dynamical studies: tracking test particles in the field of an irregular body, or computing the mutual gravitational interaction of two polyhedral bodies such as a binary asteroid system. Again, the potential must be evaluated a vast number of times, at every time step and for every particle, and the computational cost of the polyhedral method scales steeply with the number of facets—roughly with the fourth power of the number of nodes in Huré’s spherical tiling. Cutting that cost by extrapolating from coarse grids rather than brute-forcing fine ones is not a luxury; it can be the difference between a feasible simulation and an impossible one.

For real objects, however, the story changes in an instructive way. Asteroids, the most common targets of this machinery, are not smooth. Their surfaces are cratered, boulder-strewn, and rough at all observable scales, and the shape models derived from spacecraft laser altimetry and stereophotoclinometry—such as the celebrated Gaskell Eros shape model hosted on NASA’s Planetary Data System—encode that roughness explicitly. When the true surface is jagged at small scales, the smoothness assumption underlying Richardson extrapolation weakens, and Huré finds that the gain from the acceleration technique becomes modest. There is, after all, little point in extrapolating away the discretization error of a mesh if the mesh itself is faithfully representing bumps that a coarse grid would smooth over.

But this observation leads to the second major contribution of the paper: a criterion for when further mesh refinement stops buying you anything at all. If the shape model itself carries uncertainties—altimetric errors, registration errors, reconstruction noise—then at some point the error in the potential due to those shape uncertainties dominates over the error introduced by unresolved sub-grid cavities and bumps. Refining the mesh beyond that point is wasted effort, because you are computing with exquisite precision the gravity of a shape you do not actually know. For bodies that are close to spherical, Huré derives a remarkably compact criterion: the number of triangular facets T should satisfy T greater than or equal to 64D divided by 3λ, where D is the diameter of the body and λ is the typical shape error. Below this threshold, sub-grid terrain matters; above it, the shape uncertainty has taken over and extra facets are simply refining noise.

Huré illustrates the criterion with the archetype of a well-mapped asteroid: 433 Eros, the target of NASA’s NEAR Shoemaker mission in 2000 and 2001 and one of the most thoroughly characterized small bodies in the solar system. Eros is decidedly not spherical—it is a contact-binary-like elongated body some 34 kilometers long—but the analysis of its shape model, with its documented altimetric uncertainties, shows how the interplay between discretization error and shape error plays out for a real dataset. The result gives mission planners and dynamicists a principled answer to the practical question of how fine a shape model to use when computing orbits, surface slopes, or particle dynamics around such a body: fine enough to beat the shape uncertainty, but no finer.

Beyond the headline results, the paper contains a set of technical results that practitioners of the polyhedral method will find immediately useful. Huré shows that once the individual triangle contributions to the potential are computed and stored, the gravitational acceleration comes essentially for free, because each contribution’s gradient is aligned with the local facet normal and can be recovered from the potential contribution itself through an exact analytical relation—echoing the classical result of Holstein and Ketteridge from the geophysics literature. This avoids the numerical differentiation of the potential, which is always a source of additional error. The paper also works out elegant limiting cases, such as the potential at a vertex of a tetrahedron, which depends on a single triangular facet thanks to the geometry of the projection.

Taken together, the study is a reminder that classical numerical methods are not finished artifacts to be used passively but living tools that reward scrutiny. By verifying the second-order convergence of the triangular polyhedral method, grafting onto it the full machinery of repeated Richardson extrapolation, and delineating the point of diminishing returns imposed by real shape uncertainty, Huré has given the celestial mechanics community both a faster algorithm and a wiser one—a method that knows not only how to compute gravity accurately, but also when to stop.