# Small bodies global gravity inversion via the level-set method

A. Caldiero 111email<EMAIL_ADDRESS>Earth and Life Institute,
Université Catholique de Louvain, Place Louis Pasteur 3, 1348 Louvain-La-
Neuve, Belgium Reference Systems and Planetology, Royal Observatory of
Belgium, 3 Avenue Circulaire, 1180 Uccle, Belgium  S. Le Maistre  Earth and
Life Institute, Université Catholique de Louvain, Place Louis Pasteur 3, 1348
Louvain-La-Neuve, Belgium Reference Systems and Planetology, Royal Observatory
of Belgium, 3 Avenue Circulaire, 1180 Uccle, Belgium

###### Abstract

We propose an approach to infer large-scale heterogeneities within a small
celestial body from measurements of its gravitational potential, provided for
instance by spacecraft radio-tracking. The non-uniqueness of the gravity
inversion is here mitigated by limiting the solutions to piecewise-constant
density distributions, thus composed of multiple regions of uniform density
(mass anomalies) dispersed in a background medium. The boundary of each
anomaly is defined implicitly as the 0-level surface of a scalar field (called
the level-set function), so that by modifying this field the shape and
location of the anomaly are varied. The gravitational potential associated
with a density distribution is here computed via a line-integral polyhedron
method, yielding the coefficients of its spherical harmonics expansion. The
density distribution is then adjusted via an iterative least-squares approach
with Tikhonov regularization, estimating at every iteration corrections to the
level-set function, the density contrast of each anomaly, and the background
density, in order to minimize the residuals between the predicted gravity
coefficients and those measured. Given the non-convexity of the problem and
the lack of prior knowledge assumed (save for the shape of the body), the
estimation process is repeated for several random initial distributions, and
the resulting solutions are clustered based on global properties independent
of the input measurements. This provides families of candidate interior models
in agreement with the data, and the spread of the local density values across
each family is used to assess the uncertainties associated with the estimated
distributions. We present multiple synthetic tests with increasingly more
realistic settings (in terms of gravity resolution and precision, and of
shape, size and distribution of the internal heterogeneities), showing that
our method is generally able to retrieve a ground-truth mass distribution even
with noisy data. For further validation, we present an application of the
method to real data, namely the Bennu gravity coefficients measured by the
OSIRIS-REx team.

###### Contents

  1. 1 Introduction
  2. 2 Methodology
    1. 2.1 From interior density to gravity field
    2. 2.2 Density function estimation
      1. 2.2.1 Level-set method
      2. 2.2.2 Iterative least-squares solution
  3. 3 Synthetic retrievals
    1. 3.1 Inversion quality metrics
      1. 3.1.1 Reduced $\chi^{2}$
      2. 3.1.2 Density correlation
    2. 3.2 Target models
    3. 3.3 Robustness
      1. 3.3.1 Gravity resolution
      2. 3.3.2 Noise
      3. 3.3.3 Shape and size of the body
      4. 3.3.4 Shape error
      5. 3.3.5 Grid discretization
      6. 3.3.6 Size, depth, magnitude of the anomaly
      7. 3.3.7 Initial model
    4. 3.4 Uncertainty estimation
    5. 3.5 Realistic cases
      1. 3.5.1 Continuous target distribution
      2. 3.5.2 Rubble-pile target
  4. 4 Preliminary analysis of Bennu real data
  5. 5 Discussion
  6. 6 Conclusions
  7. 7 Acknowledgements
  8. 8 Data Availability

##  1 Introduction

Any information on the current interior structure of a small body can provide
insights into its formation, as well as its collisional and dynamical
evolution, all strictly related to the history of the solar system itself
(Walsh, 2018). Although ground-based methods can allow to constrain the
internal mass distribution of asteroids, either in terms of their
macroporosities (Carry, 2012), or from their spin evolution (Lowry et al.,
2014), and close-encounters with our planet could prove useful to test
candidate mass distributions against their tidal response (Dinsmore and de
Wit, 2023), _in-situ_ observations remain the one reliable way to determine
properties of an asteroid beyond the bulk. Surface geological structures,
seismic analysis, tidal dissipation, and radar sounding (Kofman et al., 2020)
can all help infer the distribution of mass within a small body (Scheeres et
al., 2015), as can the measurement of its gravitational potential. The
external gravitational potential is a direct expression of its internal mass
distribution, and can be characterized from its perturbation on the motion of
spacecraft in its vicinity. The precise reconstruction of a spacecraft orbit
from radio-tracking to retrieve such signals (radio science) requires no
additional payload other than its telecommunication system, making gravity
ubiquitous among the outputs of a space mission to any celestial body. We aim
therefore to propose an approach to retrieve candidate mass distributions
inside small bodies from their gravity field, as measured in a radio science
campaign.

Robotic exploration of the Solar System has allowed for several encounters
with small bodies, be it in opportunistic flybys or in the frame of dedicated
missions. The NEAR Shoemaker mission to asteroid Eros (Miller et al., 2002)
opened an era of small bodies exploration which is far from reaching its peak,
with highlights such as the JAXA sample-return missions Hayabusa (Fujiwara et
al., 2006) to asteroid Itokawa and Hayabusa-2 (Watanabe et al., 2019) to
asteroid Ryugu, the NASA missions Dawn to Vesta and Ceres (Konopliv et al.,
2018) and OSIRIS-REx to Bennu (Scheeres et al., 2020), and the ESA Rosetta
spacecraft rendez-vous with comet 67P/Churymov-Gerasimenko (Pätzold et al.,
2016). While, due to the fast decay of the gravity signal with distance, most
of these encounters have only allowed to determine the central gravity of the
target, and thus its total mass, a handful of missions have successfully
measured the extended gravitational potential of their target (see Figure 3),
and the growing number of dedicated small bodies missions planned for the
coming years could increase this count. As the OSIRIS-REx spacecraft, fresh
from returning to Earth samples from asteroid Bennu, prepares its approach to
Apophis, various new encounters are expected within the next decade, with NASA
probes Lucy (Levison et al., 2021) and Psyche (Zuber et al., 2022) already on
the way to their targets (the Jupyter Trojans and the metallic asteroid
Psyche, respectively). The ESA Hera mission (Michel et al., 2022) is expected
to reach the Didymos binary system in 2027, already visited by the NASA DART
spacecraft in its planetary-defense demonstration of the deflection of the
secondary body, Dimorphos. Juventas, one of the two CubeSats on-board the Hera
spacecraft (the other being Milani), is planned to land on Dimorphos and
carries multiple scientific instruments, including a gravimeter (Ritter et
al., 2022) and a ground-penetrating radar (Herique et al., 2022), that would
provide local measurements complementing the global gravity field in the
constraining of the interior structure of the asteroid. Finally, the approach
we present here is equally suited to irregular bodies like the Martian moons
Phobos and Deimos, which will be targeted by the JAXA MMX mission (Matsumoto
et al., 2021).

As is well-known, returning from the measured gravity potential to the
internal mass distribution which generates it is an ill-posed problem:
attempts to model the interior of a body from gravity measurements have to
deal with the non-uniqueness of the solution and its instability (Blakely,
1995; Michel and Fokas, 2008; Chao, 2005). Unambiguous reconstruction of the
interior therefore relies heavily on the injection of _a priori_ knowledge.
For the Earth, where constraints from other observations (such as seismic or
magnetic) are widely available, this has been performed extensively both at
global and local scales (Zheglova et al., 2017), as reflected by the wide
range of commercial and open-source tools available (e.g. Rücker et al.,
2017). The interior of the planet is generally approximated by prismatic
elements for which the density is estimated via least-squares inversion with
various degrees of regularization (Li and Oldenburg, 1998), although spectral
or Monte-Carlo methods are equally popular (Michel and Fokas, 2008).

For small bodies, however, the initial knowledge about the internal structure
may be little to none. Therefore, except in special cases where additional
observables and theoretical considerations may limit the possible models to a
small set of families (Le Maistre et al., 2019), exploring the full parameter
space with forward approaches might in general prove infeasible. Moreover,
inverse approaches might fail to provide a full picture of the sets of
possible density solutions that fit the data within the noise. The problem of
non-uniqueness is addressed in Tricarico (2013), where forward methods are
employed to explore the full space of exact solutions, as found by determining
the null space of the matrix relating the measurements to basis functions of
the density expansion. The basis functions are generally orthogonal
polynomials, which makes such an approach less suited to model interiors with
sharp contrasts between density regions (although a mixed model where
additionally the shape of an anomaly is adjusted iteratively as in Silva and
Barbosa (2006) is already discussed in Tricarico (2013)). We focus here on the
complementary case of discrete regions of constant density, which might well
approximate bodies that have undergone shattering or reaccretion. In
particular, we think that our piece-wise constant approach should be well
suited to small bodies because their low gravity most likely left their
building blocks intact during the accretion process (no material mixing,
important density jump, possible large cavities inside). Other proposed
approaches to the gravity inversion of small bodies include multiresolution
(Sorsa et al., 2020) or Markov-chain Monte-Carlo (Izquierdo et al., 2023)
bayesian inversion, matrix inversion (Ledbetter et al., 2021), and machine
learning (Izzo and Gómez, 2022) approaches. In Park et al. (2010), the body
was discretized as cubes and a density was estimated for each element, but the
solutions were highly affected by the decay of the sensitivity for deeper
elements. Takahashi and Scheeres (2014a) combined a linear inversion of the
density over a limited number of regions of the body with a forward approach
to determine the shape of these regions. Comparison of these methods with our
proposed approach will be addressed in Section 5. In recent years there have
been multiple real-case applications of gravity inversions for asteroids,
starting from the Dawn mission analyses of Vesta (Ermakov et al., 2014) and
Ceres (Ermakov et al., 2017). The gravity data of OSIRIS-REx for Bennu has
been extensively analyzed by Scheeres et al. (2020), where the polynomial
inversion of Tricarico (2013) was combined with a forward approach itself
driven by analytical modelling and morphological constraints to derive a
possible interior structure of Bennu composed of a light central core and an
equatorial ring of lower density than the mantle. From the same dataset,
Tricarico et al. (2021) were able to propose interior models under the
assumption of Bennu being a rubble-pile. As demonstrated by these efforts and
due to the non-uniqueness of the problem, inverse methods should still be
combined with forward approaches in the exploration of the solution space and
possibly in providing realistic uncertainties associated with each solution.

The proposed approach, from hereon referred to as GILA (Gravity Inversion via
Level-sets for Asteroids), relies on the use of level-set functions to
estimate the shape and location of density anomalies within the body, a
technique already well-established in many fields (Hedges et al., 2017),
including Earth local-gravity inversion (Giraud et al., 2021, and references
therein), but to the best of our knowledge not yet applied to gravity
inversion of other bodies. We strive to keep our approach free from initial
assumptions, although a synergy with other methods of interior inference is
vital to resolve the inherent degeneracy of the problem. Therefore, while in
the current analysis we limit ourselves to the sole use of gravity data and
assume no external information about the interior distribution, further
efforts will be required to inject into our approach constraints from other
observables and different gravity inversions paradigms.

An outline of GILA is presented in Section 2, with Section 3 displaying its
application in the recovery of a known ground-truth model from simulated
observables. Finally, Section 4 shows how the method fares when applied to
real data from the OSIRIS-REx mission.

##  2 Methodology

GILA figures as an extension to small bodies of methods well established in
Earth local gravity inversion. We therefore give a high-level overview of the
aspects shared with these methods, while focusing on the characteristics
implemented to make it suitable to the global gravity inversion of small
bodies.

###  2.1 From interior density to gravity field

The classical expression of the gravitational potential of a generic body
makes use of a spherical harmonics expansion (e.g. Heiskanen and Moritz,
1967):

$U(r,\theta,\phi)=\frac{G\mathcal{M}}{r}\sum_{l=0}^{\infty}\sum_{m=0}^{l}\left(\frac{r_{0}}{r}\right)^{l}P_{lm}(\cos\theta)\cdot(C_{lm}\cos
m\phi+S_{lm}\sin m\phi)$ (1)

with $G$ the universal gravitational constant, $\mathcal{M}$ the mass of the
body, $r_{0}$ the reference radius, $\left(r,\theta,\phi\right)$ the spherical
coordinates (radius, colatitude, longitude), $P_{lm}$ the (associated) fully-
normalized Legendre polynomials of degree $l$ and order $m$, and $C_{lm}$ and
$S_{lm}$ the dimensionless fully-normalized Stokes coefficients. We will use
the symbol $\mathcal{O}_{lm}$ to indicate any of $C_{lm}$ or $S_{lm}$ when an
expression holds for both sets of coefficients. The $P_{lm}$ functions and
$\mathcal{O}_{lm}$ coefficients are obtained from their unnormalized
counterparts by respectively multiplying and dividing by the normalization
factor $\sqrt{(2-\delta_{m,0})(2l+1)(l-m)!/(l+m)!}$, $\delta_{m,0}$ being the
Kronecker delta. In practical applications, the expansion is generally
truncated at a given degree $l_{\max}$.

The Stokes coefficients are typically part of the set of solve-for parameters
in an orbit determination campaign, and thus represent the way the
gravitational potential is estimated in radio science. As is well known, this
expression is well suited for bodies of nearly-spherical shape, but less so
for irregular bodies, where the series might display convergence issues when
computed within their circumscribing (Brillouin) sphere. For small bodies,
alternative representations of the gravity, such as those based on ellipsoidal
(Park et al., 2014) or interior (Takahashi and Scheeres, 2014b) harmonics,
Fast Fourier Transform (Perez-Molina et al., 2022), or neural networks (Izzo
and Gómez, 2022) might prove more accurate. Nonetheless, spherical harmonics
are still usually provided as deliverables of radio science campaigns, for
reasons ranging from legacy in orbit determination software to their relative
simplicity and direct connection to specific perturbations in the spacecraft
orbit (Scheeres et al., 2016), and because in most orbit or fly-by
configurations the spacecraft does not enter the Brillouin sphere. It is for
these reasons that we choose to employ the Stokes coefficients as observables
in our study, meaning we assume they are provided as an external input, be it
from radio science or from complementary methods.

