Heterogeneous mass distribution of the rubble-pile asteroid (101955) Bennu.

The gravity field of a small body provides insight into its internal mass distribution. We used two approaches to measure the gravity field of the rubble-pile asteroid (101955) Bennu: (i) tracking and modeling the spacecraft in orbit about the asteroid and (ii) tracking and modeling pebble-sized particles naturally ejected from Bennu's surface into sustained orbits. These approaches yield statistically consistent results up to degree and order 3, with the particle-based field being statistically significant up to degree and order 9. Comparisons with a constant-density shape model show that Bennu has a heterogeneous mass distribution. These deviations can be modeled with lower densities at Bennu's equatorial bulge and center. The lower-density equator is consistent with recent migration and redistribution of material. The lower-density center is consistent with a past period of rapid rotation, either from a previous Yarkovsky-O'Keefe-Radzievskii-Paddack cycle or arising during Bennu's accretion following the disruption of its parent body.

INTRODUCTION

The distribution of a planetary body’s mass is a fundamental physical property that defines its current physical state, which can provide insight into its past, and serves as a geophysical basis for other investigations. Exterior measurements of the gravity fields of planetary bodies have long been used to provide clues and constraints on the interior mass distributions of differentiated bodies. Most recently, NASA’s Dawn mission was able to compare measured gravity fields with the detailed shape and topography of its target bodies, the asteroid (4) Vesta and dwarf planet (1) Ceres, both large enough to have undergone differentiation, and found radial variations in density that provided a wealth of insight into the interior composition and strength of these bodies (–).

For asteroids smaller than ~10 km, which are expected to be rubble piles with no internal differentiation (), there is a more limited history of gravity field measurements and interpretations. The NEAR (Near Earth Asteroid Rendezvous) mission measured the gravity field of the ~10-km asteroid (433) Eros up to degree and order 10 (, ). A direct comparison of Eros’ measured gravity field with the constant-density field computed from its shape was found to be statistically relevant up to degree and order 6 and showed almost no heterogeneity within the body, although there were some indications of a possible surface layer of less dense material (). Furthermore, comparison of its bulk density with its meteorite analog indicated a porosity of 10 to 30% for the asteroid (). From these observations, Eros was classified as a “shattered monolith,” a body whose interior was likely fractured, which had never been fully disassembled or accreted from multiple other bodies. The Hayabusa mission to the ~320-m-diameter asteroid (25143) Itokawa found a very different situation. Although the mission was unable to measure the asteroid’s higher-order gravity field coefficients, the overall bulk density of Itokawa, together with its S-type classification, indicated a porosity of 41% (), which, combined with the visual appearance, demonstrated that this body was a rubble pile. There have been tentative constraints placed on Itokawa’s mass distribution consistent with its “head” having a higher density than its overall “body”; however, these have not been directly tested (, ) and may be explained by other interpretations ().

The OSIRIS-REx (Origins, Spectral Interpretation, Resource Identification, and Security–Regolith Explorer) mission to asteroid (101955) Bennu offers a measurement opportunity not covered by these previous missions. First, the subkilometer size of Bennu (~500 m in diameter) only compares with that of Itokawa, for which no measured gravity field exists. Eros is about an order of magnitude larger and hence several orders of magnitude more massive—at a size where it is resistant to the nongravitational forces and torques that control the evolution of smaller asteroids (). Furthermore, the protoplanets Vesta and Ceres are in a different size class, where self-gravity has shaped their overall structure by pulling these bodies into spheroidal shapes and leading to internal differentiation, owing to the failure of interior material strength. By contrast, the interior of Bennu has pressures on the order of fractions of pascals, too small to actively cause its constituent material to fail, let alone differentiate. The preliminary measurements of Bennu upon the arrival of the OSIRIS-REx spacecraft in late 2018 showed that this body is a rubble pile, given the determined porosity of ~50% based on its total mass, volume, and meteorite analogs (). Thus, the measurement of Bennu’s gravity field opens up a view into the interiors of small, rubble-pile asteroids.

RESULTS

The measured gravity field of Bennu

Owing to Bennu’s small size and low bulk density (~1.2 g cm−3) (), the signal from its gravity field is relatively modest, making it difficult to determine the field to higher spherical harmonic degrees and orders. The measured gravity field of Bennu that we present here is based on two distinct sets of observations and techniques, comprising a previously unavailable strategy for determining the spherical harmonic coefficients of the mass distribution of an asteroidal body. One approach involves direct tracking of the spacecraft using the Deep Space Network (DSN) in combination with spacecraft-based optical imaging. The other approach involves using optical observations of particles that were ejected from the asteroid surface during several mass shedding events observed by the OSIRIS-REx spacecraft (, ).

Gravity field measurement from spacecraft-based tracking

Using traditional radio science techniques, the dynamics of the spacecraft were tracked by a combination of DSN ranging and Doppler measurements with spacecraft-based optical imaging of the asteroid (see Materials and Methods). The DSN tracking was X-band up and down, through both the spacecraft high-gain antenna (HGA) and low-gain antenna (LGA). Tracking passes occurred every day for approximately 5 hours through the HGA and 5 to 15 hours through the LGA. The range measurements constrain the ephemeris of Bennu but do not contribute to the spacecraft trajectory and gravity field estimates, as they do not provide position information relative to the asteroid. The Doppler data are important because they provide measurements proportional to the acceleration of the spacecraft, although it is challenging to observe the gravitational signature of Bennu owing to the many additional accelerations to which the spacecraft is subject, including direct solar radiation pressure (SRP), thermal emissions from the spacecraft, albedo pressure from the asteroid itself (), and occasional propulsive maneuvers by the spacecraft either to change its orbit or, more frequently, to manage its rotational angular momentum (). The optical measurements are crucial, as they can constrain the spacecraft location in the Bennu body frame, achieved by identifying landmarks on the surface of the asteroid and correlating the spin dynamics and shape of the asteroid based on these measurements (). Scanning lidar measurements are also important and were used to validate the estimated models. The highest-quality spacecraft data for these gravity field measurements were obtained during the Orbital B phase of the mission (), when the spacecraft was at its minimum orbit radius about Bennu (average semimajor axis of 925 m and minimum periapsis down to 840 m) (fig. S1). However, earlier spacecraft-based gravity measurements acquired during the Orbital A and Preliminary Survey phases were also used to constrain the values and models ().

Given the altitude of the spacecraft, more than 2.5 Bennu radii in altitude, the gravity field measurements were limited in their signal. Prearrival analysis of the OSIRIS-REx mission plan indicated reliable detection of the gravity field to degree and order 3, and that the degree and order 4 gravity field coefficients would be able to be sensed but not strongly determined (). The fixed geometry of the spacecraft orbit when in close proximity to the asteroid also limited the ability to fully explore the dynamical implications of the mass distribution. Specifically, the low orbits were operated in the vicinity of a “frozen orbit,” which balanced the net effects of gravity and SRP to fix the spacecraft in an orbit with minimal variations (, ). The closest period of orbiting was allowed to oscillate around the frozen orbit, providing periods of increasing and decreasing eccentricity and oscillation of the orbit plane out of the terminator by up to 10°. This orbit orientation is beneficial from a trajectory design point of view, as it enables the qualitative orbit of the spacecraft relative to Bennu to be reliably projected into the future. As a downside, because it actively balances gravity against SRP, it is possible for some of these effects to be aliased into each other. To guard against this, the nongravitational forces were characterized at several different points of the mission, ensuring that the models used in the lowest orbit were accurate (, ).

Estimating and reconstructing the spacecraft’s trajectory required that the asteroid gravity field be estimated, in addition to the nongravitational forces on the spacecraft, spacecraft maneuvers, and other perturbations. In addition, several independent estimates of the gravity field within the team were used to validate the estimates. A difficulty in deriving a precise estimate for the spacecraft trajectory and asteroid gravity models was related to the coupled dependence of estimates of the spacecraft state and the asteroid shape for interpreting the landmark images (see Materials and Methods).

Gravity field measurement from ejected particles

The sporadic particle ejection events observed at Bennu () were fortuitous in terms of gravity field determination (). Many of the ejected particles were placed into temporary orbits about Bennu before reimpacting the surface or escaping from the asteroid. These dynamics were due to the combined interactions between the gravitational attractions and nongravitational forces acting on the particles dominated by direct and reflected SRP (). These particles were observed by the spacecraft navigation cameras and their orbits fit by combining current knowledge of the spacecraft location with observation geometry. In this orbit-fitting process, the gravitational attraction of Bennu proved to be a major model component as a result of the close interactions between the particles and the body. Substantial nonmodeled accelerations were observed; however, as the gravity field is constant over the time spans of interest, its gravity coefficients could be reliably estimated across many of the different particle tracks (, ).

The particles were ultimately much more sensitive to the asteroid’s gravity field than the spacecraft, owing to their very low altitudes above the surface. Thus, although the accelerations acting on these particles could not be as well determined, the precise optical measurements made it possible to distinguish a clear signal in the gravity field coefficients (). The particle-based gravity field estimates did not directly rely on the measured model of the asteroid shape. Thus, they also provided a check on the spacecraft-derived models.

The quoted uncertainties in the particle-derived gravity field include uncertainties of the modeled and unmodeled accelerations that were seen to be acting on these particles, as detailed in (). The limiting uncertainty of the gravity field estimates are at degree and order 9, mainly driven by the lack of detailed shape and photometric properties of the ejected particles and the sparsity of observations (Fig. 1) ().

Fig. 1

Spacecraft- and particle-based gravity fields and uncertainties as a function of degree and order.

Absolute magnitudes (A) and differences between estimated values and the constant-density gravity field (B). Motivated by this comparison, the particle field is used for all subsequent analysis.

Comparison of the spacecraft- and particle-derived gravity fields

The spacecraft- and particle-derived gravity fields are compared in Fig. 1. These plots show the root mean square (RMS) magnitude of the gravity field coefficients for each spherical harmonic degree and the formal uncertainty in these estimates. Direct comparisons (Fig. 1A) show that the spacecraft field only has statistical significance up to degree 4, whereas the particle-derived field has statistical significance up to degree 9. The RMS difference between the particle and spacecraft fields is statistically consistent within measurement errors (Fig. 1B), adding confidence in the particle field estimate.

Also shown in Fig. 1B are differences between the measured gravity fields and a constant-density gravity field computed from the asteroid shape, which indicate the level of statistical significance for interpreting density inhomogeneities in the fields. The uncertainties in the constant-density field due to shape uncertainties are much less than the estimation errors (see Materials and Methods). Here, we see that the spacecraft field only has significance at degree 2, whereas the particle field has statistical significance up to degree 3 with several degree 4 terms still having statistical significance. Table 1 shows the key components from the fields, specifically the gravitational parameter of the particle and spacecraft solutions, their zonals and uncertainties, and the computed zonals from the asteroid shape. Given this comparison, the particle gravity field is used exclusively for the rest of the analysis in this paper.

Table 1

The key gravitational parameters of the universal gravitational constant times mass (GM), normalized zonals, and their uncertainties for the particle and spacecraft field.

Also given are the constant-density (CD) zonals and their comparison to the particle field coefficients. Normalizing radius is 290 m.

	CD shape	Particle field		Particle/CD		Spacecraft field
	Value	Value	1σ	Value	1σ	Value	1σ
GM (m3/s2)	—	4.890450	9.0 × 10−4	—	—	4.894557	2.34 × 10−3
J1	4.7506 × 10−4	—	—	0	0	—	—
J2	1.875791 × 10−2	1.926101 × 10−2	5.20 × 10−5	1.027 × 100	2.78 × 10−3	1.905129 × 10−2	2.38 × 10−4
J3	−1.055026 × 10−3	−1.22194 × 10−3	7.87 × 10−5	1.158 × 100	7.46 × 10−2	−1.308456 × 10−3	4.64 × 10−4
J4	−6.438191 × 10−3	−6.496002 × 10−3	7.10 × 10−5	1.009 × 100	1.10 × 10−2	−6.785849 × 10−3	1.04 × 10−3
J5	−3.601228 × 10−6	6.728866 × 10−5	1.02 × 10−4	−1.868 × 101	2.83 × 101	—	—

Quantifying the density heterogeneity within Bennu

By comparing the difference of the measured (particle-based) and constant-density gravity fields with the particle field uncertainties (Fig. 1B and Table 1), we find that the differences between the gravity fields are larger than the uncertainty in the gravity coefficients only up to degree and order 4, and thus, we focus our density heterogeneity analysis only up to this level. Up to this degree, all the zonals for the measured gravity field are larger in magnitude than the constant-density zonals. This can be ascribed to redistribution of material away from the rotational axis of the body, analyzed in detail later. In addition, the even zonals dominate the gravity field, in keeping with the observed general north-south symmetry of the body [although there are some key differences in the asteroid shape between these hemispheres ()].

Characterizing the impact of the density inhomogeneity

Although we find evidence of inhomogeneity, with the gravity coefficients having a variation of up to several percent, we also wish to know the overall impact of this deviation on the geophysical models of Bennu. These impacts are traditionally indicated using a Bouguer map, which differences the measured gravitational acceleration from the constant-density acceleration on the surface of a sphere above the asteroid (fig. S2). This technique is widely used in planetary geophysics and is well suited to those near-spherical bodies, as the accelerations above the surface mimic the accelerations at the surface. For asteroidal bodies, these maps are not as diagnostic owing to the strong deviation of the asteroid surface from a sphere, meaning that the actual accelerations to which particles are subject on an asteroid surface are not well represented on the surface of a sphere. For example, a change in radius of 5% between two regions can result in a change in surface acceleration of 10% and can represent a markedly different environmental regime, with the same relative variation applicable for a Bouguer map (the asteroid radius varies by ±8% across Bennu’s surface). Thus, it is more useful to instead look at the differences between the measured gravity field and the constant-density reference field mapped onto the surface of the asteroid.