In turn, the normalized Stokes coefficients are related to the interior
density of the body via the integral (e.g. Jekeli, 2007):

$\begin{bmatrix}C_{lm}\\\
S_{lm}\end{bmatrix}=\frac{1}{(2l+1)\mathcal{M}}\int_{\mathcal{V}_{\mathcal{B}}}\rho{\left(\frac{r}{r_{0}}\right)}^{l}\cdot
P_{lm}(\cos\theta)\begin{bmatrix}\cos(m\phi)\\\ \sin(m\phi)\end{bmatrix}dV$
(2)

where $\mathcal{V_{\mathcal{B}}}$ indicates the total volume of the body.
Thus, given any expression for the density $\rho$ within the body (e.g. a
series of orthogonal functions), the Stokes coefficients can be computed via
Eq. 2. Similarly to Le Maistre et al. (2019), we discretize Eq. 2 by
considering the body as a collection of elements: the volume of the body
$\mathcal{V}_{\mathcal{B}}$ is assumed to be composed of a surface layer and
an interior layer. The interior layer is here generally represented by a
hexahedral grid (see Figure 1), but can be any kind of 3D mesh. The surface
layer consists of the region between the interior layer and the surface mesh.
We choose this structure to allow for a generic mesh discretization of the
interior which is independent of the surface mesh, itself assumed to be an
external input and usually composed of triangular elements. The Stokes
coefficients of the body are then computed as a weighted sum of the
contributions of each element (cells of the interior grid, plus the surface
layer), the weight being the element mass:

$\mathcal{O}_{lm}=\frac{1}{\mathcal{M}}\left[\rho_{s}\mathcal{V}_{s}\mathcal{O}_{lm}^{s}+\sum_{i=1}^{N}\rho(\bm{x}_{i})\mathcal{V}_{i}\mathcal{O}_{lm}^{i}\right]$
(3)

where $N$ is the total number of grid cells, $\mathcal{V}_{i}$ the volume of
the $i^{th}$ cell, and $\bm{x}_{i}=(x_{i},y_{i},z_{i})$ the coordinates of its
center, while the surface layer is assigned a constant density $\rho_{s}$ and
has a volume $\mathcal{V}_{s}=\mathcal{V}_{B}-\sum_{i=1}^{N}\mathcal{V}_{i}$.
The Stokes coefficients of each element, here denoted as
$\mathcal{O}_{lm}^{s}$ and $\mathcal{O}_{lm}^{i}$, are computed via the line-
integral approach of Jamet and Tsoulis (2020). Compared to other well-
established methods for the computation of the Stokes coefficients of a
generic polyhedron (Werner, 1997; Tsoulis et al., 2009), this method is non-
recursive, allowing for vectorization in the gravity degree-and-order
dimensions. While not analytical, it has sufficient accuracy at the low
degrees we consider in our study even with a coarse discretization of the line
integrals (we generally use $2^{3}$ line integration points for each edge of
the exterior mesh and the internal grid). Moreover, it is not restricted to
polyhedra with a vertex at the origin of the potential expansions, allowing
for a more immediate modelling of heterogeneities within the body, compared
for example to the approach of Takahashi and Scheeres (2014a). Overall, the
computational overhead of this polyhedron model implementation compared to the
simpler mascons approach (e.g. Tardivel, 2016) is not impractical, especially
when high accuracies, and therefore more mascons, are required.

Eq. 3 provides a linear relation between the gravity coefficients and the
density of each grid cell, itself computed as the value of the density
function at the location of the cell center ($\bm{x}_{i}$), which corresponds
to a node in the discretization grid. The density distribution within the body
is here assumed to be piecewise-constant, meaning that the elements are
grouped to form density anomalies (regions of uniform density). Hence:

$\rho(\bm{x})=\sum_{j=0}^{M}\rho_{j}\mathds{1}_{j}(\bm{x})$ (4)

where $M$ is the total number of density anomalies and $\mathds{1}_{j}$ is the
indicator function of the $j^{th}$ anomaly, which is 1 inside the anomaly and
0 outside. For better clarity, we will reserve the index $j$ to describe any
of the $M$ anomalies, and the index $i$ to refer to any of the $N$ grid
elements. The index $j=0$ refers to the background density $\rho_{0}$, for
which the indicator function is $\mathds{1}_{0}=1$ over the whole body, and
which is also the density of the surface layer (that is, we impose
$\rho_{s}=\rho_{0}$). The parameter $\rho_{j}$ then represents the density
jump with respect to $\rho_{0}$, associated with the $j^{th}$ anomaly. This
way, heterogeneous interior models can be constructed by defining $M$ domains
for the density anomalies and assigning the corresponding density jumps. As in
Li et al. (2017) and Giraud et al. (2021), we don’t allow for overlap of
anomalies, meaning that at each point inside the body only one of the
$\mathds{1}_{j}$ apart from $\mathds{1}_{0}$ can have a value of 1. This is
handled in our algorithm on a "last-in-first-out" basis, meaning that in case
of overlapping of anomalies the $\mathds{1}_{j}$ of highest $j$ index is
assigned a non-zero value, while all the other indicator functions (save for
$\mathds{1}_{0}$) are set to 0.

###  2.2 Density function estimation

The linear relation of Eq. 3 forms an under-determined system that can be
inverted to obtain a vector of local densities
$\bm{{\rho}}=\left(\rho_{s},\rho(\bm{x}_{1}),...,\rho(\bm{x}_{N})\right)$ over
all the discretization elements. However, doing so without any regularization
terms or depth-weighting of the partials will give a density solution with
anomalies concentrated at the surface of the body, where the strength of the
partials is larger (Park et al., 2010). Depth-weighting in such an
underconstrained case was explored in several Earth local-gravity inversion
studies (e.g. Li and Oldenburg, 1998), but our attempts to extend it to the
global gravity inversion of small bodies were so far unsuccessful. On the
other hand, the regularizing assumption of piecewise-constant density (Eq. 4)
allows to substantially reduce the number of parameters to estimate, at the
cost of a limitation in the range of possible models. With these assumptions,
the $M$ density jumps of each anomaly can be estimated, along with the
background density, in place of the $N+1$ densities of each discretization
element (where $M<<N$), as in Takahashi and Scheeres (2014a) for instance.
This linear estimation problem is generally over-determined, as long as the
number of zones $M$ is selected to be below $(l_{\max}+1)^{2}$ (Scheeres et
al., 2000), and can be solved via weighted least-squares. However, since the
indicator functions $\mathds{1}_{j}$ are fixed in the estimation, the solution
will strongly depend on their value, and therefore on the initial subdivision
of the body.

####  2.2.1 Level-set method

In order to reduce the influence of the initial zones subdivision on the
density solution, the shape of the anomalies is here estimated along with
their density contrast. The shape derivative is computed via the level-set
method, as in Giraud et al. (2021). The boundary of each anomaly is defined
implicitly as the 0 level-set of a scalar function
$\phi:\mathbb{R}^{3}\mapsto\mathbb{R}$ (Figure 1). Clearly there are infinite
level-set functions $\phi$ providing the same anomaly boundary, but commonly
$\phi$ is taken to be the signed-distance function from the boundary, in order
to improve the numerical stability of iterative methods (Osher and Fedkiw,
2003). We compute the signed-distance via the _scikit-fmm_
222https://github.com/scikit-fmm/scikit-fmm Python implementation of the fast-
marching method (Sethian, 1996).

The indicator function of the $j^{th}$ anomaly is then defined as:

$\mathds{1}_{j}(\bm{{x}})=H(\phi_{j}(\bm{{x}}))$ (5)

where $\phi_{j}(\bm{x})$ is the level-set function specific to the $j^{th}$
anomaly, and in our case $\bm{x}$ takes the discrete values of the centers of
the grid elements. $H$ is the Heaviside (or step) function, so that:

$H(\phi_{j})=\begin{cases}1&\text{if }\phi_{j}\geq 0\text{ (inside the
anomaly)}\\\ 0&\text{otherwise (outside the anomaly)}\end{cases}$ (6)

Figure 1: Level-set representation of the density. Left: values along the
$Y=0$ plane of 2 level-set functions ($\phi_{1}$ and $\phi_{2}$), obtained as
the signed-distance function of the Stanford bunny and the Utah teapot. The
black cubes indicate the narrow-band region. Right: density model produced by
these 2 level-set functions, assuming as density contrasts
$\rho_{1}=\qty{600}{kg/m^{3}}$ and $\rho_{2}=\qty{-200}{kg/m^{3}}$ , and as
background density $\rho_{0}=\qty{1000}{kg/m^{3}}$ .

By combining Eqs. 3, 4, and 5, we therefore obtain an expression of the Stokes
coefficients in terms of the density jump of each anomaly and of the values of
the corresponding level-set function at the grid nodes:

$\mathcal{O}_{lm}=\frac{1}{\mathcal{M}}\left[\rho_{s}\mathcal{V}_{s}\mathcal{O}_{lm}^{s}+\sum_{i=1}^{N}\sum_{j=0}^{M}\rho_{j}H(\phi_{j}^{i})\mathcal{V}_{i}\mathcal{O}_{lm}^{i}\right]$
(7)

where, for ease of notation, $\phi_{j}^{i}$ denotes $\phi_{j}(\bm{x}_{i})$,
that is the value of the $j^{th}$ level-set function at the location of the
$i^{th}$ grid element. In order to generalize our conclusions to bodies of any
shape and size, we normalize the signed-distance function, dividing it by the
minimum grid resolution.
For practical applications where an analytical derivative is advantageous, a
differentiable approximation of the Heaviside function is generally preferred,
such as:

$H(\phi_{j})=\begin{cases}1&\text{if }\phi_{j}>\tau\\\
\frac{1}{2}+\frac{\phi_{j}}{2\tau}+\frac{1}{2\tau}\sin{\left(\frac{\pi\phi_{j}}{\tau}\right)}&\text{if
}\left|\phi_{j}\right|\leq\tau\\\ 0&\text{if }\phi_{j}<-\tau\end{cases}$ (8)

where $\tau$ is a distance defining the maximum distance from the anomaly
boundary where the level-set function is allowed to vary (black cubes in
Figure 1). This process is known as _narrow-banding_ (Adalsteinsson and
Sethian, 1995). The parameter $\tau$ is usually taken to be a small multiple
of the grid resolution, expressed as the minimum grid element size
$\Delta\bm{x}$: here nominally $\tau=1.5\min(\Delta x,\Delta y,\Delta z)$.
From the continuous approximation of Eq. 8, it is easy to compute the partial
derivatives of the anomaly indicator function with respect to the values of
the level-set function at the grid nodes. However, this approximation proved
to be less stable in our iterative optimization, possibly because of vanishing
partials for cells near the edge of the narrow-band increasing the condition
number of the observation matrix. For this reason, while still opting for the
narrow-banding philosophy due to its computational efficiency, we use the non-
differentiable expression for the Heaviside function (Eq. 6), with approximate
partials coming from the signed-distance property of $\phi_{j}$:

$\frac{\partial\mathds{1}_{j}}{\partial\phi_{j}}=\frac{1}{\left|\phi_{j}\right|+\xi}$
(9)

$\xi$ being a small positive term that limits the value of the partials when
$\phi_{j}\mapsto 0$. This expression reflects the fact that if the value of
$\phi_{j}^{i}$ changes by $-\phi_{j}^{i}$, the border of the $j^{th}$ anomaly
shifts to the center of the $i^{th}$ cell, meaning that $\mathds{1}_{j}^{i}$
varies by $\pm 1$. Small values of $\xi$ provide a better approximation of the
numerical partials, but on the other hand lead to high condition numbers of
the observation matrix. Here we take $\xi=0.1$, which balances these two
effects.

The partials of the Stokes coefficients with respect to the level-set values
of each anomaly and at the location of each grid element are then computed via
the chain rule:

$\frac{\partial\mathcal{O}_{lm}}{\partial\phi_{j}}=\frac{\partial\mathcal{O}_{lm}}{\partial\rho}\rho_{j}\frac{\partial\mathds{1}_{j}}{\partial\phi_{j}}$
(10)

where, from Eq. 3 and as in Scheeres et al. (2000):

$\frac{\partial\mathcal{O}_{lm}}{\partial\rho}=\frac{1}{\mathcal{M}}\sum_{i=1}^{N_{j}}\mathcal{V}_{i}\mathcal{O}_{lm}^{i}$
(11)

with $N_{j}$ the number of grid elements contained in the $j^{th}$ anomaly.

####  2.2.2 Iterative least-squares solution

Estimating the shape of the anomaly along with its density jump leads to a
non-linear problem. The target parameters for each anomaly are its density
jump $\rho_{j}$ and the nodal values of the associated level-set function,
i.e. the set $\\{\phi_{j}^{i}\\}_{i=1,...,N}$, although in our case $i$ only
takes the index-values of the grid elements within the narrow-band of the
anomaly. We minimize the residuals between the observed Stokes coefficients
and those computed from Eq. 7 in a least-squares sense. Given the partial
derivatives defined above, equation 7 is linearized in a neighborhood of the
current values of the target parameters. If $\bm{y}$ is the vector of
$P\leq(l_{\max}+1)^{2}$ residuals and
$\bm{t}=\left(\rho_{0},...,\rho_{M},\phi^{1}_{1},...,\phi^{N}_{1},\phi^{1}_{2},...,\phi^{N}_{M}\right)$
the vector of $K$ solve-for parameters, then:

$\bm{y}=J\bm{\delta t}+\bm{\varepsilon}$ (12)