Although mapping the gravity field to the asteroid surface is more useful for irregularly shaped bodies, it is compromised by the well-known convergence issues of a spherical harmonic series expansion when close to or within its circumscribing sphere (). A spherical harmonics gravity field does not converge to the true gravity field at and below the circumscribing sphere about an asteroid. This representation requires higher degrees and order for accuracy as the sphere is approached and diverges from truth (physically) when within the sphere. In contrast, the constant-density gravity field computed from a polyhedron model has no divergence, as it is exact; however, it does not accurately represent the heterogeneous mass distribution. There are numerous approaches to dealing with this issue, but these do not provide a unique model for representing the true surface gravity field (, ). Here, we use the exact potential for the constant-density component of the body and only the spherical harmonics expansions for the non–constant-density contributions (computed as the difference between the measured coefficients and the constant-density coefficients), as proposed in (). This approach, which we call the “split” gravity field, still has some error owing to the physical divergence of the spherical harmonic expansion, yet it can reduce this error as described below. To be statistically significant, we only include terms up to the fourth degree and order.

The relative errors between the constant-density spherical harmonic accelerations and those computed from the exact constant-density field are less than or equal to ±15% (fig. S3, top). The difference between the measured spherical harmonic coefficients and the split gravity field approach is shown to be the same, validating the new representation (fig. S3, bottom). Thus, the error in the surface acceleration due to evaluating the spherical harmonics expansion on the asteroid surface is estimated to be less than 15% overall. The RMS-summed coefficients of the differenced gravity field up to degree and order 4 have a relative magnitude of less than 5%. Multiplying these two contributions indicates that the error due to the spherical harmonic field divergence should be less than 1% everywhere. Thus, using the split gravity field technique, we can generate a modified Bouguer map evaluated directly on the asteroid surface with an estimated error less than 1% (Fig. 2A). We can also make direct computations of slope deviations across the surface between the constant-density and measured gravity field at the same order or accuracy (Fig. 2B).

Fig. 2

Surface variations due to density variations up to degree and order 4.

(A) Difference between measured and constant-density surface accelerations shown on the shape model and in terms of surface latitude and longitude, measured in terms of percent variation (top) and microgals (bottom). (B) Geopotential slope variations using both views.

The acceleration variations due to the nonconstant internal density distributions up to degree and order 4 vary between ±3% over the surface of Bennu (Fig. 2A). A comparison of the surface topography, as defined by the radius from the Bennu center of mass (fig. S4), does not show any strong correlation between gross surface topography and the density variations. Specifically, it does not show mass anomalies in the vicinity of the larger craters along the Bennu equator, which may be expected if these were impact craters formed in the strength regime, resulting in compaction and lithification of the bedrock around the crater (); however, these anomalies may not show up in such a low degree and order field.

The surface slope variations due to nonconstant internal density vary by ±1° across the surface (Fig. 2B). The slope variations are correlated with the surface accelerations, with increasing acceleration yielding a lower slope. This map indicates a few regions where systematic changes in surface slope exist (as compared to the constant-density model) and can be investigated using surface geology techniques, such as developed in ().

Geophysical environment of Bennu given its measured gravity

Using the measured gravity field and our method for mapping the accelerations and potential to the surface, we can compute more definitive surface geophysical maps than were previously possible. The maps shown here account for the fluctuations in internal density, unlike previously published results, which relied on the assumption of constant density ().

The parameters that are most sensitive to the local gravity field are those which should be reevaluated. These are the near-surface dynamical environment, the rotational Roche lobe, and the slope map to determine whether the previously identified transition in slope at the Roche lobe persists with these more accurate data (). Although the constant-density field yields qualitatively consistent results with the current estimated field, specific differences warranting discussion are detailed in the next paragraphs. Other aspects of the surface geophysical environment, such as escape and return speeds and surface acceleration and geopotential trends, are similar enough to the previously reported constant-density shape that they need not be explicitly computed herein.

We first compute all of the co-orbiting equilibrium points in the Bennu-fixed frame (Fig. 3). We verify that the general placement of these equilibria is similar to the constant-density case, although their stability properties are modified. The constant-density models in () and the current constant-density model show that one of the center equilibrium points is stable, meaning that at longitudes in this area, particles could be trapped in low orbit for hundreds of days. However, our computation here using the measured gravity field shows that all of the center manifolds are unstable. This has implications for the dynamics of particles lofted from the surface: Lofted particles are immediately subject to a chaotic dynamical environment across all longitudes. Such a chaotic environment can transport them far from the lofting site, implying a stronger chaoticity for the motion of lofted particles than was expected based on constant-density models. The characteristic exponents of the equilibria provide a measure of the rapidity with which lofted particle trajectories diverge from each other and what their characteristic oscillation times are in the body frame. For the saddle points, the characteristic instability time is on the order of an hour, with their center manifold oscillation periods on the order of 3.5 to 4.5 hours. For the center points, the instability is longer, on the order of 2 to 6 hours, with their oscillation periods on the order of 6 hours.

Fig. 3

The energetics and dynamics associated with close motion about Bennu.

(A) Bennu’s rotational Roche lobe, along with locations and type of its co-orbital equilibrium points. (B) Select trajectories lofted from the surface that fly close to the center manifolds.

To find the rotational Roche lobe () of Bennu (the minimum energy surface that separates the Bennu equator from space, trapping loose material to the surface), we identify the minimum energy value across the different equilibrium points and compute the constant geopotential surface at this value of energy (Fig. 3), which is qualitatively similar to the constant-density computations (). The slopes over the surface of Bennu and the intersection of the Roche lobe with the surface (the “fence”), as computed from the measured gravity field (Fig. 4A), show an even more clearly evident transition in slope than indicated by previous work (). This assessment reinforces the geophysical importance of the rotational Roche lobe for the migration, trapping, and redistribution of surface material. Figure 4B shows the longitudinally averaged slopes and Bennu radius as a function of latitude, quantifying the observed transition seen in Fig. 4A. Here, the average slope uniformly decreases across the Roche lobe transition in latitude.

Fig. 4

Surface slopes on Bennu.

(A) Slopes mapped to the Bennu surface (color scale) showing the intersection of the rotational Roche lobe (the fence; black lines) with the surface. (B) Longitudinally averaged slopes (top) and radius (bottom) as a function of Bennu latitude. The numbers shown in the plots are the average values within and outside of the Roche lobe.

Interpretation of the observed density inhomogeneities in terms of interior structure

Given the measured gravity field, a key question is what constraints can be placed on the actual mass distribution within the asteroid. To address this challenge, we develop hypotheses for mass distributions and then compare these with the measured gravity field using the statistically significant fourth degree and order field (Fig. 1). The higher degree and order terms represent an expansion of degrees of freedom that goes beyond the most meaningful contributions. Because Bennu is a rubble-pile asteroid, the contribution of boulders across a range of size distributions will generate signals at higher degrees and orders, requiring careful consideration to better understand. The goal with the current analysis is to identify the lowest-order systematics in the density distribution.

Considering this relatively low degree and order gravity field, the hypotheses that we consider are relatively simple. We derive a set of hypotheses based on the observed gravity coefficients and theories for mass movement in rubble-pile asteroids. Then, we use two different numerical techniques to test the robustness of the conclusions from the analytical model. The geophysical basis for these theories is laid out in () and is motivated by the past events that could have formed the distinctive shape of Bennu. These are either the surface migration of material to the equator (), a past period of rapid spin rate that caused an internal failure (, ), or the reimpact of a past co-orbital companion onto its surface ().

Analytical assessment

From an analytical point of view, we first view the simplest aspect of the measured gravity field as compared to the constant-density field: the zonal coefficients. Given the global spheroidal shape of Bennu, this analytical assumption is also consistent with the hypothesis of there being a longitudinal symmetry in the mass distribution. Thus, we first look at the zonals to gain clues to the overall mass distribution and then progress to more sophisticated models that can use the full gravity field. Given our constraints, we consider the comparisons between the measured and constant-density values of the zonals J1, J2, J3, and J4.

If an otherwise constant-density body has a central core with a density different than the bulk density, then there is a clear signature in the measured gravity coefficients (). If we define the difference between the measured and constant-density shape zonal coefficients as , then their differences, scaled by the constant-density shape zonal , will all be scaled by the same quantity, −ΔMC/M, at every order l, where ΔMC is the relative mass of the spherical core and M is the total body mass. If the central core has the same overall bulk density as the body, then ΔMC = 0; if the core has a lower density, then ΔMC < 0 and the measured zonal coefficients will be uniformly larger than the constant-density coefficients; and if it has a larger mass, then ΔMC > 0 and the measured zonal coefficients will be uniformly smaller than the constant-density coefficients. All of the measured zonals (except J1, whose case is discussed below) are greater than the constant-density zonals (Table 1), consistent with a possible central mass anomaly in general and particularly with a lower-density core. A lower-density center is one prediction from geophysical models of rubble-pile evolution, as described in (), particularly if the asteroid has had a past rapid rotation. Thus, as one component of our model, we assume a central core, modeled as a sphere centered on the spin axis, but with a possible displacement zC from the equatorial plane. This single anomaly is not able to fully capture the observed density variations, however, requiring additional density components.

The second component of our analytical model is tied to the morphology of Bennu’s shape, with its equatorial bulge that shows up as an average higher radius centered on the equator [Fig. 4B, bottom; (, )]. Also, there is strong evidence for the migration of material toward the equator based on the slope transition at the Roche lobe (described above) and from recent detailed geological analysis of high-resolution images (). A simple geometrical model of the equatorial bulge can be a toroid of central radius RR and inner radius RI, centered at the origin with an overall latitude δR. The gravitational potential of a torus is known in closed form but is computationally difficult (). A simple, yet accurate, approximation is to model it as a simple ring (i.e., let RI → 0 while keeping its mass constant) (, ). This yields a set of formulae for the zonal gravity coefficient of a ring with relative mass ΔMR, defined similarly to the relative mass of the core given above. Independent of the core, one could also expect the equator to be more dense to accommodate the larger values of the measured zonals; however, the two components must be evaluated in tandem to produce a clear picture of the result.

These two components give us four parameters to fit for our four zonal coefficients: ΔMC, ΔMR, zC, and δR (assuming that RR ∼ R0). With this two-component model, the measured zonal coefficient are characterized using the shape zonal coefficient and the four free parameters aswhere l is the degree of the zonal coefficient and R0 is the normalizing radius of the gravity field. The first term accounts for the constant-density contribution, the second term is due to the ring, and the third term is due to the center of mass. Assuming that the offsets of the ring and central core from the equator are small (δR ≪ 1 and zC ≪ R0), the equations decouple into two sets, one involving J2 and J4 from which we can solve for the relative masses and one involving J1 and J3 from which we can solve for the displacements away from the Bennu equator (see Materials and Methods). This occurs as the odd Legendre polynomials go to zero as δR, zC → 0.

Solving these conditions using values from Table 1, we find ΔMR/M = (− 2.35 ± 1.50) × 10−3 and ΔMC/M = (− 5.23 ± 1.69) × 10−2, predicting both an underdense equator and center. The uncertainties are computed from the coupled uncertainties for the two zonal coefficients. The uncertainties are 1σ, implying a weak detection for the underdense equatorial bulge and a stronger detection for the underdense center. We find modest out-of-plane offsets of δR = (0.082 ± 0.052) rad (4.7∘ ± 3∘) and zC = (− 5.88 ± 1.90) m, corresponding to the underdense region of the equator having a displacement into the northern hemisphere and the underdense center having a displacement into the southern hemisphere. Mapping these into geometric regions, we can compute the mass as a function of relative density and the geometric size of a region. For the ring, this corresponds to an inner toroidal radius of 6 to 12 m, depending on whether it is a region of no density or a region with 75% of the bulk density, respectively. The overall ring displacement from the equator is 21 m. Similarly, the central core will have a radius of 108 to 245 m, for a region of no density to one that has 95% of the bulk density, respectively (see Materials and Methods). These regions are shown in fig. S5, where we see that they yield physically plausible distributions that lie within the overall asteroid shape. Given this guidance, we use two numerical approaches that account for the full gravity field to fourth degree and order to more precisely probe such a distribution.

Numerical analysis

We test the robustness of these results by using numerical techniques that can accommodate a general gravity field and fit density variations to match the difference between the measured and constant-density values. We use two distinct techniques to evaluate this: one that provides an exact, but nonunique, density distribution that can account for the difference and one that models the density variations using delineated constant-density regions within the measured Bennu shape and accounting for the gravity field covariance.

The global gravity inversion (GGI) technique (, ) allows one to explore the range of interior mass density distributions compatible with the observed shape and gravity field of arbitrary bodies. It develops the density distribution within the body as a power series in the coordinates, with the coefficients chosen to exactly account for the difference between the measured and constant-density shape gravity coefficients. Because the number of coefficients in the power series expansion is much larger than the number of constraints used in this fit (the 25 coefficients of the fourth degree and order gravity field), there are many unconstrained degrees of freedom in the density fit, with each degree of freedom representing a continuum of possible values—albeit constrained by the data. Thus, there are an infinite number of possible density distributions that can fit the data, although penalty functions can be used to constrain the space of possible solutions (see Materials and Methods).