$\bm{\delta t}$ being a small correction to the parameters, and $J$ the
$P\times K$ Jacobian matrix, where the term $J_{pk}$ is computed from Eq. 11
or Eq. 10 based on whether the $k^{th}$ parameter is a density or a level-set
correction, respectively. The measurements errors $\bm{\varepsilon}$ are
generally correlated, but the associated covariance matrix
$\Omega\in\mathbb{R}^{P\times P}$ on the Stokes coefficients can be assumed to
be known, as it is an output of a radio-science campaign. These correlations
can be accounted for in a generalized least-squares approach (e.g. Björck,
2015): if $\Omega$ is assumed to be symmetric positive definite, it has a
unique Cholesky factorization $\Omega=WW^{T}$, with $W\in\mathbb{R}^{P\times
P}$ nonsingular. By pre-multiplying both sides of Eq. 12 by $W^{-1}$, we then
transform the system so that the elements on the left-hand side have
associated errors which are uncorrelated and have the same variance. The
whitened residuals vector is then $\bm{y^{\prime}}=W^{-1}\bm{y}$ and the
corresponding observation matrix $J^{\prime}=W^{-1}J$. In the following we
drop the prime symbol, so that $\bm{y}$ indicates residuals with uncorrelated
and homoscedastic errors and $J$ their observation matrix. In the simulation
tests of Section 3 we assume no correlations in the synthetic
$\mathcal{O}_{lm}$, so that the problem is reduced to a weighted least-
squares, where $W=\operatorname{diag}(\sigma_{1},...,\sigma_{P})$ is a
diagonal matrix with the standard deviations of the measurement errors as
elements.

The objective function ($\Psi$) we wish to minimize is then that of ordinary
least-squares, with additional penalty terms:

$\Psi(\bm{\delta t})={\lVert\bm{y}-J\bm{\delta
t}\rVert}_{2}^{2}+\lambda^{2}{\lVert\bm{\delta
t}\rVert}_{2}^{2}+\Psi_{c}(\bm{\delta t})$ (13)

where ${\lVert\cdot\rVert}_{2}$ is the 2-norm. The second term in the
expression of $\Psi$ is the Tikhonov regularization term (e.g. Björck, 2015),
which penalizes large corrections on the parameters. It is necessary in order
to render the method robust with respect to noise in the data, and its weight
$\lambda$ can be assigned empirically, or via generalized cross-validation
(GCV). The last penalty term is a generic quadratic constraint, where
available. For example, one could search for solutions close to a given
reference (_a priori_) density distribution. Then, if $\bm{t^{\prime}}$ are
the values of the solve-for parameters in the reference model, and $\bm{t}$
their values at the current iteration, the penalty term could be:

$\Psi_{c}(\bm{\delta t})=\nu^{2}{\lVert\bm{t^{\prime}}-\bm{t}-\bm{\delta
t}\rVert}_{2}^{2}$ (14)

where $\nu$ is the weight of the constraint, which in this case favours
corrections bringing the current parameter values closer to those of the
reference model. In practice, these quadratic penalties can be treated as
supplemental observations and added as extra rows to Eq. 12. We use a Gauss-
Newton method to iteratively adjust the target parameters, so that at each
iteration the correction $\bm{\delta t}$ estimated from the linearized system
is simply added to the vector of parameters. The system is inverted using the
SciPy implementation of the LSQR algorithm
333https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsqr.html,
which can handle the sparse matrices resulting from the addition of penalties.
Alternatively, the _RidgeCV_ function of the scikit-learn Python library
(Pedregosa et al., 2011) is used to solve the Tikhonov-regularized problem by
means of GCV.

More often than not, the objective function turns out to be non-convex.
Similarly to restarts in gradient-descent optimization, we mitigate the
problem of early convergence to local minima by periodically amplifying the
estimated $\bm{\delta t}$ so that the maximum level-set correction is larger
than the narrow-band parameter $\tau$. For the same reason, and because we
assume no prior knowledge about the interior of the body, we nominally start
from initial models where the different (usually 3 or more) level-set
functions take positive values all over the body.

Figure 2 shows a schematic of GILA. The orange section (left-hand side of the
figure) is the core part of the method and consists of the iterative least-
squares inversion, which provides a single density solution given the inputs.
The latter are represented in blue, in the top-right region, and include: the
observables, namely the Stokes coefficients $\mathcal{O}_{lm}$; the shape
model of the body, in the form of a triangular mesh; eventual quadratic
constraints, which would help restrict the range of possible solutions; a set
of initial density distributions, as running the inversion algorithm with
different _a priori_ is advised given the non-convexity of the problem. For
the moment, GILA does not support running multiple inversion loops in
parallel. A single run of the iterative least-squares until convergence
requires on average 3 minutes on a laptop with CPU Intel i7-12700H (14 cores,
up to 4.70 GHz) for an interior grid of 50 elements per dimension, while the
time for the computation of the $\mathcal{O}_{lm}$ (which are stored once
computed) depends on the resolution of the mesh and the degree of the gravity,
but is itself at most a few minutes in the cases considered here. If several
initial models are used, the outputs will be a set of solutions which may fall
into different regions of the solution space. In order to distinguish between
mass distributions which agree with the data within the noise but are widely
different from each other, we use global quantities independent from the input
measurements, such as the moments of $\rho(\bm{x})$ (see Eq. 18). Then, as
displayed in the green section of the figure, the set of solutions can be
subdivided into a small number of interior families via a clustering
algorithm, and each family can be characterized on the basis of its average
distribution and the spread of the density values at each grid element.

Figure 2: Pipeline of the inversion algorithm

##  3 Synthetic retrievals

In this section, we assess the ability of the method to retrieve a known
density distribution from synthetic gravity measurements. In each case, the
known density model is used for the computation of synthetic observed Stokes
coefficients $\mathcal{O}_{lm}$ according to Eq. 3. Then, the iterative least-
squares algorithm is run until either the reduced chi-squared of the residuals
(as defined in Section 3.1.1) reaches a minimum threshold of 0.1 and the model
has stopped evolving (after a warm-up of 500 iterations), or a set limit of
maximum iterations is reached. Generally, for the first 100 iterations the
density contrasts are fixed, and only the background density and the level-set
parameters (meaning the shapes of the anomalies) are adjusted. This appears to
help avoid early convergence of the solution, although it makes the iterations
more dependent on the initial contrasts assigned to each density anomaly.
We start with noise-free measurements, in order to gauge how the method is
robust with respect to the non-uniqueness of the solution and the non-
convexity of the objective function. We then present more realistic scenarios
where the synthetic data are contaminated with Gaussian noise, according to
different noise profiles, with the target of testing the behaviour of the
method with respect to another characteristic of ill-posed problems, namely
their instability.

The absolute uncertainty of the Stokes coefficients measured in a radio
science campaign is usually higher for high degrees, since their lower signal-
to-noise ratio means their effects on the noisy trajectory reconstruction are
harder to decorrelate than for low degrees. This could however not always be a
fair prediction. For example, when the coefficients are obtained via optical
tracking of orbiting dust particles they might display a flat uncertainty
profile as a function of the degree (Chesley et al., 2020). Still, the
standard deviations of the synthetic measurements errors are here assigned
following a noise profile varying as a function of the harmonic degree and
based on the exponential expression of Tricarico (2013), itself derived from
the uncertainty spectrum of the gravity coefficients estimated for Eros using
NEAR radio-tracking data (Miller et al., 2002). The noise spectrum is then
computed as:

$\sigma(l)=\alpha
10^{\beta(l-l_{\max})}\sqrt{\frac{\mathcal{S}_{l_{\max}}}{2l_{\max}+1}},$ (15)

where $l_{\max}$ is the maximum harmonic degree of the measured (or simulated)
Stokes coefficients and $\mathcal{S}_{l_{\max}}$ is their power spectrum (sum
of their squares) at the degree $l_{\max}$. The parameter $\alpha$ indicates
instead the relative strength of $\sigma(l_{\max})$ compared to
$\mathcal{S}_{l_{\max}}$. We therefore talk of $100\%$ noise to indicate a
case where $\alpha=1$, meaning that most of the coefficients of degree
$l_{\max}$ are consistent with 0. Similarly, 50%,10%, 1% noise values
correspond to $\alpha$ values of 0.5, 0.1, and 0.01, respectively. The
coefficient $\beta$ is the slope parameter, indicating how fast the
uncertainty decreases for degrees lower than $l_{\max}$. Figure 3 shows the
power spectrum of measured coefficients and their uncertainties for 5
different real datasets, as well as the $\beta$ value approximating their
noise profile. Degree-1 values and uncertainties are not provided in the real
datasets because assumed to be known and equal to 0, meaning that the
potential is expressed in a frame with its origin at the center of mass of the
body. The uncertainties resulting from Eq. 15 tend to underestimate the GM
uncertainty, which is why we multiply the degree-0 $\sigma$ by a factor of 10.
The simulated noise profile also underestimates the errors of the Dawn
datasets at low degrees. As in Tricarico (2013), in all our simulations we
assume $\beta=1/3$, which is roughly the value fitting the NEAR measured
uncertainties, and sits between those of the OSIRIS-REx spacecraft-field
uncertainties and those of the Dawn datasets.

Figure 3: Power spectrum of measured coefficients (solid lines) with
corresponding uncertainties (dashed lines) for the Eros gravity field of NEAR
(Miller et al., 2002), the two OSIRIS-REx gravity solutions for Bennu
(Scheeres et al., 2020), and the two Dawn gravity solutions for Vesta
(Konopliv et al., 2014) and Ceres (Konopliv et al., 2018). The dotted lines
represent the noise profile obtained from Eq. 15, assuming $l_{\max}$ to be
right before the eventual intersection between the solid and dashed lines. The
values of $\beta$ are computed as the average slope of the dashed lines over
the range of degrees where they show log-linearity.

###  3.1 Inversion quality metrics

The estimated density model is evaluated in terms of its accuracy in both
fitting the synthetic measurements within the given noise and reproducing the
ground-truth density distribution. Specifically, we test two different quality
metrics of the inversion: the reduced $\chi^{2}$ statistic of the whitened
measurements residuals and the Pearson correlation between the estimated and
the true density distributions. We look at the correlation because it appears
to be less model-dependent compared to other measures of the similarity
between interior models, such as the RMS of the density errors at each grid
element with respect to the ground-truth, or the relative error on the total
mass of the anomalies. However, of all these metrics only the reduced
$\chi^{2}$ is applicable to real data analysis, and while it can provide an
indication of good fit, the corresponding solution is not guaranteed to
correctly represent the true interior of the body.

####  3.1.1 Reduced $\chi^{2}$

The reduced $\chi^{2}$ statistic, denoted as $\chi_{P}^{2}$, is computed as
the ratio of the weighted residuals sum-of-squares (RSS) and the number of
measurements (P), that is:

$\chi_{P}^{2}=\frac{\sum_{p}{y_{p}^{2}}}{P}$ (16)

Values of $\chi_{P}^{2}$ close to 1 indicate a good fit consistent with the
given measurement weights, while $\chi_{P}^{2}>>1$ suggests an inaccurate fit
and $\chi_{P}^{2}<<1$ an overfitting of the data (where the model also fits
the noise). From simulations we observed that the latter case, that of
overfitting, is very rare for GILA, due to the values of the objective
functions being limited by the model error. This error is generally related to
the inability of the inversion model to reproduce exactly the ground-truth
with its parameters, but in our case an even larger contribution is due to the
fact that the iterative inversion algorithm is rarely able to reach exactly
the global minimum, but instead convergence happens at best in a neighborhood
of the optimum. This is justified by the non-convexity of the objective and
possibly by the partials of Equation 9 approximating less well the gradient of
the objective function around a minimum.

In order to predict the effect of this inaccurate convergence on the final
residuals, we introduce the concept of model resolution. The grid
discretization of the interior implies a non-dense image of the Jacobian
matrix $J$, since the smallest possible variation in the density model (e.g.
expanding an anomaly by a single grid element) results in a finite variation
in the predicted measurements. We characterize this model resolution by
looking at the columns of $J$ corresponding to the level-set parameters. For
each row of the matrix, the maximum of its absolute values over these columns
represents the minimum resolution of the model for the corresponding
measurement. This is because the level-set partials in $J$ are the
measurements variations per unitary change of the level-set function, and
given our normalization of $\phi_{j}$ by the grid-size, this unitary change is
a change of 1 grid element. Trivially, the model resolution thus computed
decreases with the density contrast of the estimated anomalies (see Eq. 10).
It also decreases with the decrease of the grid elements’ size, itself tied
however to higher computation times and a higher degeneracy of the problem.
Figure 4 shows that with 50 grid elements per dimension, the model resolution
is higher in magnitude than the uncertainty of the measurements for low
degrees, and drops quickly for higher degrees, while the noise profile
increases. The same behavior is observed for 100 grid elements per dimension,
but this time the intersection between the model resolution and the data noise
happens at a lower gravity degree. This trend and the gravity degrees at which
the model resolution intersects the synthetic uncertainties are similar for
bodies of different shapes and sizes, although the absolute values of the
model resolutions are body-dependent.
In simulations, we often test GILA on noise profiles where for lower degrees
the measurements $\sigma$ might be smaller than the model resolution.
Therefore, even converged models which differ by as little as 1 grid element
from the truth might show large values of $\chi_{P}^{2}$. To better assess the
convergence of the solution, we replace the measurements weights with the
model resolution for data points where the latter is larger than the input
measurements uncertainty. Therefore, the reduced $\chi^{2}$ statistic becomes
a measure of the agreement between the predicted and the observed measurements
when both the data uncertainty and the model resolution are taken into
account. Of course, for real data applications, the $\chi_{P}^{2}$ with only
the data uncertainty should also be considered. Until further improvements to
GILA, if the residuals are limited by the model resolution the options are to
increase the grid resolution (where feasible, see effect of doubling the
resolution in Figure 4) or to accept a data fit limited by small errors in the
model.