Applying this approach to the Bennu data up to degree and order 4, while evaluating a penalty function that encourages density variations in the center and equator, yields, as an example, the density distribution patterns in Fig. 5, showing the density variations across a section through the center of mass. These are randomly generated from different seed density distributions, thus exploring a range of possibilities. Some of the solutions match the overall model developed analytically. This does not prove the correctness of the model, but it does show that it can be made consistent with the measured field.

Fig. 5

Randomly generated feasible density distributions computed using the GGI technique.

Images show slices through the y-z plane. Those marked with a yellow star have similar characteristics to the analytical model of an underdense center and equatorial bulge, those with a blue star do not share the geometric characteristics, and those with a piebald star have some of the features. These represent a sampling of possible density distributions that are consistent with the shape and gravity field but are not individually statistically significant.

An alternate method was also applied that delineates regions of the asteroid and evaluates different density values to best fit the measured gravity field (, ). In our implementation of this approach, we account for the covariance of the gravity field estimates to assess the statistical significance of the different density estimates (see Materials and Methods). We use this technique to evaluate the hypothesis that Bennu’s gravity field can be explained by it having an underdense center and an underdense equatorial region. Specifically, we model Bennu’s potential using the three components shown in Fig. 6A, which consist of a spherical mass concentration at the origin to capture the center, an outer layer at radii larger than 245 m and comprising the entire equatorial bulge region, and an inner layer for everything else. For this analysis, we use the full measured gravity field in conjunction with the covariance, and thus, the fits are only sensitive up to degree and order 4, drawing different gravity coefficient values from the uncertainty distribution. We find that the hypothesized model is consistent with the measured gravity field. The inner core has a mass deficit of 6 to 16% of the total mass (more than the 4% deficit found analytically), and the equatorial bulge has a density that is 5 to 12% lower than the bulk density, leading to the middle layer having a density of 8 to 17% more than the bulk density (Fig. 6B). These estimates cannot be directly compared to the analytical estimate, as here a specified volume of the asteroid is marked off as the equatorial region. Also, the gravity field is sampled from its uncertainty estimate; however, the consistency with the simple analytical model is apparent.

Fig. 6

Three-component density model and fits, accounting for the gravity field covariance.

Constant-density components (A) and candidate density and mass distributions (B) to fit the hypothesized variation, computed accounting for gravity field coefficients, and their associated covariance.

From the analytical and numerical assessments, we have multiple lines of evidence consistent with this simple mass distribution hypothesis. This is relevant, as it provides insight into the simplest component of a shape, its center, and the most prominent global feature of Bennu’s shape, its equatorial bulge. The lower densities within the center of Bennu and in its equatorial region both inform our interpretation of the current and past geophysical evolution of Bennu.

DISCUSSION

Interpretations of the hypothesized density distributions

From the density distribution analysis, we can interpret different hypotheses of the past evolution of Bennu and its characteristic shape. Given the underdense regions at the equator and center, we focus on both the migration of surface material and a past rapid spin rate, as both of these are needed to understand the current distribution, as the surface and interior is structurally stable at its current spin period of 4.3 hours (, ). We discount the hypothesis that the equatorial bulge was formed by a past co-orbital satellite that fell back to the surface () because of the lack of supporting evidence.

Characteristics of the equatorial bulge

One of the defining features of Bennu is its so-called equatorial bulge, a region of higher radius (Fig. 4B). Our analysis shows that this region does not exist in isolation, however, and is physically tied to the rotational Roche lobe and a clear break in the surface slope structure. This transition to lower slopes within the Roche lobe was previously hypothesized to relate to recent movement of material as the asteroid spun at an increasing rate (). Furthermore, morphologic evidence of migration of material into the Roche lobe in the recent past has been identified in OSIRIS-REx images, similarly consistent with being driven by the increasing spin rate (). Our identification of lower density in the bulge region is consistent with these other observations and indicates that this migrating material is settling into a more porous state. It is well known that landsliding cohesive powders experience dilation, with up to a 25% increase in porosity [see section 4.6 in ()]. Such a dilation of material would be exacerbated by the geophysical environment on Bennu’s surface, as the surface accelerations are at a minimum in the equatorial region, meaning that the material is migrating into an environment where natural compaction forces are decreasing. Figure 3B shows the acceleration over the surface, as viewed from the north pole, where it can be seen that the magnitude of surface acceleration is decreased by 50% or more as compared to the midlatitudes and polar regions. Thus, as material migrates into the Roche lobe, it also experiences a lower weight and less compaction. In addition, as noted in (), any material within the lobe that is disturbed will be redistributed yet likely trapped within this region—settling back to the surface with very low impact speeds on the order of 1 cm/s or less.

In addition to the identification of lower density and the transport and settling of material into the equatorial region, a higher thermal inertia has also been detected in this area (). However, this goes against typical trends of higher thermal inertia and lower porosity. One possibility is that the asteroid material migrating into the equatorial region is more porous; however, surface materials that are less durable may have been removed by impact phenomenon, leaving a stronger and less porous (and hence higher thermal inertia) covering in this region, as discussed in (). A different possibility is that less-dense micrometer- to millimeter-scale surface particles have been preferentially removed from the surface by electrostatic lofting () due to a lower surface acceleration, thus affecting its surface properties but keeping the subsurface more porous.

Independent of these interpretations, a relevant question is whether this loose material is “primordial” or has been created more recently since Bennu’s transit into the inner solar system. If the observed particle ejection events are tied to Bennu’s current location in the solar system, then the creation of surface material by these mechanisms may have increased more recently. On the basis of the observed particle ejection rate () (and independent of a specific hypothesis), the depth of material over Bennu’s surface was estimated to be accumulating at a rate of 10 cm per 1 million years (Ma) (). Redistribution of this material into the Roche lobe caused by the increasing spin rate owing to the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect () could then explain the currently observed lower slopes and density deficit.

Estimates of a near-Earth residence time of 1.75 Ma () for Bennu would suggest that these formation rates do not produce sufficient material to fully account for the bulge, however, which would require a uniform covering across the entire asteroid of 5 m followed by migration into the Roche lobe to account for its observed size. Also, this leaves the underdense center of Bennu unexplained.

Characteristics of Bennu’s center

Although there is a connection between the current characteristics of Bennu’s equatorial bulge and its geophysical evolution, the tie between its lower-density center and its current state is not as clear. The formation of a lower-density interior has been speculated as one result of a period of fast spin rate (, ), which leads to a failure within the asteroid’s center and creation of an equatorial bulge due to migration of material from inside the asteroid. This is also tied to a finite, but weak, overall cohesive strength of the body. It is estimated that an interior Bennu strength of 1 Pa would require a spin period of 3.5 hours to fail (). However, such a spin period would also lead to loss of surface material and a shaping of the regolith trapped in the equatorial region with a predictable pattern (). The current observations do not clearly show such a signature, leading us to conclude that the underdense center may be a much older geophysical feature of the body, predating the current distribution of surface material toward the equatorial bulge. This indicates that the underdense center formed either at an earlier epoch or when Bennu was originally formed by accretion following the catastrophic disruption of its parent body (, ). This failure mode would also naturally cause a reshaping of the body at the equator, where internal motion and the lower gravitational and rotational acceleration would lead to an oblate shape similar to the averaged profile (Fig. 4B). We evaluate three different possible pathways for such a past fast spin rate to occur.

YORP cycles

The first hypothesis that we consider is that the original failure that led to the underdense center occurred at some point after the formation of Bennu, during a previous YORP cycle. Our current understanding of how YORP cycles work can allow a given body to go through multiple periods of faster and then slower spin (). The maximum spin rate can vary from cycle to cycle but would consist of a period of fast spin followed by a deceleration and a period of slow spin and perhaps tumbling. During these cycles, we would expect motion of surface material toward or away from the equator. Motion of material is toward Bennu’s equator at the current epoch (, ), and, for comparison, the motion of material on (162173) Ryugu is away from its equator at its current longer spin period of 7.6 hours (). If such an oscillation occurs, then it would be expected that there would be regions on the surface that serve as a net trap for material and that they would have older surfaces than other regions that material can move over more easily. These disparities in age, inferred from space weathering expression, have not been found ().

Should Bennu be subject to these periods, it is possible that during one cycle, the spin rate was fast enough to cause an internal failure that created the underdense center and equatorial bulge. However, it apparently was not fast enough to cause a larger-scale global deformation of the body, as has been seen for Ryugu and which caused its longitudinal asymmetry () [although see () for possible longitudinal asymmetries in Bennu]. During one of these periods of more rapid spin rates, it is also possible that material could have been fissioned from Bennu’s equator, leaving a distinctive “divot” along the equator (). An argument against this overall hypothesis is the lack of apparent evidence for the repeated migration of material across the surface of Bennu. Also, given the short YORP time scale of Bennu [1.5 Ma in its current orbit and likely 10 times longer in the main belt ()], the delicate balance between spinning fast enough to cause the interior to fail yet not so fast as to cause a catastrophic failure over multiple cycles would seem to be a statistical anomaly.

YORP equilibrium

Another possibility is that Bennu may have been trapped in a spin-equilibrium between the nominal YORP and tangential YORP effects, as outlined in (), perhaps in the main belt before its migration to the inner solar system at an estimated 1.75 Ma ago (). Once in such an equilibrium, the system could remain there indefinitely until disturbed by a change in its heliocentric orbit (such as injection into the inner solar system) or a perturbation to its spin state by a planetary flyby. This could be consistent with Bennu having YORP cycles early in its history (when its YORP time scale would have been longer, over 10 Ma for its current shape) but subsequently becoming trapped in an equilibrium, stopping the cycling of its spin state, making the trapped regolith less likely to be observed and statistically allowing for fewer cycles and hence a more plausible single fast spin rate that caused its interior failure. This would be consistent with the migration of loose surface regolith more recently, since the transition of Bennu into the inner solar system, and is consistent with its 1.5-Ma YORP time scale.

Early formation

An alternate possibility is an early formation of the inner underdense region and the sculpting of the original equatorial bulge in the aftermath of the catastrophic disruption and reaccretion that created Bennu (, , ). It is well known that the angular momentum of an accreting mass distribution will exert control over how the material settles and the morphology of its final shape. This idea has been shown to apply to stars, explaining the population of binary star systems, Kuiper belt objects (), and rubble-pile asteroids (). The scenario here is that a nominally uniform distribution of material accretes with a given angular momentum, causing the system to spin faster as material accretes on the surface. If the angular momentum is too high, then the system cannot form as a single body and must split into a binary system. Below this critical level of angular momentum, however, the body may still spin fast enough to drive material away from the central region and create an underdense interior and an oblate shape. This implies that the aggregate is in a quasi-equilibrium state and thus that the forces at play will balance between gravitational attraction and centripetal acceleration. The validity of such a model for Bennu requires additional theoretical work to evaluate whether this hypothesized process will produce such an underdense region in the interior of a small, rubble-pile asteroid during its accretion phase.

We conclude that the mass distribution of Bennu, based on analysis of the measured gravity field, has signatures of both currently active processes and of an ancient event. Several lines of evidence show that the equatorial bulge has aspects that are of a recent origin and tied to the migration of material into this region. However, these cannot explain the total size of the bulge nor the underdense central region. Those are consistent with a period of rapid rotation in the past, when Bennu was presumably in the main asteroid belt. These shape features are either from a YORP-induced rapid spin rate or date to Bennu’s accretion in the aftermath of its parent body’s disruption.

MATERIALS AND METHODS

Gravity field RMS representations

It is usual to summarize the values that define a full gravity field by using an RMS computation for every degree. If the gravity field coefficients at a degree n are represented as {C, S}, where m = 0,1, …, n, then the RMS of the gravity field at this degree is . If instead a comparison between two gravity fields is being made, then the RMS of the differenced values {ΔC, ΔS} is taken. Last, if the covariance matrix of the gravity coefficients are denoted as P with the diagonal terms being the variance of each gravity coefficient, or , then the RMS for a given degree is computed as .

Shape model

This study made use of the v42 shape model developed using stereophotoclinometry (SPC) () and slightly higher fidelity than the model presented in () and used in (). Detailed assessment of this model using image-derived limb information and range measurements collected in spring of 2019 by the OSIRIS-REx Laser Altimeter () indicate that it has an RMS uncertainty that ranges from 0.5 to 0.7 m. This RMS error folds in an overall scale error and small-scale surface variations. For instance, images reveal that many intermediate boulders that are about 1 m or so in size are underrepresented in the model.

Assuming an uncertainty of 0.5 m in each vertex, the corresponding uncertainties in the gravity field coefficients were computed for orders 2 to 8, following the method published in (). The uncertainties are uniformly less than the estimation uncertainties in the particle field and are 28% of the estimation uncertainties at order 2, 17% at order 3, and less than 10% at all higher orders. Thus, we do not consider these uncertainties in our analyses.

Gravity field estimation

The Bennu gravity field measurement carried out by the OSIRIS-REx mission involved several teams each using distinct combinations of software tools and data processing techniques. The Radio Science teams were based at the University of Colorado (CU) in the Colorado Center for Astrodynamics Research and at the Jet Propulsion Laboratory (JPL). The navigation and flight dynamics teams were represented by the KinetX Inc. team in Simi Valley, CA and in residence at Lockheed Martin’s Waterton Campus in Denver and by an independent Goddard Space Flight Center team. The gravity estimates and other fitting data from each team were compared against each other and found to converge to the same values within the expected errors. The specific values quoted in the paper are from the JPL spacecraft and particle estimates.

Spacecraft gravity field solution