Figure 4: Comparison between model resolution and measurements accuracy for a
Bennu-shaped body. The blue shaded areas are envelopes of the model
resolutions for 100 random density distributions, using a grid of 50 (darker
area) or 100 (lighter area) elements per dimension. Solid and dashed orange
lines indicate the RMS of the synthetic coefficients and corresponding errors,
respectively, for $l_{\max}=5$ (darker lines) and $l_{\max}=10$ (lighter
lines). The measurements values are computed for a homogeneous distribution,
and their uncertainties follow the noise profile of Equation 15 with
$\alpha=1$ (i.e. in the 100% noise case) and $\beta=1/3$.

####  3.1.2 Density correlation

The Pearson correlation between the ground-truth and retrieved density
distribution is computed as (e.g. Lee Rodgers and Nicewander, 1988):

$\rho_{CORR}=\frac{\sum_{i}{\left(\rho_{i}^{(est)}-\overline{\rho}^{(est)}\right)\left(\rho_{i}^{(sim)}-\overline{\rho}^{(sim)}\right)}}{\sqrt{\sum_{i}{{\left(\rho_{i}^{(est)}-\overline{\rho}^{(est)}\right)}^{2}}\sum_{i}{{\left(\rho_{i}^{(sim)}-\overline{\rho}^{(sim)}\right)}^{2}}}}$
(17)

where, $\bm{\rho}^{(est)}$ is the vector of estimated density values over each
element grid, and $\bm{\rho}^{(sim)}$ the ground-truth density values
interpolated over the inversion grid, while $\overline{\rho}^{(est)}$ and
$\overline{\rho}^{(sim)}$ are their mean values.

###  3.2 Target models

We mostly limit our assessment of the method to three reference models,
denoted with numbers 1,2,3 with increasing level of complexity. Specifically:
Model 1 is composed of a single off-centered positive anomaly of regular
shape, namely a triangular prism, with a density contrast of 600kg/m^3; Model
2 has a single positive anomaly of 600kg/m^3 centered at the origin and of
complex shape, based on the Stanford dragon; Model 3 has multiple anomalies of
complex shapes and different density contrast, in both sign and magnitude (2
positive ones of 600kg/m^3, and 2 negative ones of -200kg/m^3, and
-600kg/m^3). All three models are piecewise-constant but approximating a wide
range of scenarios where the GILA could be applied. Degenerate cases, such as
those with concentric shells, are not discussed here as their severe non-
uniqueness hinders any possible conclusions about the performance of GILA.

Figure 5 shows the three reference models and the corresponding retrieved
models at convergence when fitting synthetic Stokes coefficients up to degree
11 for a Bennu-shaped body. These plots are meant to show the performance of
the model in a purely ideal case, assuming a 1% noise profile for the
measurements weights but no noise perturbation on the measurements and
consequently no regularization term in the objective function, which is only
composed of the data misfit. Moreover, the shape model used in these
simulations is the same as that of the inversions. As density and shape
parameters are highly correlated, even at convergence the detected anomalies
have density jumps which are not exactly those of the true model, while still
agreeing with the data at the noise level, as suggested by the value of the
metric $\chi_{P}^{2}\sim 1$. Nevertheless, Figure 5 proves that GILA is able
to retrieve number, sign, and approximate shape of multiple density anomalies
without injection of any prior knowledge about the ground truth. This
convergence to the true model, while seemingly clashing with the non-
uniqueness of the true inversion, is justified by the reference models being
composed of isolated anomalies with relatively high density contrast, by the
high resolution of the gravity measurements, and by the forward and inverse
discretization grids being identical.

Figure 5: Three synthetic retrievals of the density models at the bottom right
of each subpanel, corresponding to Models 1, 2, and 3, from top to bottom. The
line plots represent the values of $\chi_{P}^{2}$ over the least-squares
iterations, with red circles indicating specific steps of the convergence
displayed in the 3D plots, in the same order. The right-most point and its
related 3D plot correspond to the solution at convergence. The 3D plots are
obtained via the marching-cubes algorithm (as implemented in the scikit-image
Python library). The density colorscale is shared between all 3D plots of a
given subpanel. The 3D meshes employed for reference Model 3 are Spot and Bob
(courtesy of Prof. Keenan Crane), #3DBenchy (http://3dbenchy.com/), and
Suzanne (Blender Foundation).

###  3.3 Robustness

In this section, we build towards an assessment of the GILA over more
realistic inversion scenarios, and at the same time attempt to delineate the
limits of applicability of the method. To this end, we modify one or more of
the scenarios of Figure 5 by one property at a time, and check how this
impacts the inversion quality metrics. The inversion settings varied here are:
the degradation of the synthetic measurements, both in their resolution
(Section 3.3.1) and level of noise contamination (Section 3.3.2); the shape
and size of the body (Section 3.3.3); the error on the body shape (Section
3.3.4); the resolution of the internal discretization grid (Section 3.3.5);
the size, depth, and density contrast of the anomaly (Section 3.3.6); the
initial density model in the iterative fit (Section 3.3.7).

####  3.3.1 Gravity resolution

Here we test the impact of the gravity field resolution (meaning the maximum
degree of the observed Stokes coefficients) on the accuracy of the retrieved
interior model. The maximum angular resolution that can be expected from a
truncated spherical harmonics series can be as a rule of thumb taken to be
that of the half-wavelength of its maximum degree, hence given by
$180^{\circ}/l_{\max}$ (e.g Jekeli, 2007). Therefore, a degree-11 gravity
would be unable to discern between anomalies at a resolution below about
$16^{\circ}\times 16^{\circ}$ in latitude and longitude. As for the radial
resolution, it is virtually limited by our interior grid discretization,
although degree-depth relationship can be found in the literature proving that
lower degrees are more sensitive to deeper anomalies (Goossens and Smith,
2023). We therefore expect models obtained with lower degrees to have a larger
uncertainty associated with the estimated anomalies shape. Similarly, since
the number of measurements decreases when $l_{\max}$ decreases while the
number of parameters remains roughly the same, we expect that at each
iteration the null space of the locally linearized problem has larger
dimensions for lower $l_{\max}$, meaning that a wider range of corrections
could provide the same decrease in the objective function.

Figure 6 shows both visually and through the density correlation the behaviour
of the GILA single-run solutions with the decrease of the gravity field
resolution, for the three reference models. Apart from the measurements
resolutions, all settings are maintained from the cases shown in Figure 5,
meaning we still consider a Bennu-shaped body and measurements unperturbed by
noise. The $\chi_{P}^{2}$, not shown here, are all close to 1.

Figure 6: Density solutions (3D plots) and corresponding Pearson correlations
(bar plots) with the true model as a function of the maximum gravity degree of
the synthetic noise-free measurement set. The density colorscale is shared
between all 3D plots of a given subpanel

It can be seen that for all three models, degree 7 provides a quality of the
reconstruction comparable to that obtained from higher degrees. Degree 5 seems
to be the threshold for an accurate detection of these models, save for Model
2, which compared to the other two has heterogeneities mostly concentrated
near the center of the body.
From these plots, we can empirically define a correlation threshold of 0.2 to
indicate a successful reconstruction, with values higher than 0.6 representing
an accurate retrieval of the target model. However, as the 3D plots show, even
for solutions with correlations lower than 0.2 the estimated density
distribution may not be extremely different from the truth, at least in terms
of location and sign of the anomaly.

####  3.3.2 Noise

In this section we examine the effect of data noise on the estimated density
distribution. To this end, the measurements are perturbed with zero-mean
Gaussian noise, sampled randomly from the profile of Eq. 15. Small
perturbations in the measurements due to noise may result in vastly different
least-squares solutions $\bm{\delta t}$, given the inherent instability of the
problem. In a singular-value-decomposition (SVD) formulation of the least-
squares solution, this effect can be traced to small singular values of the
observation matrix $J$ (compared to the highest singular-value) that divide
the residuals vector $\bm{y}$, greatly amplifying the errors in the data
(Björck, 2015). To mitigate this instability, we add a classical Tikhonov
regularization term to the least-squares objective function, minimizing the
norm of the correction $\bm{\delta t}$ (second term in Equation 13). We select
an empirical regularization parameter of $\lambda=3$ for both density and
level-set corrections. This solution is equivalent to the classical least-
squares solution with _a priori_ knowledge, where the _a priori_ central
values of all parameter corrections in $\bm{\delta t}$ are 0 and their _a
priori_ covariance is a diagonal matrix with diagonal elements equal to
$\lambda$, corresponding to _a priori_ formal errors of $1/\sqrt{\lambda}$.
Figure 7 shows the dependence of the solutions on the values of the
regularization parameter, justifying our choice of $\lambda$. Each row in the
plot corresponds to solutions for the Model 1 ground-truth, two different
bodies (Bennu and Eros), and 3 different levels of noise perturbation (0%, 1%,
and 100%). For $\lambda<<1$ the solution becomes equivalent to the least-
squares solution without regularization. The figure also reports the
generalized cross-validation solutions for each case, and their agreement with
the $\lambda=3$ solutions, coupled with the heavier computation efforts of the
GCV, leads us to prefer the use of a fixed regularization parameter.

Figure 7: Correlations of simulation solutions as a function of the
regularization settings, for different noise levels, bodies, and gravity
resolutions. The first column is the classical Tikhonov regularization for
$\lambda$ in the range [$10^{-7},10^{4}$], with the vertical dashed line
covering values for $\lambda=3$. The second column is the Tikhonov solution
with $\lambda$ selected by generalized cross-validation.

The effects of data noise using this regularization approach are presented in
Figure 8, which shows the evolution of the density correlations for different
values of $\alpha$ in Eq. 15 (corresponding to measurement noise levels of 0%,
1%, 10%, and 100%), different resolutions of the synthetic gravity
measurements ($l_{\max}$ of 3, 7, and 11), and for the 3 target models of
Section 3.2. Otherwise, all simulation settings are the same as in that
section. The set of points labeled as $\lambda=0$ are the solutions using
noise-free data and no regularization penalty, which are in general different
from the Tikhonov 0%-noise solutions due to the damping of high-frequency
contributions coming from the regularization term. These $\lambda=0$ solutions
are only shown for the noise-free case, given their poor handling of noisy
data (see Figure 7). For noise levels higher than 0%, each point in the plot
represents a solution obtained from measurements sampled randomly within the
noise profile, meaning each point is obtained from a different random
perturbation of the simulated measurements. The degradation of the retrieved
model coming from the contamination of noise in the data is considerable, with
most of the Model 1 and 2 solutions in the degree-3 case displaying 0 or
negative correlations. However, even with 100% noise, for degrees 7 and 11 the
estimated density distributions preserve a strong correlation (>0.2) with the
ground truth.

Figure 8: Density correlations of GILA solutions as a function of the level of
data noise contamination, with $\lambda=0$ representing the non-regularized
inversions. The points at each noise level represent $\rho_{CORR}$ of
solutions for different sets of measurement perturbations sampled within the
given noise profile. For each noise level, the box-plots bounds mark the first
and third quartile, while their internal line is the median. The length of the
whiskers is 1.5 times the interquartile range. For each gravity degree, the
lines connect the median values of the density correlations, and the shaded
regions cover their range of values.

####  3.3.3 Shape and size of the body

So far, we have limited our analysis to a Bennu-shaped body. Now we examine
the performance of the model for bodies with different shapes and sizes. The
settings are the same as in Section 3.3.2, with Model 1 as the ground truth
and varying gravity resolution and level of noise perturbation, but the input
shape model of both the forward and inverse computation is different. Figure 9
shows the correlations for the interior densities estimated using shape models
of Eros, Kleopatra, Phobos, and comet 67P/Churyumov-Gerasimenko, as a function
of the maximum degree of the input gravity measurements and for 4 different
noise levels between 0% and 100%. The size of these bodies in terms of maximum
radius ranges from about 2km for comet 67P to about 60km for Kleopatra,
compared to the 290m radius of Bennu. This size difference is not expected to
impact the performance of GILA given the normalization of the level-set
functions. Indeed, for degrees 5 and higher the correlation of the retrieved
model with the ground truth is independent of the body. For degree 3, the
reconstruction accuracy is considerably worse for Phobos, as was the case for
Bennu (Figure 8). We think this could be a consequence of their rounded shapes
compared to others of more irregular shape, for which the anomaly would tend
to be closer to the surface, where the sensitivity of gravity measurements is
higher.

Figure 9: Density correlation between the Model 1 ground-truth and GILA
solutions for simulations over different bodies (one per panel). The points at
each noise level represent solutions for different sets of measurement
perturbations sampled within the given noise profile. For each gravity degree,
the box-plots bounds mark the first and third quartile, while their internal
line is the median. The length of the whiskers is 1.5 times the interquartile
range. For each noise level, the lines connect the median values of the
density correlations, and the shaded regions cover their range of values. The
3D mesh plots show the shape of the body.

####  3.3.4 Shape error

In all previous examples, we assumed a perfect knowledge of the shape, meaning
that the same surface mesh used for the data simulation was also used in the
inversion. In reality, the shape model of a body has a certain uncertainty
associated with it, both coming from the reconstruction errors and from the
polyhedral approximation. For Phobos, the 1$\sigma$ error on the shape model
can reach 40m, corresponding to about 0.3% of its radius (Willner et al.,
2014). For Bennu, the shape error relative to the radius is of similar
magnitude, at about 0.25% (Asad et al., 2021). Here, we simulate an uncertain
knowledge of the shape of the body by perturbing the location of the mesh
vertices by a random offset along their normals. The offset for each point is
sampled from a normal distribution of given standard deviation. The perturbed
model is then used in the inversion, fitting data generated with the original
shape model. The results are shown in Figure 10 for a Model 1 ground-truth and
noise-free measurements, and shape-offset standard deviations ranging between
0.001 and 1% of the body radius. As a perturbation on the shape model becomes
effectively a perturbation on the residuals, Tikhonov regularization is here
necessary to stabilize the inversion.