To model the OSIRIS-REx spacecraft trajectory about Bennu, we used JPL’s Mission Analysis, Operations, and Navigation Toolkit Environment (). For gravity estimate purposes, we considered the Orbital B phase of the mission, particularly the data arc between the M5B maneuver on 28 June 2019 and the M1C maneuver on 6 August 2019. This interval is well suited for gravity estimation because it was maneuver free and conducted in close proximity to Bennu (~1-km orbit radius). The Sun-Earth-Bennu angle was approximately 40° at the start of this period. The spacecraft was tracked at least 8 hours per day throughout this period.

Among the DSN radiometric data (), we processed two-way X/X Doppler using a 60-s compression rate. The Doppler measurements were corrected for ionospheric and tropospheric effects and weighted on the basis of the internal noise level of each pass, with an uncertainty floor set to 2.8 mHz. Ranging data were not integral to the gravity field measurement but were included as they provide information on the Bennu ephemeris (), with minimal constraints on the spacecraft motion about Bennu.

Optical navigation measurements () are derived from a landmark-tracking process using SPC (). Spacecraft imagery from all cameras was used to construct a global shape model of the asteroid, consisting of overlapping networks of landmarks at variety of scales that increased with improving picture resolution (). Optical observations for orbit determination consist of the camera-frame coordinates of the center of all of the landmarks that are visible in a given image. For the purpose of gravity science in the science orbits, all observations come from NavCam 1, part of the TAGCAMS (Touch and Go Camera System) suite (). However, they are registered against a shape model derived from all imaging data. We weighted the landmark location in the images at 5 pixels and assumed a pointing error with an SD of 0.5 pixels per image. The pointing error is applied to all landmarks in the image and does not change their relative geometry as seen from the camera. Additional information on the navigation approach and analysis is given in (, , ).

Owing to the large volume of optical navigation data, to balance the relative weight of Doppler and optical data, we downsampled the set of landmarks to 108, uniformly distributed on the surface of Bennu. This choice was based on prediction tests: When fitting subsets of the data arc, we found that increasing the number of landmarks to more than 108 did not provide meaningful prediction improvements for the data not yet included in the fit. This downsampling also helps mitigate a potential bias caused by optical navigation data, which favor a larger Bennu mass than either a Doppler-only solution or the particle solution, regardless of the orbital mission phase and the specific treatment of the data.

The force model included the 8 × 8 gravity field of Bennu; the gravitation of the Sun, planets, Pluto, and the Moon, including relativistic effects, based on JPL’s planetary ephemeris DE424 (); the gravitational acceleration of 25 main belt perturbers (); and SRP based on a 10-plate model (6 plates for the bus and 4 for the solar arrays) with the spacecraft reconstructed attitude kernels and desaturation burns from telemetry as distributed by JPL’s Navigation and Ancillary Information Facility (). We also included stochastic polynomial accelerations of the order of 3 × 10−12 km/s2 to account for missing terms and mismodeling in the model such as thermal radiation from the Bennu surface or higher-order SRP terms ().

In our orbit determination process, we estimated the OSIRIS-REx initial state, Bennu’s GM and 4 × 4 gravity field [initial values from a constant-density shape model and a priori uncertainties as discussed in ()], impulsive ΔV for each desaturation burn with an a priori uncertainty of 0.5 mm/s, stochastic accelerations in 3-hour batches, areas of the 10 plates in the SRP model (with +Y and −Y solar arrays having the same front area and same rear areas and the +Y and −Y plates of the bus having the same area), landmark locations (a priori uncertainty 50 cm), offset between Bennu center of mass and landmark origin, Bennu pole, prime meridian, and rotation rate. We also had the following consider parameters: Earth orientation parameters, DSN station locations, and media calibrations. All of these data and model uncertainties provided the formal uncertainties shown in Fig. 1.

Particle gravity field solution

We derived a particle-based gravitational field, truncated at degree and order 9, by simultaneously estimating the gravitational field for a number of well-observed particles. This served the dual purposes of constraining the gravitational field of Bennu for follow-on geophysical studies and facilitating more reliable orbit estimation for the many particles that did not have a solid gravitational signal, i.e., those having few detections or short data arcs or both. The detection, tracking, and modeling of the ejected particles are described in (), the underlying dataset for which is available in ().

The final gravity field estimate () was based on tracking particles with more than 30 detections and with observational arcs either more than 6 hours in length or covering more than 80% of the particle lifetime. The latter constraint allowed well-observed suborbital particles to be included in the gravity estimate. This led to a set of 20 particles, along with their relevant particulars such as arc length and number of detections. These particles were ejected from various locations on the surface of Bennu, as detailed in ().

The initial values for the gravity field spherical harmonic coefficients were based on the constant-density model. Uncertainties in these coefficients were based on the modified Kaula rule derived for Bennu (). Degree 1 terms, reflective of offsets between the center of mass and center of volume, were zeroed and not estimated, implying that the estimated gravity field was based at the Bennu center of mass. Leveraging the higher-fidelity modeling from optical imaging, the assumption that Bennu is in simple rotation about its maximum moment of inertia was enforced. Because of this, the C21and S21 gravity harmonics terms were zeroed and not estimated. When estimating the gravity field, the rotation axis orientation was allowed to vary, with constraints from spacecraft radio science. However, the particle tracking data are not strongly sensitive to spin axis orientation, and there was no appreciable deviation from the values.

Bennu’s geophysical environment computations and supporting results

The methods and supporting documentation on how the geophysical environment items of surface acceleration, slope, equilibrium points, and the rotational Roche lobe were computed are presented in ().

Internal density modeling approaches

The different methods used for modeling the internal density distribution of Bennu are given in what follows. First, the supporting analytical analysis for the simple model is given, followed by a more detailed explication of the two numerical approaches.

Analytical assessment of a core and ring

If an otherwise constant-density body has a spherical core with a density different than the bulk density, there is a clear signature in the gravity coefficients. Scheeres et al. () show that the difference between the non–constant-density zonals and the constant-density zonals has a simple formwhere ΔMC is the relative mass of the spherical core, M is the total body mass, is the constant-density zonal as computed from the shape, and ΔJ is the measured non–constant-density gravity zonal minus the constant-density shape zonal. Thus, if the core has a lower mass (meaning ΔMC < 0), then the non–constant-density zonal will be uniformly larger than the constant-density zonal. Given that all of the zonals in Table 1 are greater than constant density, this is consistent with a possible lower-density core.

If the central core is given a displacement along the z axis, along the principal moment of inertia, then its additional contribution to the zonal term is

The next observation is tied to the morphology of Bennu’s shape, with its most distinctive feature being the equatorial bulge. A simple model of this shape can be a toroid around the equator, for which the gravitational potential is known (). However, an accurate approximation to the computationally difficult torus is to model it as a simple ring, which has been shown to be accurate for evaluating the gravitational potential outside of the torus—which is what is needed for this application (). Then, a simpler model is to use the classical result for a narrow ring, given in ().

Applying this model, the contributions of an equatorial ring to the zonal gravity coefficients can be found aswhere ΔJ is the change in the normalized zonal coefficient, RR is the radius of the torus, δR is the latitude of the torus, and P( − ) are the Legendre polynomials. Assuming that δR ≪ 1, , , , and .