Overall, GILA solutions show strong correlation with the ground-truth even in
the presence of moderate error on the shape. As in Section 3.3.3, the overall
behaviour is better for the more irregular bodies, with high correlations even
at 1% mesh perturbation levels and down to degree-3 gravity. For Phobos, as
for Bennu, the 0.3% level appears to be at the edge of detectability for this
simulation case.

Figure 10: Influence of shape error on simulated retrievals for Model 1 and
noise-free measurements. Each panel displays a heatmap of correlation values
between the estimated density distribution and the ground-truth for a specific
body shape. The x-axis represents the 1$\sigma$ mesh perturbation applied
before the inversion, and the y-axis the maximum gravity degree of the
synthetic measurements.

####  3.3.5 Grid discretization

Since in the previous cases the interior grid used in simulation was the same
as that used in the inversion (same cell size), the target anomaly could in
theory be perfectly represented by the inversion model. Of course, in a real
case the interior discretization would be another source of error in the
reconstruction of the real anomaly. We find however this error to be
negligible for GILA, as shown in Figure 11. This heatmap shows the correlation
between the inverted model and the ground truth as a function of the
resolution of the grid used for the computation of synthetic measurements (in
terms of number of grid cells per dimension) and that used in the inversion.
The settings are those of Section 3.2, with Model 1 as a ground-truth. The
correlation is strong for all the cases, with minimum values around 0.6 for
coarser models in both simulations and estimation.

Figure 11: Heatmap of correlation between reconstructed and ground-truth
model. On the x-axis is the resolution of the interior grid used for the
inversion, and on the y-axis that of the grid used in simulation, both in
terms of number of cells per dimension.

####  3.3.6 Size, depth, magnitude of the anomaly

In this section we explore the range of detectability of an anomaly in terms
of its size, the magnitude of its density contrast, and its depth within the
body. To this end, we consider ground-truth distributions based on Model 1,
where this time the prismatic anomaly is randomly modified in the location of
its center (along the x-axis, with x-coordinate in the range [0, $r_{0}$]),
its density contrast (in the range [$0.1\rho_{0}$, $1.1\rho_{0}$]), and its
size (scaling factor in the range [$0.06r_{0}$, $r_{0}$]), all 3 parameters
being sampled from uniform probability distributions. The solutions
corresponding to each random target model are shown in Figure 12, for 3 levels
of gravity resolution and against the density contrast, size, and distance
from the body center of the prismatic anomaly in the true model. The synthetic
measurements are here perturbed by 100% zero-mean Gaussian noise. The size of
the anomaly is expressed in terms of its volume relative to
$\mathcal{V}_{\mathcal{B}}$, and the skewness towards lower values of its
distribution (while the sampling of the scaling factor was uniform) is due to
the prism not being fully contained inside the body for higher scaling
factors.
The curves represent order-2 polynomial regressions over the points in blue,
grey points being excluded from the fit. These grey points correspond to
ground-truth models where the volume of the anomaly relative to
$\mathcal{V}_{\mathcal{B}}$ is lower than $1/{l_{\max}}^{2}$. This threshold,
represented by grey boxes in the figure, is the fraction of the volume of a
sphere occupied by a tesseroid of $180^{\circ}/l_{\max}$ in latitude and
longitude, which is the half-wavelength spatial resolution of the
degree-$l_{\max}$ gravity field (Section 3.3.1). Accordingly, the second
column in Figure 12 shows that the density reconstruction becomes worse below
this threshold, with most solutions showing near-zero correlations with the
ground truth at low gravity degrees. The accuracy of the solution also worsens
for large sizes of the target anomalies, possibly because there the bulk
density (which we use as an initial value for the background density,
$\rho_{0}$) is further away from the background density of the true model.
Independently of the size and depth, the accuracy of the retrieved model is
higher for higher density contrasts of the target anomaly, although for high
resolutions of the gravity successful retrievals are common even for density
contrasts of 10% of the background density. The dependence of the correlation
on the depth of the anomaly, here represented by its distance from the center
of the Bennu-shaped body, is lower for high degrees. This is presumably
because due to the exponential noise profile, a higher resolution of the
gravity leads also to a lower relative uncertainty of the low-degree
coefficients, which are more sensitive to deep anomalies. Even at degree 3,
GILA seems to be able to retrieve anomalies deep inside the body, as long as
their size is above the gravity spatial resolution.

Figure 12: Scatter plot of density correlations of GILA solutions with their
respective ground-truth with a prismatic anomaly of random size, depth, and
density contrast. The x-axes represent: the density contrast of the target
anomaly relative to the background density $\rho_{0}$ (first column); the size
of the anomaly as a percentage of $\mathcal{V}_{\mathcal{B}}$ (second column);
the distance of the target anomaly from the center of the body (closely
related to its depth for a Bennu-shaped body) and relative to the radius of
the body (third column). The blue curves represent degree-2 polynomial
regressions over the blue points, the shaded areas being their 1$\sigma$ error
bounds. The grey points, excluded from the fit, represent solutions where the
volume of the target anomaly is lower than $1/{l_{\max}}^{2}$ times the body’s
volume (grey boxes in the second column). The rows represent the maximum
gravity degree of the synthetic Stokes coefficients, which are contaminated
with 100% Gaussian noise.

####  3.3.7 Initial model

The non-convexity of the objective function means that the solution is
strongly dependent on the starting point of the iterative algorithm, that in
our case is the initial density model. In a realistic scenario, the core
algorithm of GILA should be run several times starting from a wide range of
density models. This should provide a wider view of the space of possible
solutions, although a complete exploration of the solution space is hindered
by both computational efforts and biases in the estimation process which could
make some solutions more likely to appear than others. The aim of this section
is to provide an example of how changing the initial model can affect the
density solution. The target model is inspired by those proposed for Bennu by
the OSIRIS-REx team (Scheeres et al., 2020), which were obtained from the
inversion of the optical-tracking $\mathcal{O}_{lm}$ dataset and supported by
analytical and geophysical considerations. These models are composed of a
spherical core with a mass deficit of 6 to 16% of the total mass, an
equatorial bulge with a density 5 to 12% lower than the bulk density, and a
middle layer with a density 8 to 17% higher than the measured bulk density of
1190kg/m^3. Consequently, in this simulation campaign we consider as ground-
truth a model with a spherical cavity at the center of radius 116m (40% of the
total radius) and an equatorial ring with a negative density jump of
-400kg/m^3, for a background density of 1370kg/m^3, 15% higher than the bulk
density. We simulate Stokes coefficients up to degree 10, and contaminate them
with 100% noise using as $\beta$ in Eq. 15 the usual value of 1/3 to be
consistent with the previous sections, although Figure 3 shows that for the
Bennu particle field $\beta$ is closer to 0 (see Figure 15 and Sec. 3.4 for
discussion relative to $\beta$).
The synthetic measurements are then inverted multiple times, assuming in each
run a random initial model, based on variations on the generic initial model
employed throughout the previous sections and shown in Figure 5. The spheres
making up each of the 3 types of anomalies are now randomized in both their
number and locations, while their radius is still fixed to 20% of the body
radius. Moreover, the initial background density $\rho_{0}$ is itself sampled
from a uniform distribution between 50% and 150% of the bulk density, and the
3 initial density contrasts are sampled from a uniform distribution in the
range [-$\rho_{0}$, $\kappa\rho_{0}$] ($\kappa=1$ here), allowing for initial
models with void regions (large macroporosity).

The results are shown in Figure 13. For all the distributions the reduced
$\chi^{2}$ value is around 1, indicating a good fit with respect to the data
and model errors. As the distribution of the density correlations appears to
be multimodal, we separate the full set of solutions into 2 cluster, based on
a _k_ -means clustering444https://scikit-
learn.org/stable/modules/generated/sklearn.cluster.KMeans.html of the values
of the (normalized) principal moment of inertia along $z$ over all the density
models. The moment of inertia $I_{zz}$ is a global property of the density
distribution which is independent of the gravity data used to produce the
interior model, since only the difference between any 2 principal moments of
inertia is constrained by the gravity (Le Maistre et al., 2019). For this
reason, it can be used to differentiate between 2 models which both fit the
gravity measurements at the same level. More generally, the multiple solutions
of the gravity inversion can be grouped based on the $l^{th}$-order moments of
the density function, defined as (Jekeli, 2007):

$\mu_{abc}=\int_{\mathcal{V}_{\mathcal{B}}}x^{a}y^{b}z^{c}\rho(\bm{x})dV$ (18)

with $a,b,$ and $c$ integer indices such that $a+b+c=l$. The moments of order
$l$ are related to the Stokes coefficients of degree $l$, but not all the
$(l_{\max}+1)(l_{\max}+2)(l_{\max}+3)/6$ moments up to order $l_{\max}$ are
unambiguously determined by the ${(l_{\max}+1)}^{2}$ Stokes coefficients up to
degree $l_{\max}$ (Tricarico, 2013) and thus constrained by the gravity
measurements.

Figure 13 shows a representative model for each cluster, obtained as the
average of the density values at each cell over all solutions. The best set of
models, shown in blue, includes about 93% of the 146 solutions, and has a mean
correlation of about 0.1. While the value itself is low (lower than the
identified 0.2 threshold), it can be seen from the average model that in
general the solutions in this group do present a lighter outer region at the
equator and negative anomalies near the center, where the void core is
supposed to be. The remaining solutions, shown in orange, have a low moment of
inertia and negative correlation, both justified by the presence of positive
anomalies deep within the body which appear in the average solution for the
cluster. Overall, the concentration of most solutions around a single cluster
may suggest an incomplete exploration of the solution space, to be ascribed
either to a range of initial models not diverse enough, or to considerable
biases in the inversion process.

Figure 13: Density solutions for Bennu’s ring-core model from degree-10, 100%
noise data. Bottom left: distributions of the quality metrics and of the
normalized moment of inertia $I_{zz}$, with colors indicating the 2 families
obtained from clustering of the $I_{zz}$ values. Top right: averages of the
density distributions in each cluster (blue and orange boxes), all represented
as isosurfaces and slices along the $X=0$, $Y=0$, and $Z=0$ planes, and the
ground-truth model (grey box). The value of $\rho_{0}$ reported above the
average solutions is the median of the average density distribution, and the
density levels of the 6 isosurfaces are selected from _k-means_ clustering of
the deviations from this $\rho_{0}$ throughout the body. The colorbar is
relative to the three 3D plots.

###  3.4 Uncertainty estimation

We characterize here possible approaches to the estimation of errors
associated with a retrieved density distribution. The non-convexity of the
problem and the under-constrained approach we chose mean the formal errors
provided by the least squares method are not meaningful quantities. A linear
estimation of the density fixing the shapes at convergence could provide error
bars for the density contrasts of each anomaly, but these will be specific for
the estimated shape and not representative of the high correlation between
anomaly shape and density jump, so probably still too optimistic. We therefore
appraise the uncertainty associated with each solution by analyzing the spread
of different runs converging to the same interior family (cluster), similarly
to Dinsmore and de Wit (2023). In simulation, where the ground-truth is known,
these uncertainties can be compared to the true errors to assess their
reliability. A set of representative solutions can be generated by starting
from different initial models, as in Section 3.3.7, or by randomly perturbing
the observed measurements within their associated noise, or by a combination
of the two approaches.

Figure 14 shows the average solution and corresponding errors for simulation
runs with the Model 1 ground-truth, assuming no noise and 3 levels of gravity
resolution. Each plot is obtained from statistics over about 200 converged
solutions, each starting from a different random initial model as described in
Section 3.3.7. As shown in that section, we additionally cluster the full set
of solutions based on the values of $I_{zz}$, using 2 cluster centers. Most of
the solutions (95%) are contained in the same cluster, which has high
correlation with the ground-truth model. We then compute the uncertainty
($\sigma_{\rho}$) as the standard deviation of the distribution of density
values of each cell over the set of solutions in a given cluster, while the
average solution is given by the mean of these distributions. We plot the
ratio between the density contrast of the average solution (difference between
the cell density and the median density over the body) and the computed
uncertainties, thus representing the statistical significance of any detected
heterogeneity. We see that inside the target anomaly the density jump is
statistically significant, with the errors generally about 60% of the
deviation. The confidence of the average solutions is lower at the edges of
the target anomaly, where the density contrasts are generally consistent with
0. This is explained by the average solution coming from models where the
retrieved anomaly is smaller and of larger density jump compared to the truth,
or larger and with a lower density jump. Near the center of the true anomaly,
all these solution will have a positive density jump, while around the border
the effect of the denser and the lighter anomalies will cancel out. The last
column in Figure 14 shows the ratio of the absolute values of the true errors
($\Delta\rho=\rho^{(est)}-\rho^{(sim)}$) and the density standard deviation.
We see again that near the center of the anomaly the uncertainties are larger
than the true errors, while at the borders they tend to slightly underestimate
the true errors in the solution.

Figure 14: Average models and associated density uncertainties for multiple
synthetic retrievals based on the Model 1 ground truth for Bennu. Only results
relative to the largest cluster are shown here. The rows indicate the
$l_{\max}$ of the simulated noise-free $\mathcal{O}_{lm}$. The first column
shows the average solutions. The second column shows the significance of the
detected density anomaly, as the ratio of the estimated deviation from the
background density and the computed density standard deviation. The third
column shows the reliability of the estimated uncertainties, as the ratio of
the true errors (average density model minus the ground-truth density) and
$\sigma_{\rho}$. Values larger than 3 in this column are set to 3. The 3D mesh
plots all represent isovalues of the quantity surrounded by slices along the
$X=0$, $Y=0$, and $Z=0$ planes

Similar statistics can be computed for the set of solutions relative to the
scenario of Section 3.3.7 with the ring-core model. Figure 15 shows the
average models and uncertainties obtained from the ring-core ground-truth,
with $l_{max}=10$ and for 3 different measurement noise settings: noise-free
measurements, 100% noise with $\beta=1/3$ (same as Figure 13), and 100% noise
with $\beta=0.03$ (approximating the OSIRIS-REx particle field profile). Each
solution within a set is obtained by starting from a random initial model, but
the values of the observables are the same in each set (same measurement
perturbation), since the aim of Section 3.3.7 was to gauge the influence of
the initial density distribution. For each of these cases, the statistics are
relative to the cluster with highest mean correlation. With perfect
measurements GILA is able to retrieve the under-dense ring and central core
with statistical significance. The ring is also present in the more realistic
$\beta=1/3$ case and with a density contrast on average larger than the
corresponding uncertainty, but the central under-dense region corresponding to
the core in the true model is less well defined. The $\beta=0.03$ average
model presents the inner lighter regions and a small portion of the equatorial
ring, but both anomalies are almost consistent with 0, meaning that a fully
homogeneous distribution could not be ruled out based on these results.

Figure 15: Average models and associated density uncertainties for multiple
synthetic retrievals based on Bennu’s ring-core model of Section 3.3.7. The
rows indicate the average model and its significance, computed as the ratio of
its deviation from its median ($\rho_{0}$) and the standard deviation of the
computed density. The first column was obtained from noise-free measurements,
while for the second and third columns the synthetic measurements were
perturbed following a 100% noise profile with $\beta=1/3$ and $\beta=0.03$,
respectively. The 3D mesh plots all represent isovalues of the quantity
surrounded by slices along the $X=0$, $Y=0$, and $Z=0$ planes. In the second
row, the isovalues are only shown for values larger than 0.5.

###  3.5 Realistic cases

The two hypotheses of piece-wise constant distribution and low number of
distinct anomalies may appear to severely limit the range of applicability of
GILA in real estimation cases. However, we show in the following simulations
that the method could converge to a reasonable approximation of the true model
even when the ground-truth does not respect one of these two assumptions.
Specifically, we treat here the case where the true density distribution is a
smooth function, and the case of a rubble-pile body, where the larger rocks
and voids could be considered anomalies dispersed in a background of finer
material.

####  3.5.1 Continuous target distribution

A smooth variation of the density within a small body could for example be
formed by compression of a porous material following an impact. This is the
characteristic of one of the interior families proposed by Le Maistre et al.
(2019) for Phobos, the "Porous Compressed" (PC) model. Based on the values of
Table 2 in that paper, we construct a member of the PC family by placing right
under the Stickney crater an ellipsoid of semi-major axes ($7.5$, $7.5$,
$3.5$) km in body-fixed Cartesian coordinates. The density contrast inside the
anomaly is then obtained by scaling the signed-distance function from its
border so that the maximum value is 200kg/m^3. Therefore, the density contrast
decreases smoothly away from this peak down to 0kg/m^3 at the border of the
ellipse, where the total density is the same as that outside the anomaly,
namely 1848kg/m^3. In practice, inside the anomaly the density is not exactly
smooth due to the interior grid discretization, for which however the
resolution is relatively high, at 100 cells per dimension. Figure 16(a) shows
the resulting ground-truth model. Hence, from this distribution we produce
synthetic measurements of $l_{\max}=5$ and 100% noise level. The corresponding
GILA distribution estimated over a grid of 50 cells per dimension is shown in
Figure 16(b), where the initial model for this single run of the algorithm is
the same as in Section 3.2. There is indeed a positive density anomaly found
inside the original ellipsoid, although off-centered, and smaller with about
twice the density jump. On the other hand, there are also multiple negative
anomalies on the opposite side of the body, with a density contrast of roughly
-200kg/m^3. To gauge the statistical significance of all these anomalies, we
perform multiple inversions from random models as in Section 3.4. The average
of the 200 solutions is shown in Figure 16(c), where the positive anomaly
below Stickney (+X, -Y direction) is clearly visible, although its density
contrast is up to 10 times the maximum density contrast of the anomaly in the
true model. Figure 16(d) shows the ratio of the density contrasts in the
average model and their spread, computed as the standard deviation
$\sigma_{\rho}$ of the density at each cell across the solutions, as in
Section 3.4. The density jump of the detected anomaly is about 1.4
$\sigma_{\rho}$ away from 0, making it statistically significant.

((a)) Smooth density distribution for the Phobos porous compressed model

((b)) Estimated distribution from a degree-5, 100%-noise gravity

((c)) Average of 200 density solutions from the same degree-5, 100%-noise
gravity and random initial model.

((d)) Density contrast of the average solution with respect to its median
value ($\rho_{0}$), divided by the corresponding standard deviation of the 200
solutions at each cell.

Figure 16: Forward model (a), single-run solution (b), average solution (c)
and its significance (d), for the Phobos porous compressed simulation case.
The 3D mesh plots in the top row share the same colorbar. The 3D mesh plots in
the bottom row all represent isovalues of the quantity surrounded by slices
along the $X=0$, $Y=0$, and $Z=0$ planes. For panel (d), the isovalues are
only shown for values larger than 0.5.

####  3.5.2 Rubble-pile target

Formation and evolution models suggest that most of the asteroids with sizes
between 200m and 10km are rubble-piles (Walsh, 2018). We therefore test here
the ability of the GILA to retrieve an approximate density distribution when
the target body is a rubble-pile. We still assume that the body is composed of
a uniform background density, which could model the finer particles and
corresponding voids while allowing the shape to be perfectly filled. In this
matrix of constant density are dispersed random blocks of larger size ($>1\%$
of the body radius). We model the polydisperse rubble pile interior with a
power-law size-frequency distribution (SFD) for the rubble and for the voids,
as in Tricarico et al. (2021). We select an SFD index of -2 for both the
rubble and the voids distributions, so as to generate models with a larger
fraction of bigger boulders within the body than what is generally observed on
the surface of rubble-piles, where this index is closer to -3 (Tricarico et
al., 2021). Both rocks and voids are modelled as polyhedra, with a random
number of vertices between $8$ and $20$ whose positions are sampled from the
vertices of a spherical mesh of radius given by the SFD. Figure 17(a) shows a
randomly generated rubble-pile target of Eros with both rubble, of density
1500kg/m^3, and voids, dispersed in a matrix of density
$\rho_{0}=\qty{1000}{kg/m^{3}}$. This ground-truth model is used to generate a
set of $\mathcal{O}_{lm}$ coefficients with $l_{\max}=10$. There is no
interior discretization grid in this forward computation, instead the total
heterogeneous $\mathcal{O}_{lm}$ are computed by summing the
$\mathcal{O}_{lm}$ of each polyhedral element, weighted by its excess mass, to
the $\mathcal{O}_{lm}$ of the body mesh weighted by the mass of a homogeneous
body with density $\rho_{0}$.

According to the results in Figure 12, the estimation may be inaccurate for
elements with volume smaller than $1/l_{\max}^{2}$ times the body volume,
based on the half-wavelength resolution of the gravity coefficients. Figure
17(b) shows only the components of the same ground-truth model with relative
size larger than $1/l_{\max}^{2}$, which can be thought of as a low-pass
filter acting on the density distribution. Indeed, only these larger elements
are detected by GILA, as can be seen in Figure 17(c), presenting the density
distribution retrieved from the degree-10 gravity perturbed by 100% noise. As
in Section 3.5.1, the statistical significance of the estimated anomalies is
tested by exploring multiple solutions obtained from the same dataset but
starting from different initial models. Clustering was not performed here, the
distribution of the $I_{zz}$ being close to unimodal. The averaged solution is
plotted in Figure 17(d), where the two negative anomalies and an overdense
region in the central part are still discernible. However, Figure 17(e) shows
that while the negative anomalies have associated errors which are less than
half their density contrast, a large part of the positive anomaly at the
center has error bars covering their full density contrast, making it
consistent with an absence of heterogeneity. Still, GILA appears to be able to
retrieve approximations to complex density models, within the resolution of
the gravity field.

((a)) Random rubble-pile distribution for an Eros-shaped body.

((b)) Ground-truth model when only elements with relative volume larger than
$1/l_{\max}^{2}$ are plotted.

((c)) Estimated distribution from a degree-10, 100%-noise gravity

((d)) Average of 200 density solutions from the same degree-5, 100%-noise
gravity and random initial model.

((e)) Density contrast of the average solution with respect to its median
value ($\rho_{0}$), divided by the corresponding standard deviation of the 200
solutions at each cell.

Figure 17: Forward model (a) and its low-pass-filtered visualization (b),
single-run solution (c), average solution (d) and its significance (e), for
the Eros-shaped rubble pile simulation case. The 3D mesh plots in the top row
share the same colorbar. The 3D mesh plots in the bottom row all represent
isovalues of the quantity surrounded by slices along the $X=0$, $Y=0$, and
$Z=0$ planes. For panel (e), the isovalues are only shown for values larger
than 0.5.

##  4 Preliminary analysis of Bennu real data

We test here the behavior of the method on real measured Stokes coefficient.
Specifically, we use the set of Bennu Stokes coefficient estimated by the
OSIRIS-REx team from optical tracking of ejecta particles in orbit around the
asteroid (Chesley et al., 2020). The particle field solution includes
coefficients up to degree 10, which are also consistent with the lower-
resolution ($l_{\max}=$4$$) spacecraft-radio-tracking estimates (see Figure
3). As discussed in Scheeres et al. (2020), only coefficients up to degree-4
are sensitive to heterogeneities within the body. This means that the
difference between the central values of these low-degree coefficients and a
constant-density field simulated from the available shape model are still
statistically significant given the estimated uncertainties. However, we
expect GILA to handle the non-significance of the higher degrees given the
inclusion of weights in the least squares, which is why the full set of
coefficients is used here in the inversion. On the other hand, although the
published particle field has an associated covariance matrix $\Omega$, we here
set all correlations to 0 and only use the diagonal part, as GILA was
validated on this simpler scenario of independent measurements. Solutions with
the inclusion of the full covariance matrix by Cholesky factorization of
$\Omega$ (Section 2.2.2) show unstable convergence history, possibly
expressing the need for more specific regularization which will be the object
of future work. The effect of this approximation on our solution is hard to
estimate, but however large it may be it leads to an overestimation of the
information content of the data.

Figure 18 shows the density distribution estimated from the Bennu particle
field, using the SPC v42 shape model (Asad et al., 2021), and starting from
the generic initial distribution shown in Figure 5. It mainly features 3
negative anomalies, 1 at the South pole and 2 at mid-latitudes around
$0^{\circ}$ and $250^{\circ}$ longitude in the Northern hemisphere, and a
concentration of positive anomalies around $120^{\circ}$ longitude in the
Southern hemisphere. The heatmap in the bottom panel represents the Bouguer
anomalies (difference between the measured and constant-density gravity
accelerations) over the surface of the body, as shown in Figure 2 of Scheeres
et al. (2020). As discussed in Scheeres et al. (2020), given that only the
heterogeneous component of the Stokes coefficients is used for this plot, the
error due to the evaluation of the spherical harmonic expansion at the surface
(thus within the Brillouin sphere) is estimated to be below 1% everywhere.
Both the measured and the constant-density sets of coefficients are truncated
at degree 4, since above that degree their difference is smaller than the
uncertainties of the measured coefficients. The Bouguer anomalies thus
computed show a strong correlation with the estimated density anomalies.
Figure 18 also shows the convergence history, with a final $\chi_{P}^{2}$ of
0.34, and the RMS of the residuals for each degree, all below the RMS of the
measurements uncertainties. It is worth noting that for this dataset the GILA
model resolution, as defined in Section 3.1.1 (see Figure 4), is always below
the measurement noise (see Figure 3), meaning that the measurements
uncertainties shown here and the corresponding weights appearing in the
computation of $\chi_{P}^{2}$ are strictly the uncertainties of the particle
field.

The inversion was repeated 500 times using multiple random initial models and
perturbing the measurements within their noise, in order to explore the
solution space and check for overfitting. Statistics for the resulting models
are shown in Figure 19, including the distribution of their center-of-mass,
which was unconstrained here since the degree-1 coefficients were not
estimated in the generation of the particle field, but set to 0. As the
distribution of $I_{zz}$ and center-of-mass coordinates are close to unimodal,
no clustering was performed here. The solutions show high variance, but
present the same significant features as the model in Figure 18, namely the 2
negative anomalies in the Northern hemisphere and that at the South Pole, and
the positive anomalies in the Southern hemisphere.

Figure 18: Density estimation from the OSIRIS-REx Bennu particle field
coefficients. RMS of the residuals and the measurement uncertainties for each
degree (a). Evolution of the $\chi_{P}^{2}$ metric over the 500 iterations
(b). Density model at convergence (c). Estimated density distribution in
spherical coordinates, with the associated surface Bouguer anomalies both
mapped onto the surface of the body (see Figure 2 of Scheeres et al. (2020))
and also projected on the $r=0$ plane (d). The density colorbar is shared
between panels c and d. Figure 19: Density solutions from the Bennu particle
field. distributions of reduced $\chi^{2}$, $I_{zz}$, and location of center
of mass (a). Average density distribution (b) and its significance, computed
as the ratio of its deviation from its median ($\rho_{0}$) and the standard
deviation of the computed density (c). The 3D mesh plots all represent
isovalues of the quantity surrounded by slices along the $X=0$, $Y=0$, and
$Z=0$ planes. For panel c, the isovalues are only shown for values larger than
0.5.