Given measurements of the zonals, we can balance them against a simple constant-density shape plus an interior core and equatorial ring, including the offset from the equator signified by zC and the ring latitude δR.

There are four of these equations, for l = 1,2,3,4, and four unknowns, ΔMC, ΔMR, zC, δR. As stated, these equations are nonlinear; however, by assuming that the offset of the ring and core from the equator is small, i.e., δR ≪ 1 and zC/R0 ≪ 1, the equations can be decoupled and reduced to linear equations. Each of the equations can be expressed uniquely as

These equations can be solved in sequence. First, the l = 2,4 equations can be solved for ΔMC and ΔMR, as they are independent of the other parameters. Then, using these solutions, we can solve l = 3 for δR, and last, we can solve for zC from l = 1. In the following, we also take RR ∼ R0, removing that ratio from the equations.

The equations for J2 and J4 can be arranged as

Inserting values for these measured and computed quantities yields a simple linear set of equations, which can be solved to find ΔMR/M ∼ − 2.35 × 10−3 and ΔMC/M ∼ − 5.23 × 10−2, corresponding to an underdense equator and core. To find the uncertainties, we solve the linear equation defined above for the mass deviations as a function of the measured zonal gravity coefficients. Then, we can directly compute the mean and covariance of the mass uncertainties.

With these values, the equation for J3 can be solved for δR = 0.082 rad (4.7∘). Using the mean radius of the equatorial region, 259 m, this corresponds to an offset of 21.24 m above the equatorial plane. Substituting these into the equation for J1 then yields zC/R0 = − 0.0203, which for R0 = 290 m yields an offset zC = − 5.88 m. See fig. S5 for the implied placement geometry of these components. The uncertainties in these solutions are analyzed from the individual equations described above. A first-order analysis shows that and .

It is instructive to interpret these solutions in terms of the geometry and density of the heterogeneities. The volume of a toroidal ring of main radius RR and inner ring radius RI will be , and the relative mass change given a relative density ΔσR = σR − σB, where σR is the density of the toroidal ring region and σB is the bulk density of the asteroid, is . The relative mass of the core can be represented as . Last, the total mass of the asteroid is . Because both of the masses are negative, this implies that σR < σB and σC < σB. Thus, we model ΔσR = − (1 − fR)σB and ΔσC = − (1 − fC)σB. If we take the equatorial radius to equal the mean equatorial radius of 259 m, then as a function of the relative density fraction, we can solve for the radius RI to see whether this makes geophysical sense. If the torus has zero density, then this corresponds to fR → 0, and if it has the same density, then fR → 1. We find that . Thus, given the solved value for ΔMR/M, we find that RI varies over 5.8 → 12 m as fR = 0 → 0.75. Similarly for the core, we find and RC = 108 → 171 m as fC = 0 → 0.75. These relative sizes are illustrated in fig. S5, showing that these solutions yield physically possible distributions.

Numerical analysis: GGI technique

The GGI technique (, ) allows one to efficiently explore the range of interior mass density distributions compatible with the observed shape and gravity field of arbitrary bodies. By exploiting the fact that the gravity of a body depends linearly on its density, we can transform the direct problem into a linear algebra problem. If we expand the density function in power series of the Cartesian coordinates, i.e., , where r0 is a simple normalization length, then the Stokes coefficients {C, S} of the expansion in spherical harmonics of the gravitational field can be calculated via simple volume integrals (, ). The inverse problem becomes then a nontrivial matrix inversion problem, and the density coefficients take the form , where is a nominal solution obtained by direct inversion, the q vectors are an orthonormal basis of the null space, and sq are free factors. This guarantees that the gravity field generated by the body is equal to the one observed and makes the exploration of all the possible interior structure solutions an efficient process. Last, the space of solutions is typically sampled using a simulated annealing Monte Carlo minimization approach () by carefully constructing a cost function that drives the search toward solutions with a given set of desirable features. The cost function is typically obtained empirically by a simple trial-and-error approach. In the specific case of the solutions obtained in Fig. 5, the cost function iswhere the first term is the ratio of the global SD of the density over the average density and the second term is the ratio of the local density ρ over the average density calculated at points j that are only from two specific regions: points with radius r < 0.1 km and points with r > 0.2 km and z < 0.1 km. The first region is simply the central region of Bennu, while the second region is approximately the equatorial ring of Bennu. Because these terms are positive, by minimizing fcost, we find solutions where the density in these regions is lower than average.

Numerical analysis: Polyhedron technique

An alternate method delineates regions of the asteroid and evaluates different density values to best fit the measured gravity field, as explored in general in (, ). The total gravitation potential can be split into N polyhedral components as , where M is the total mass and m is the mass contribution of each component. The U are expanded as a spherical harmonic expansion with a common origin and reference radius, each with coefficients . Thus the total field potential is itself a spherical harmonic expansion, with coefficients equal to the weighted sum of the components’ coefficients . In this analysis, we prescribe a model for Utotal that uses a combination of constant-density polyhedral layers and spherical mass concentrations (mascons) and generate their respective coefficients. We then use Markov chain Monte Carlo (MCMC) to determine the model parameters needed to fit the observed gravity harmonics. We use an implementation of Goodman’s (, ) affine invariant ensemble sampler. Here, we are primarily interested in estimating each m, the mass contribution of each component, so .

The log likelihood function that is used to drive the MCMC is computed using the observed minus computed residuals of the gravity coefficients , where R is the estimated covariance matrix of the observed coefficients as determined via orbit determination. This likelihood function assumes that the observed coefficients follow Gaussian statistics and that this is consistent with the assumptions used in the orbit determination process. We constrain each m such that . This inversion can be formulated as a least squares problem, but we elected to use MCMC as (i) it gave us the flexibility to explore different numbers and combinations of U without having to deal with issues associated with linearization, (ii) it is trivial to add parameters to the estimated state , such as the position of each of the MassCons, and (iii) it can capture non-Gaussian statistics.

We use this technique to evaluate the hypothesis that Bennu’s gravity field can be explained by it having an underdense center and an underdense equatorial bulge. To do this, we model Bennu’s potential using the three components shown in Fig. 6, which consists of a spherical mass concentration at the origin to capture the core, an outer layer comprising the equatorial bulge, and an inner layer for everything else.

We run the MCMC with an ensemble size of 15 chains for 10,000 steps. The samples are initialized in a small region around the constant-density solution. It takes approximately 200 steps to randomize the trials; however, we discarded the first 2000, and then we thinned each chain by only keeping every 10th step, leaving a total of 12,000 independent samples. We find that the central core has a mass deficit of 6 to 16% of the total mass, the equatorial bulge has a density that is 5 to 12% lower than the bulk density, and the middle layer has a density of 8 to 17% more than the bulk density (Fig. 6).