While agreeing with the gravity measurements (within their formal errors and
neglecting correlations) and correlating well with the heterogeneous component
of the potential, these density distributions are different from the models
proposed by (Scheeres et al., 2020) and described in Section 3.3.7. Their
ring-core models fit the data equally well, yet GILA was never able to
retrieve such a configuration. In fact, all solutions in Figure 19 are similar
and essentially belonging to the same family of interior, which could again
indicate an incomplete exploration of the solution space. After all, the
possibility of retrieving a ring-core model was already tested in simulations
and proven to be challenging with settings close to that of the real OSIRIS-
REx data ($\beta=0.03$ in Figure 15). On the other hand, a ring-core model for
Bennu is strongly supported by geophysical considerations. The ring
corresponds to the intersection of the Roche lobe with the surface of Bennu,
where loose material is expected to accumulate (Scheeres et al., 2016). The
underdense core could as well be justified by failure modes of rubble-pile
bodies undergoing spin-up (Scheeres et al., 2016; Zhang et al., 2022). The
density distributions retrieved by GILA were instead free of any physical
assumptions. However, the method supports the addition of quadratic
constraints, as explained in Section 2.2.2. We therefore add such constraints
to favor the presence of an equatorial ring of low density, which is likely
for a top-shaped asteroid, by forcing the level-set function of one of the
input anomalies to be close to that of the toroidal anomaly used in the
ground-truth model of Figure 13. As mentioned, this model is sampled from the
set of density distributions proposed by (Scheeres et al., 2020). The density
contrast of this anomaly and the background density $\rho_{0}$ are also
constrained to be close to the density contrast of the ring and the background
density of that same reference model. We select as weights of the level-set
constraints (see Eq. 14) $\nu=$3$$, while for the density constraints the
weight is 20 times larger. Figure 20 shows the corresponding GILA density
solution, starting as usual from the generic initial distribution of Figure 5.
The $\chi_{P}^{2}$ is still close to 1, indicating a good fit of the data, as
shown also by the RMS of the residuals for each gravity degree, with only
degree-2 errors being larger than the measurements uncertainties (by a factor
of about 2.5). The estimated model does present a lighter equatorial ring, as
well as an inner region of low density, which was not enforced explicitly
(although somewhat implicitly by constraining the background density to be
larger than the bulk density). The density contrast of the ring is 35% smaller
in magnitude than that of the ring in the reference model (-400kg/m^3), and
the bulk density 13% smaller than the value of the constraint (1370kg/m^3).
Therefore, the constraints added to the inversion do not appear excessively
tight.

Figure 20: Density estimation from the OSIRIS-REx Bennu particle field
coefficients, with constraints promoting the presence of an underdense ring at
the equator. RMS of the residuals and the measurement uncertainties for each
degree (a). Evolution of the $\chi_{P}^{2}$ metric over the 500 iterations
(b). Density model at convergence in Cartesian coordinates (c), and in
spherical coordinates along with the surface Bouguer anomalies (d).

##  5 Discussion

Overall, simulations have proved that GILA is able to retrieve a target model
in ideal conditions of noise-free data and perfect shape knowledge,
notwithstanding the non-uniqueness of the gravity inversion. Typical target
models where the method fails are those where the degeneracy is higher, such
as the case of concentric shells (not shown here), which would require
adapting the method. Similarly, we don’t expect the method to correctly
retrieve size and density jump of perfectly spherical anomalies, as their
contribution to the gravity field is that of a point mass, although generally
our tests still show a good retrieval of the total mass of the anomaly in such
cases.

We tested the robustness of the inversion approach with respect to various
possible deviations from such ideal settings, primarily a lower resolution of
the observed gravitational potential and presence of noise on the measurements
and the shape of the body. For the more irregular bodies the performance of
GILA was consistently better than for more spherical bodies like Bennu and
Phobos. For these latter two, the range of applicability of the method seems
to be right at the edge of what might actually be expected in many real-
mission scenarios, meaning a gravity field of degree 3 and uncertainty on the
shape of about 0.5% of the body radius. Nonetheless, we stress here that low
values of the correlation metric used in this study may not mean a completely
wrong model compared to the truth, as seen for example in Figure 6. These
limits can therefore be considered conservative, and in any case it is
complicated to give quantitative expectations in terms of predicted
performance in such an ill-posed problem.

Still, a single density solution accounts for very little in the context of
gravity inversion (unless strong constraints are added), since as of yet GILA
has no way to provide realistic uncertainties associated with the estimated
parameters. A possible approach, given the too-optimistic uncertainties output
by the least squares, would be to compute the uncertainties on the converged
solution via Markov-Chain Monte Carlo stochastic sampling about the estimated
parameters, as in Galley et al. (2020). Even then, uncertainties on a single
solution fitting the measurements would not be representative of the other
families of interior equally agreeing with the data within the noise. We
therefore deem it necessary to use the inversion algorithm of GILA multiple
times with different inputs, to then extract statistics from the sets of
solutions thus generated. The proposed approach, including a clustering of the
solutions, and still computationally feasible for a number of solutions in the
order of thousands, would in theory not only provide information on the
uncertainties associated with a single family of solutions, but also an
exploration of the range of possible interior distributions. While this
exploration of the solution space is bound to be limited by computational
constraints, we find that in the cases presented here it is excessively
narrow. In most of our applications of the clustering approach, we detected
less than 3 significantly different families of interior, although generally
the true solution was included in one of those. This was confirmed by visual
inspection, meaning that using higher-order moments of the density function
would have meant separating solutions very similar to each other. This is
certainly a shortcoming for a method ideated for an early exploration phase of
the interior, when additional constraints are not yet available. Hence,
additional efforts are required in removing biases from the random inputs and
from the algorithm itself.

Another important limitation of the method in its current state is the
residuals being limited by the model noise more often than the data noise. In
applications over very precise measurements, this would mean reducing the
information provided by the data and increasing the degeneracy of the problem.
As the model error was a limiting factor in many of the tests shown here, the
conclusions here derived could be seen as conservative. That is, assuming that
improvements on GILA aimed at reducing these convergence issues (such as more
accurate partials, better optimization algorithm, or forward optimization
starting from the converged model) are effectively possible.

It is challenging to quantitatively compare the performance of GILA to any of
the other approaches to the gravity inversion of small bodies mentioned in
Section 1, nor has this kind of comparison with different approaches over the
same test cases been performed as of yet. As mentioned, the piecewise-constant
assumption certainly makes this method more suitable than others for non-
layered bodies which are indeed composed of distinct anomalies, be it because
of fracturing or reaccretion. Nevertheless, as shown in Section 3.5.1, GILA
can still retrieve solution that approximate simple smooth distributions. With
respect to the forward approach to the shape determination of Takahashi and
Scheeres (2014a), the direct estimation of the anomaly shape as part of the
least squares inversion leads perhaps to a more immediate solution, but at the
cost of dealing with the pitfalls of non-convex estimation. On the other hand,
methods such as that of Tricarico (2013) are better suited than GILA if the
density distribution is smooth, as well as providing a more complete view of
the space of possible solutions, while the power of neural networks as
universal function approximators makes the method of Izzo and Gómez (2022)
able to retrieve virtually any kind of interior distribution. In any case,
when dealing with bodies about which little to nothing is known, synergy among
different estimation approaches is essential, as demonstrated in Scheeres et
al. (2020).

As for the application of GILA to OSIRIS-REx data, the flat residuals show
good convergence properties of the method in real applications. The
disagreement between the distribution output by GILA and the more plausible
one proposed by Scheeres et al. (2020) is realistically explained by the lack
of physical constraints in our inversion, along with our apparent incomplete
exploration of the solution space and the inherent non-uniqueness of the
problem. This is confirmed by the convergence to something closer to the ring-
core model as soon as direct constraints are added.

##  6 Conclusions

We have presented GILA, a novel gravity inversion algorithm extending methods
established in Earth geodesy to the gravity inversion of small irregular
celestial bodies. The interior of the body is assumed to be divisible into a
small number of regions where the density is uniform. The shape of each of
these mass anomalies is represented implicitly by a level-set function, which
is adjusted along with the density contrasts and the background density to fit
the measured gravitational potential via iterative least-squares. We have
tested the algorithm over different synthetic scenarios, increasingly more
realistic, and found it to provide reasonable approximations of the true model
in most of the cases. The non-uniqueness of the problem is addressed by
performing several inversions with different initial conditions. Yet one of
the main limitations of the method in its current state is perhaps the limited
exploration of the solution space it provides, which is essential for
asteroids or comets, where most of the time little is known about the
interior. As shown in the processing of Bennu real data, this can be mitigated
by providing additional information to the method. This would imply extending
GILA to the processing of different kinds of measurements, such as local
gravimetric data or observations about the dynamical state of the body. If not
by joint inversion, inputs from other spacecraft instruments or from theory
should be added as constraints. For now, the only type of constraints tested
are those directly on the target shape and density, but extension to general,
non-quadratic constraints of different nature would be needed to increase the
reliability of GILA’s outputs in real scenarios.

##  7 Acknowledgements

This work relies heavily on open-source software, particularly the Python
libraries _NumPy_ and _PyVista_ , the latter employed for mesh manipulation
and some of the plots. Most of the figures presented here were generated using
the _Plotly_ library, along with _matplotlib_ and _seaborn_. We thank the
developers for making their work easily accessible to the community. We also
thank Daniel Scheeres, Matthias Noeker, and Grégoire Henry for their
suggestions and feedback, as well as Véronique Dehant, Özgür Karatekin, Tim
Van Hoolst, and François Massonnet, who supervise the project. This work is
funded by the French community of Belgium, within the frame of a FRIA grant.

##  8 Data Availability

The code for GILA, which was used to generate the data presented in this
paper, is available on GitLab 555https://gitlab-as.oma.be/sbm/gila.

## References

  * Adalsteinsson and Sethian (1995) D. Adalsteinsson and J. A. Sethian.  A Fast Level Set Method for Propagating Interfaces.  _Journal of Computational Physics_ , 118(2):269–277, May 1995.  ISSN 0021-9991.  doi:10.1006/jcph.1995.1098.
  * Asad et al. (2021) M. M. A. Asad, L. C. Philpott, C. L. Johnson, et al.  Validation of Stereophotoclinometric Shape Models of Asteroid (101955) Bennu during the OSIRIS-REx Mission.  _The Planetary Science Journal_ , 2(2):82, April 2021.  ISSN 2632-3338.  doi:10.3847/PSJ/abe4dc.
  * Björck (2015) Å. Björck.  _Numerical Methods in Matrix Computations_ , volume 59 of _Texts in Applied Mathematics_.  Springer International Publishing, Cham, 2015.  ISBN 978-3-319-05088-1 978-3-319-05089-8.  doi:10.1007/978-3-319-05089-8.
  * Blakely (1995) R. J. Blakely.  _Potential Theory in Gravity and Magnetic Applications_.  Cambridge University Press, Cambridge, 1995.  ISBN 978-0-521-41508-8.  doi:10.1017/CBO9780511549816.
  * Carry (2012) B. Carry.  Density of asteroids.  _Planetary and Space Science_ , 73(1):98–118, December 2012.  ISSN 0032-0633.  doi:10.1016/j.pss.2012.03.009.
  * Chao (2005) B. F. Chao.  On inversion for mass distribution from global (time-variable) gravity field.  _Journal of Geodynamics_ , 39(3):223–230, April 2005.  ISSN 0264-3707.  doi:10.1016/j.jog.2004.11.001.
  * Chesley et al. (2020) S. R. Chesley, A. S. French, A. B. Davis, et al.  Trajectory Estimation for Particles Observed in the Vicinity of (101955) Bennu.  _Journal of Geophysical Research: Planets_ , 125(9):e2019JE006363, 2020.  ISSN 2169-9100.  doi:10.1029/2019JE006363.
  * Dinsmore and de Wit (2023) J. T. Dinsmore and J. de Wit.  Constraining the interiors of asteroids through close encounters.  _Monthly Notices of the Royal Astronomical Society_ , 520(3):3459–3475, April 2023.  ISSN 0035-8711.  doi:10.1093/mnras/stac2866.
  * Ermakov et al. (2017) A. I. Ermakov, R. R. Fu, J. C. Castillo-Rogez, et al.  Constraints on Ceres’ Internal Structure and Evolution From Its Shape and Gravity Measured by the Dawn Spacecraft.  _Journal of Geophysical Research: Planets_ , 122(11):2267–2293, 2017.  ISSN 2169-9100.  doi:10.1002/2017JE005302.
  * Ermakov et al. (2014) A. I. Ermakov, M. T. Zuber, D. E. Smith, et al.  Constraints on Vesta’s interior structure using gravity and shape models from the Dawn mission.  _Icarus_ , 240:146–160, September 2014.  ISSN 0019-1035.  doi:10.1016/j.icarus.2014.05.015.
  * Fujiwara et al. (2006) A. Fujiwara, J. Kawaguchi, D. K. Yeomans, et al.  The Rubble-Pile Asteroid Itokawa as Observed by Hayabusa.  _Science_ , 312(5778):1330–1334, June 2006.  doi:10.1126/science.1125841.
  * Galley et al. (2020) C. G. Galley, P. G. Lelièvre, and C. G. Farquharson.  Geophysical inversion for 3D contact surface geometry.  _GEOPHYSICS_ , 85(6):K27–K45, November 2020.  ISSN 0016-8033.  doi:10.1190/geo2019-0614.1.
  * Giraud et al. (2021) J. Giraud, M. Lindsay, and M. Jessell.  Generalization of level-set inversion to an arbitrary number of geologic units in a regularized least-squares framework.  _GEOPHYSICS_ , 86(4):R623–R637, July 2021.  ISSN 0016-8033.  doi:10.1190/geo2020-0263.1.
  * Goossens and Smith (2023) S. Goossens and D. E. Smith.  Gravity degree–depth relationship using point mass spherical harmonics.  _Geophysical Journal International_ , 233(3):1878–1889, June 2023.  ISSN 0956-540X.  doi:10.1093/gji/ggad036.
  * Hedges et al. (2017) L. O. Hedges, H. A. Kim, and R. L. Jack.  Stochastic level-set method for shape optimisation.  _Journal of Computational Physics_ , 348:82–107, November 2017.  ISSN 00219991.  doi:10.1016/j.jcp.2017.07.010.
  * Heiskanen and Moritz (1967) W. Heiskanen and H. Moritz.  _Physical Geodesy_.  Series of Books in Geology. W. H. Freeman, 1967.  ISBN 978-0-608-30923-1.
  * Herique et al. (2022) A. Herique, D. Plettemeier, and W. W. Kofman.  JuRa: The Juventas Radar on Hera to fathom Dimorphos.  2022:NH12C–0296, December 2022.
  * Izquierdo et al. (2023) K. Izquierdo, V. Lekić, and L. G. J. Montési.  An Object-Oriented Bayesian Gravity Inversion Scheme for Inferring Density Anomalies in Planetary Interiors.  _Earth and Space Science_ , 10(7):e2023EA002853, 2023.  ISSN 2333-5084.  doi:10.1029/2023EA002853.
  * Izzo and Gómez (2022) D. Izzo and P. Gómez.  Geodesy of irregular small bodies via neural density fields.  _Communications Engineering_ , 1(1):1–12, December 2022.  ISSN 2731-3395.  doi:10.1038/s44172-022-00050-3.
  * Jamet and Tsoulis (2020) O. Jamet and D. Tsoulis.  A line integral approach for the computation of the potential harmonic coefficients of a constant density polyhedron.  _Journal of Geodesy_ , 94(3):30, 2020.
  * Jekeli (2007) C. Jekeli.  _Potential Theory and Static Gravity Field of the Earth_ , volume 3.  January 2007.  doi:10.1016/B978-044452748-6.00054-7.
  * Kofman et al. (2020) W. Kofman, S. Zine, A. Herique, et al.  The interior of Comet 67P/C–G; revisiting CONSERT results with the exact position of the Philae lander.  _Monthly Notices of the Royal Astronomical Society_ , 497(3):2616–2622, September 2020.  ISSN 0035-8711.  doi:10.1093/mnras/staa2001.
  * Konopliv et al. (2014) A. S. Konopliv, S. W. Asmar, R. S. Park, et al.  The Vesta gravity field, spin pole and rotation period, landmark positions, and ephemeris from the Dawn tracking and optical data.  _Icarus_ , 240:103–117, September 2014.  ISSN 0019-1035.  doi:10.1016/j.icarus.2013.09.005.
  * Konopliv et al. (2018) A. S. Konopliv, R. S. Park, A. T. Vaughan, et al.  The Ceres gravity field, spin pole, rotation period and orbit from the Dawn radiometric tracking and optical data.  _Icarus_ , 299:411–429, January 2018.  ISSN 0019-1035.  doi:10.1016/j.icarus.2017.08.005.
  * Le Maistre et al. (2019) S. Le Maistre, A. Rivoldini, and P. Rosenblatt.  Signature of Phobos’ interior structure in its gravity field and libration.  _Icarus_ , 321:272–290, 2019.  ISSN 0019-1035.  doi:10.1016/j.icarus.2018.11.022.
  * Ledbetter et al. (2021) W. G. Ledbetter, R. Sood, J. Keane, and J. Stuart.  SmallSat swarm gravimetry: Revealing the interior structure of asteroids and comets.  _Astrodynamics_ , 5(3):217–236, September 2021\.  ISSN 2522-0098.  doi:10.1007/s42064-020-0098-1.
  * Lee Rodgers and Nicewander (1988) J. Lee Rodgers and W. A. Nicewander.  Thirteen Ways to Look at the Correlation Coefficient.  _The American Statistician_ , 42(1):59–66, February 1988.  ISSN 0003-1305.  doi:10.1080/00031305.1988.10475524.
  * Levison et al. (2021) H. F. Levison, C. B. Olkin, K. S. Noll, et al.  Lucy Mission to the Trojan Asteroids: Science Goals.  _The Planetary Science Journal_ , 2(5):171, August 2021.  ISSN 2632-3338.  doi:10.3847/PSJ/abf840.
  * Li et al. (2017) W. Li, W. Lu, J. Qian, and Y. Li.  A multiple level set method for three-dimensional inversion of magnetic data.  _GEOPHYSICS_ , 82:1–92, June 2017.  doi:10.1190/geo2016-0530.1.
  * Li and Oldenburg (1998) Y. Li and D. W. Oldenburg.  3-D inversion of gravity data.  _GEOPHYSICS_ , 63(1):109–119, January 1998.  ISSN 0016-8033, 1942-2156.  doi:10.1190/1.1444302.
  * Lowry et al. (2014) S. C. Lowry, P. R. Weissman, S. R. Duddy, et al.  The internal structure of asteroid (25143) Itokawa as revealed by detection of YORP spin-up.  _Astronomy & Astrophysics_, 562:A48, February 2014.  ISSN 0004-6361, 1432-0746.  doi:10.1051/0004-6361/201322602.
  * Matsumoto et al. (2021) K. Matsumoto, N. Hirata, H. Ikeda, et al.  MMX geodesy investigations: Science requirements and observation strategy.  _Earth, Planets and Space_ , 73(1):226, December 2021.  ISSN 1880-5981.  doi:10.1186/s40623-021-01500-6.
  * Michel et al. (2022) P. Michel, M. Küppers, A. C. Bagatin, et al.  The ESA Hera Mission: Detailed Characterization of the DART Impact Outcome and of the Binary Asteroid (65803) Didymos.  _The Planetary Science Journal_ , 3(7):160, July 2022.  ISSN 2632-3338.  doi:10.3847/PSJ/ac6f52.
  * Michel and Fokas (2008) V. Michel and A. S. Fokas.  A unified approach to various techniques for the non-uniqueness of the inverse gravimetric problem and wavelet-based methods.  _Inverse Problems_ , 24(4):045019, July 2008.  ISSN 0266-5611.  doi:10.1088/0266-5611/24/4/045019.
  * Miller et al. (2002) J. Miller, A. Konopliv, P. Antreasian, et al.  Determination of Shape, Gravity, and Rotational State of Asteroid 433 Eros.  _Icarus_ , 155(1):3–17, January 2002.  ISSN 00191035.  doi:10.1006/icar.2001.6753.
  * Osher and Fedkiw (2003) S. Osher and R. Fedkiw.  _Level Set Methods and Dynamic Implicit Surfaces_ , volume 153 of _Applied Mathematical Sciences_.  Springer, New York, NY, 2003.  ISBN 978-1-4684-9251-4 978-0-387-22746-7.  doi:10.1007/b98879.
  * Park et al. (2014) R. S. Park, A. S. Konopliv, S. W. Asmar, et al.  Gravity field expansion in ellipsoidal harmonic and polyhedral internal representations applied to Vesta.  _Icarus_ , 240:118–132, September 2014.  ISSN 0019-1035.  doi:10.1016/j.icarus.2013.12.005.
  * Park et al. (2010) R. S. Park, R. A. Werner, and S. Bhaskaran.  Estimating Small-Body Gravity Field from Shape Model and Navigation Data.  _Journal of Guidance, Control, and Dynamics_ , 33(1):212–221, January 2010.  ISSN 0731-5090.  doi:10.2514/1.41585.
  * Pätzold et al. (2016) M. Pätzold, T. Andert, M. Hahn, et al.  A homogeneous nucleus for comet 67P/Churyumov–Gerasimenko from its gravity field.  _Nature_ , 530(7588):63–65, February 2016.  ISSN 1476-4687.  doi:10.1038/nature16535.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, et al.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Perez-Molina et al. (2022) M. Perez-Molina, A. Campo-Bagatin, and N. Trógolo.  FFT gravity field calculation method and super-ellipsoid generated field.  Technical Report EPSC2022-995, Copernicus Meetings, July 2022.
  * Ritter et al. (2022) B. Ritter, Ö. Karatekin, J. A. Carrasco, et al.  Measuring gravity with the GRASS instrument on the Hera mission.  pages EPSC2022–1115, September 2022.  doi:10.5194/epsc2022-1115.
  * Rücker et al. (2017) C. Rücker, T. Günther, and F. M. Wagner.  pyGIMLi: An open-source library for modelling and inversion in geophysics.  _Computers and Geosciences_ , 109:106–123, 2017.  doi:10.1016/j.cageo.2017.07.011.
  * Scheeres et al. (2000) D. J. Scheeres, B. Khushalani, and R. A. Werner.  Estimating asteroid density distributions from shape and gravity information.  _Planetary and Space Science_ , 48(10):965–971, 2000.  ISSN 0032-0633.  doi:10.1016/S0032-0633(00)00064-7.
  * Scheeres et al. (2015) D. J. Scheeres, D. Britt, B. Carry, and K. A. Holsapple.  Asteroid Interiors and Morphology.  In P. Michel, F. E. DeMeo, and W. F. Bottke, editors, _Asteroids IV_. University of Arizona Press, 2015.  ISBN 978-0-8165-3213-1.  doi:10.2458/azu_uapress_9780816532131-ch038.
  * Scheeres et al. (2016) D. J. Scheeres, S. Hesar, S. Tardivel, et al.  The geophysical environment of Bennu.  _Icarus_ , 276:116–140, September 2016.  ISSN 00191035.  doi:10.1016/j.icarus.2016.04.013.
  * Scheeres et al. (2020) D. J. Scheeres, A. S. French, P. Tricarico, et al.  Heterogeneous mass distribution of the rubble-pile asteroid (101955) Bennu.  _Science Advances_ , 6(41):eabc3350, October 2020\.  ISSN 2375-2548.  doi:10.1126/sciadv.abc3350.
  * Sethian (1996) J. A. Sethian.  A fast marching level set method for monotonically advancing fronts.  _Proceedings of the National Academy of Sciences of the United States of America_ , 93(4):1591–1595, February 1996.  ISSN 0027-8424.
  * Silva and Barbosa (2006) J. B. C. Silva and V. C. F. Barbosa.  Interactive gravity inversion.  _GEOPHYSICS_ , 71(1):J1–J9, January 2006.  ISSN 0016-8033, 1942-2156.  doi:10.1190/1.2168010.
  * Sorsa et al. (2020) L.-I. Sorsa, M. Takala, P. Bambach, et al.  Tomographic inversion of gravity gradient field for a synthetic Itokawa model.  _Icarus_ , 336:113425, 2020.  ISSN 00191035.  doi:10.1016/j.icarus.2019.113425.
  * Takahashi and Scheeres (2014a) Y. Takahashi and D. J. Scheeres.  Morphology driven density distribution estimation for small bodies.  _Icarus_ , 233:179–193, 2014a.  ISSN 0019-1035.  doi:10.1016/j.icarus.2014.02.004.
  * Takahashi and Scheeres (2014b) Y. Takahashi and D. J. Scheeres.  Small body surface gravity fields via spherical harmonic expansions.  _Celestial Mechanics and Dynamical Astronomy_ , 119(2):169–206, June 2014b.  ISSN 1572-9478.  doi:10.1007/s10569-014-9552-9.
  * Tardivel (2016) S. Tardivel.  The Limits of the Mascons Approximation of the Homogeneous Polyhedron.  In _AIAA/AAS Astrodynamics Specialist Conference_ , AIAA SPACE Forum. American Institute of Aeronautics and Astronautics, September 2016.  doi:10.2514/6.2016-5261.
  * Tricarico et al. (2021) P. Tricarico, D. J. Scheeres, A. S. French, et al.  Internal rubble properties of asteroid (101955) Bennu.  _Icarus_ , 370:114665, December 2021.  ISSN 0019-1035.  doi:10.1016/j.icarus.2021.114665.
  * Tricarico (2013) P. Tricarico.  Global gravity inversion of bodies with arbitrary shape.  _Geophysical Journal International_ , 195(1):260–275, 2013.  ISSN 0956-540X.  doi:10.1093/gji/ggt268.
  * Tsoulis et al. (2009) D. Tsoulis, O. Jamet, J. Verdun, and N. Gonindard.  Recursive algorithms for the computation of the potential harmonic coefficients of a constant density polyhedron.  _Journal of Geodesy_ , 83(10):925–942, October 2009.  ISSN 0949-7714, 1432-1394.  doi:10.1007/s00190-009-0310-9.
  * Walsh (2018) K. J. Walsh.  Rubble Pile Asteroids.  _Annual Review of Astronomy and Astrophysics_ , 56(1):593–624, 2018.  doi:10.1146/annurev-astro-081817-052013.
  * Watanabe et al. (2019) S. Watanabe, M. Hirabayashi, N. Hirata, et al.  Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile.  _Science (New York, N.Y.)_ , 364(6437):268–272, April 2019.  ISSN 1095-9203.  doi:10.1126/science.aav8032.
  * Werner (1997) R. A. Werner.  Spherical harmonic coefficients for the potential of a constant-density polyhedron.  _Computers & Geosciences_, 23(10):1071–1077, December 1997.  ISSN 0098-3004.  doi:10.1016/S0098-3004(97)00110-6.
  * Willner et al. (2014) K. Willner, X. Shi, and J. Oberst.  Phobos’ shape and topography models.  _Planetary and Space Science_ , 102:51–59, November 2014\.  ISSN 0032-0633.  doi:10.1016/j.pss.2013.12.006.
  * Zhang et al. (2022) Y. Zhang, P. Michel, O. S. Barnouin, et al.  Inferring interiors and structural history of top-shaped asteroids from external properties of asteroid (101955) Bennu.  _Nature Communications_ , 13(1):4589, August 2022\.  ISSN 2041-1723.  doi:10.1038/s41467-022-32288-y.
  * Zheglova et al. (2017) P. Zheglova, P. Lelievre, and C. Farquharson.  Multiple level-set joint inversion of traveltime and gravity data with application to ore delineation: A synthetic study.  _GEOPHYSICS_ , 83:1–101, October 2017.  doi:10.1190/geo2016-0675.1.
  * Zuber et al. (2022) M. T. Zuber, R. S. Park, L. T. Elkins-Tanton, et al.  The Psyche Gravity Investigation.  _Space Science Reviews_ , 218(8):57, October 2022\.  ISSN 1572-9672.  doi:10.1007/s11214-022-00905-3.