Mass and Shape Determination of (101955) Bennu Using Differenced Data from Multiple OSIRIS-REx Mission Phases

The Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer (OSIRIS-REx) mission collected a sample from the rubble-pile asteroid (101955) Bennu for return to Earth. For the successful Touch And Go sample acquisition maneuver, the shape and mass of the asteroid needed to be known precisely. Here we use a combination of radiometric, image landmark, and laser altimetry data to determine Bennu’s mass, shape, and orientation simultaneously and to verify existing models thereof. Our shape determination consists of estimating a scale factor and three frame rotation angles that apply to both the global digital terrain model (GDTM) and the landmark coordinates. We use a data type called image constraints, where we take the difference of the observation of the same landmark in images taken at two different times. We analyze data from two phases of the OSIRIS-REx mission, Orbital B and Recon B, and show that interphase image constraints greatly reduce interdependencies between estimated parameters for mass, GDTM scale, and biases on the altimetry data. This results in an improved solution for the mass and shape relative to considering a single mission phase. We find Bennu’s gravitational parameter GM to be 4.89256 ± 0.00035 m3 s−2, and we find a scale factor of 1.000896 ± 0.00036 for the altimetry-based GDTM. Using the scaled volume, this results in a bulk density of 1191.57 ± 1.74 kg m−3 , which is within the uncertainties of previous analyses but more precise.


Introduction
The Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer (OSIRIS-REx) spacecraft was launched on 2016 September 8 and arrived at the ∼500 m diameter rubble-pile asteroid (101955) Bennu on 2018 December 3 (Lauretta et al. 2019a), with the primary goal of collecting a sample from Bennu's surface (Lauretta et al. 2017). For almost 2 yr after arrival, OSIRIS-REx collected measurements to characterize Bennu and its microgravity environment using a suite of scientific instruments and spacecraft tracking data (Lauretta et al. 2017. The instruments most relevant to the work presented here include the cameras of the OSIRIS-REx Camera Suite (OCAMS; Rizk et al. 2018;Golish et al. 2020); the cameras of the Touch And Go Camera System (TAGCAMS; Bos et al. 2018), including two fully redundant cameras to support optical navigation and natural feature tracking; and the OSIRIS-REx Laser Altimeter (OLA; Daly et al. 2017).
When OSIRIS-REx arrived at Bennu, it became clear that the asteroid is very different in certain respects from expectations based on pre-encounter knowledge, with surface properties such as albedo variation and roughness being beyond the spacecraft design specifications for sample acquisition (Lauretta et al. 2019a. Global digital terrain models (GDTMs) produced from OCAMS data mapped Bennu's top-like shape (bulging at the equator, more narrow at the poles), confirming that it is a rubble pile formed from reaccumulated fragments of a larger progenitor and has undergone past periods of fast spin ). Analysis of OLA-based GTDMs suggested that a hemispherically asymmetric distribution of boulders established at Bennu's formation was followed by a partial wedge-like disruption . Shape data combined with gravity data were used to characterize Bennu's geophysical environment and parameters such as mass and bulk density; these analyses indicated that Bennu's surface slopes have been changing in its recent past and that Bennu's interior density distribution is heterogeneous, with an underdense center and equatorial bulge (Scheeres et al. , 2020a. Bennu's rotational acceleration, which had been detected from lightcurve analysis , was shown to be most plausibly due to the Yarkovsky-O'Keefe-Radzievskii-Paddack effect (Hergenrother et al. 2019), that is, the absorption of sunlight that is reradiated thermally on an asteroid of irregular shape (Rubincam 2000).
One of the most fascinating finds of the spacecraft encounter was that Bennu is actively ejecting particles of rock from its surface (Lauretta et al. 2019b), something that was hypothesized to be possible (but not observed) prior to arrival at Bennu from polarimetric observations (Cellino et al. 2018). This activity at Bennu proved to be serendipitous for gravity field determination because the observed particles were much closer to Bennu's surface than OSIRIS-REx, increasing the sensitivity to variations in the gravity field (for details of the radio science experiment that was initially conceived for OSIRIS-REx, see McMahon et al. 2018). The particles were detected and tracked in images (Liounis et al. 2020), from which the particle trajectories could be determined (Hergenrother et al. 2020;Leonard et al. 2020). This trajectory analysis was used in gravity field determination that resulted in models with a higher resolution than possible using only spacecraft tracking data Chesley et al. 2020). The effects of the particle ejections on Bennu's rotation and orbit were found to be negligible (Scheeres et al. 2020b).
The selection of a suitable sample collection site was complicated by the unexpected roughness of Bennu's surface . Nonetheless, on 2020 October 20, the spacecraft performed a successful Touch And Go (TAG) sample collection maneuver, and evidence suggests that the collected material is on the order of hundreds of grams, well beyond the mission requirement of 60 g (Lauretta & the OSIRIS-REx TAG Team 2021). OSIRIS-REx will deliver its sample to Earth in 2023 September.
The successful TAG execution required precise knowledge of Bennu's mass and shape, among other geophysical parameters. The spacecraft had to autonomously navigate down to Bennu's surface, requiring uncertainties below 30 cm for the GDTM's overall scale. For this reason, it was necessary to account for any discrepancies between results from different data (such as between the OCAMS and OLA GDTMs) or between parameter estimates from different phases of the mission. For example, analyzing data from different phases of the mission yielded inconsistent estimates of the asteroid's gravitational parameter GM (the product of the gravitational constant G and the body's mass M). The mass loss due to the aforementioned particle ejections is small at ∼ 10 4 g per orbit (Hergenrother et al. 2020) and cannot explain the GM differences. Eventually, it was determined that the problem of inconsistent GM estimates was tied to the difficulty of estimating additional parameters that are highly correlated with GM. This means that it may be difficult to obtain good individual estimates for certain parameters, because their estimate may strongly depend on the value of another parameter with which they are strongly correlated. Their combined effect may be well determined, but we need good estimates for each individual parameter to characterize the asteroid. For example, GM is linked to OSIRIS-REx's orbital period and semimajor axis through a well-known relationship, Kepler's third law, , where P is the orbital period and a e is the semimajor axis. So, a strong determination of GM can be accomplished by a strong determination of both the period and semimajor axis. However, we found that most standard OSIRIS-REx orbit solutions gave a strong determination of the period only and thus could determine only the factor a GM e 3 accurately. In this paper, we describe our research designed to disentangle GM, a e , GDTM scale, and other parameters, with the goal to obtain individual estimates for these parameters for a robust determination of Bennu's mass and shape. The methods described here were developed in collaboration with the OSIRIS-REx navigation and altimetry teams to help in resolving discrepancies that could affect TAG. However, the final results described in this paper were determined independently of the operational time line and did not directly influence TAG operations. Other papers in this Focus Issue will describe the operational response to such discrepancies in preparation for TAG and will assess the final GDTMs produced by the mission.
Perhaps the most important aspect of our approach is the use of a single large solution for orbital and geophysical parameters that spans two mission phases. This also allows us to simultaneously estimate mass and shape parameters and thereby to determine the asteroid's bulk density. Other aspects include (1) the use of altimeter ranges from OLA, with special attention given to off-nadir ranges; and (2) the use of a differential camera data type. This data type has two time tags, which are used to tie different phases together.
This paper is structured as follows. In Section 2, we discuss the data types that we use in our analysis. We present our methodology in Section 3. We show the results of our analysis, including estimates of the geophysical parameters and correlations between them, as well as residual statistics for the data, in Section 4. We evaluate our estimates using independent data, and we show the results of that analysis in Section 5. We discuss our results in Section 6, and we end with conclusions and an outlook in Section 7.

Data Types
In addition to using direct image data acquired by the spacecraft and radiometric data (both Doppler and range, as is standard in planetary geodesy; see Thornton & Border 2000;Moyer 2000) from the Deep Space Network (DSN), we also employ two data types that are used less frequently: (1) a differential data type called image constraints, and (2) altimetric data. We discuss the direct image data and image constraints in Section 2.1 and the altimetric data in Section 2.2.

Landmark Data and Image Constraints
The image data that we use were collected by the navigation camera NavCam 1 of TAGCAMS (Bos et al. 2018) on board OSIRIS-REx. The image data consist of two different types, called landmark (direct) data and image constraints. Direct landmark data use landmarks on the surface and the direction vector from the camera to the landmark to constrain the spacecraft's orbit. Landmarks are fixed points on the body's surface that can be identified in (multiple) images, and they are determined, prior to our use here, through shape modeling algorithms, e.g., stereophotoclinometry (Gaskell et al. 2008) in the case of OSIRIS-REx (e.g., Barnouin et al. 2019Barnouin et al. , 2020. Using a camera model (which describes, among other things, focal length, principal point offset, and distortion), the vector from landmark to spacecraft can be transformed into pixel space, which is the actual measurement. Such a data type has been used often in navigation and in the determination of the gravity field of celestial bodies, especially asteroids (e.g., Konopliv et al. 2002Konopliv et al. , 2014Mazarico et al. 2014b). Landmarks identified in images are directly compared to the prior knowledge of the landmark coordinates, and they thus depend on a priori knowledge of the asteroid shape. Because this measurement depends on the landmark position, it is sensitive to the instantaneous orientation of the central body.
The image constraints measurement (Centinello et al. 2015;Mazarico et al. 2017) relies on the fact that a landmark's location in two images taken at different times should be the same in the proper frame (although not in the images themselves). In this case, landmark pixel coordinates (line and sample) from a direct landmark observation are translated into a pointing vector. This measurement type is illustrated in Figure 1. If at times 1 and 2 we have positions P of the camera on board the spacecraft (a 3 × 1 vector) and unit vectorsv (also 3 × 1) from camera to the landmark, in the body-fixed frame, then the straight-line distance d between the two directions from camera to landmark can be found as follows. For times 1 and 2, one can define a line segment that includes P and has a direction ofv, with a running variable t to describe any point on that line:ˆˆ( The distance d between any two points is d = Pr 2 − r 1 P, and the shortest distance can be found by noting that this occurs for points r 1 and r 2 that are on a line perpendicular to both unit vectors. Hence, we can define a unit vector 2 , and the shortest distance is simply the projection of any distance on this line: Using Equation (1), this becomes Because the unit vectors are perpendicular by definition to the unit normal vectorn, the dependency on t 1 and t 2 vanishes, and d becomes (Centinello et al. 2015) In practice, d will not be zero because of errors on, for example, spacecraft position and attitude, and d can be minimized in a least-squares sense. This measurement type is not dependent on the landmark location or the shape of the body, because all that is required are the directions to the landmarks (Centinello et al. 2015), as can also be seen from Equations (1)-(4), where the exact shortest-distance points (as given by knowledge of the values of t 1 and t 2 ) on the line segments are not needed. Because of this, this measurement type is not sensitive to the instantaneous orientation of the body, but it is sensitive to changes in orientation, as it uses images taken at different times. Also, this measurement is not as geometrically strong as a direct landmark measurement because more combinations of satellite positions together with pointing may allow for a similar distance d discrepancy (Mazarico et al. 2017). However, because it depends on images at two distinct times, it can connect two spacecraft states separated in time. This will prove to be especially advantageous for our analysis, because it will allow a more direct connection between parameters estimated from different mission phases, as we will explain in more detail in Section 3. We will use image constraints from two images taken within one mission phase (intraphase image constraints) and from images taken in different mission phases (interphase image constraints).

Altimetry Data: Off-nadir Data and Altimetry Bias
Laser altimetry data have long been used for geodetic analyses for Earth (e.g., Luthcke et al. 2000Luthcke et al. , 2002Luthcke et al. , 2003Luthcke et al. , 2005 and planetary applications (e.g., Rowlands et al. 1999;Lemoine et al. 2001;Neumann et al. 2001). In particular, altimetric crossovers are based on the assumption that two spacecraft groundtracks that intersect should essentially measure the same topography, in the absence of changes on the surface. Discrepancies can then be used to update the spacecraft position and pointing. These crossovers, as presented in Rowlands et al. (1999), have been used, for example, to improve the orbits of lunar spacecraft (Rowlands et al. 2009;Goossens et al. 2011;Mazarico et al. 2012). In Mazarico et al. (2010), altimetric crossovers were expanded to swath crossovers that use the unique cross-track (the direction perpendicular to the orbital plane) information present in multibeam altimetry data. We refer to them hereafter as altimetry constraints. Altimetry constraints and crossovers also do not depend on prior knowledge of the shape of the body.
Alternatively, the direct altimetry measurement uses the altimetry ranges to derive a topographic profile for a given spacecraft orbit and pointing, and it directly compares these profiles to an existing topography model. Thus, this measurement type relies on previous knowledge of the topography. Recently, with the improved topography data from the Lunar Orbiter Laser Altimeter instrument (Smith et al. , 2017 on board the Lunar Reconnaissance Orbiter (LRO) spacecraft (Chin et al. 2007), this data type has been used for orbit determination of LRO . This data type also greatly improved the orbit determination of the Japanese Selenological and Engineering Explorer mission (Goossens et al. 2020a). It was also used to improve the trajectory of the Japanese Hayabusa2 spacecraft in the vicinity of asteroid (162173) Ryugu . We use direct altimetry in our analysis, as well as altimetry constraints to evaluate our results.
Altimetry data are often taken when the spacecraft is nadir pointing. This poses a problem when simultaneously estimating a scale factor on the GDTM and the body's gravitational parameter GM. This is illustrated in Figure 2. A nadir altimetry measurement cannot distinguish between an error in knowledge of the spacecraft's semimajor axis a e and an error in knowledge of the shape of the asteroid: a higher orbit would still be the same distance away from a larger body as a lower orbit to a smaller body. As stated before, the semimajor axis of an orbit is related to the GM of the body through Kepler's third law, and this means that altimetry cannot readily distinguish between a change in GDTM scale factor or a change in GM. We dealt with this issue in two ways: (1) the analysis of altimetry in two mission phases that have different semimajor axes, and (2) the use of off-nadir altimetry.
The high orbit versus small asteroid problem plays out differently in mission phases that have different semimajor axes. The effect of GDTM scale is the same in different mission phases when the altimetry is nadir pointing. However, the trade-off between GM and semimajor axis is different across mission phases where the semimajor axis is not the same. That is because in Kepler's third law only the semimajor axis is raised to the third power. So, a single solution for geophysical parameters that uses data from two orbital configurations has an advantage.
In addition, even within a single mission phase, the problem plays out differently when off-nadir altimetry is included in the solution. To illustrate this, we can compute the remaining bias in altimetry that results from a different orbit or a different GDTM scale (Figure 2). In Figure 3, we show the altimetry bias versus off-nadir angle for either a 30 cm GDTM decrease or a 30 cm orbit increase. For simplicity in this case, we modeled the asteroid as a sphere and assumed a circular orbit, as this has closed-form solutions, but the principle holds for a more complex shape. We also include a 30 cm instrument bias in the figure (the OLA instrument bias is discussed in Sections 3 and 4.1). When the spacecraft is nadir pointing, the altimetry bias is 30 cm for both. But as the off-nadir angle becomes larger, the effect is different between the two cases. This means that, in principle, data at different off-nadir angles can distinguish between the effects from a GDTM scale factor discrepancy and a difference in orbital height (and, as a consequence, GM). However, the differences in residual altimetry bias are small, which would also mean that the orbit needs to be determined very precisely.

Data Analysis and Parameter Estimation Method
In this section, we first describe the overall orbit determination procedure that is the basis of our estimation of geophysical parameters. We then discuss how the various data types are used in our analysis. We discuss the modeling of both measurements and forces, and then we describe how we estimate parameters such as the asteroid's GM, shape, and orientation, among others. Details about OSIRIS-REx data processing for orbit determination and parameter estimation can also be found in Leonard et al. (2019), and this work shares many of the modeling aspects discussed there.
The precise determination of a spacecraft's orbit in the vicinity of a body is the basis of determining that body's geophysical parameters such as mass. The forces acting on the spacecraft can be inferred from position knowledge, and this results in estimates of geophysical parameters of interest. The standard procedure is to numerically integrate the equations of motion using a start model for the forces acting on the spacecraft, as well as knowledge of the spacecraft's position and velocity at the initial epoch. From the spacecraft's position and modeling of the measurements, computed measurements can be compared with the actual observations. Discrepancies between the two, called residuals, are then used to infer changes to the initial force, measurement, and spacecraft parameters. This requires the simultaneous numerical integration of variational equations that describe the sensitivity of the state with respect to the estimation parameters. From the variational equations, one can build a matrix of partial derivatives of the measurements with respect to the estimation parameters, which can then be solved in a weighted batch leastsquares procedure that minimizes the residuals (e.g., Montenbruck & Gill 2000;Tapley et al. 2004). Our processing software, GEODYN II (Pavlis & Nicholas 2017), uses weighted batch least-squares, and it has been used for the determination of geophysical parameters of Earth and other bodies in the solar system (e.g., Luthcke et al. 2013;Mazarico et al. 2014a;Genova et al. 2016Genova et al. , 2019Goossens et al. 2020b).
The orbit is determined over a continuous span of time called an arc. Arc length is often a balance between sensitivity with respect to certain geophysical parameters, which favors long arcs, and the buildup of errors owing to the mismodeling of forces, which favors short arcs. For low-orbiting spacecraft around the terrestrial planets, we typically use arcs of several days, but for OSIRIS-REx in orbit around Bennu, relatively long arcs are possible. From the various mission phases of OSIRIS-REx we select the Orbital B and Recon B phases for our analysis (see  for a detailed time line of the mission phases flown by OSIRIS-REx). Orbital B was specifically designed for gravity and topography mapping. The spacecraft's distance from Bennu's center during this mission phase was about 1 km, in a closed orbit. During Recon B, OSIRIS-REx orbited at a safe distance from Bennu except when performing targeted observations of candidate sample collection sites. The spacecraft distance to Bennu's center for our Recon B arc varied between 0.9 and 1.5 km. Additional topographic data were collected using OLA during Recon B, including off-nadir data designed to help address the scale and GM issue as discussed in Section 2.2. Our arcs are as follows: for Orbital B, the arc is from 2019 June 28, 00:00:00 Coordinated universal Time (UTC) until 2019 August 28, 00:00:00 UTC, and for Recon B, it is from 2020 January 23, 00:00:00 UTC until 2020 February 8, 16:30:00 UTC. Additionally, we use data from phase Orbital A for evaluation; Figure 2. Schematic illustration of the different components involved in an altimetry measurement. An altimetry measurement taken from nadir pointing cannot distinguish between the effects of a higher orbit and a smaller shape. The bias from off-nadir data, however, can provide different information to decorrelate the two (see also Figure 3). the times for the Orbital A arc are listed in Section 5. The spacecraft distance to Bennu's center during the Orbital A phase was between 1.5 and 2 km.
We use Doppler, sequential range, image (landmarks and image constraints), and altimetry data in our analysis. The data weights and the total number of data points are listed in Table 1. Most of these weights are above their typical noise levels to account for possible systematic signals in the data. Sequential range and Doppler data are X-band data (with an uplink frequency of around 7.9 GHz), and their weights correspond to about 2.8 m for range and 0.5 mm s −1 for Doppler (see Moyer 2000 for an explanation of sequential range, and for conversion factors from range units, RU, and from Doppler data, expressed in Hz, to m and m s −1 ). Our Doppler data have a count interval of 60 s, for which the noise level is 0.1 mm s −1 (McMahon et al. 2018). The range data weight is similar to that used in an earlier analysis, whereas our Doppler weight is half of that of earlier analysis . We focus on using image constraints in our analysis, because, as explained, they are not dependent on the locations of the landmarks and as such not dependent on a priori knowledge of the shape model that we aim to verify. We construct our image constraint data set from a database of landmarks. It would be too computationally intensive to use all possible combinations, so we limit ourselves using several criteria. For intraphase constraints, the two images have to be at least a half day apart so that the viewing geometry has changed enough to prevent image constraints with sight lines that are close to parallel. This is not an issue for interphase landmarks, but there not all landmarks are usable because of resolution differences that can exist between images. For both phases, for common landmarks, we targeted around 50 images per landmark to have a good sample of possible geometries. In our analysis there are typically 30-70 images per landmark. We are especially interested in the performance of the interphase image constraints. All image constraints are weighted at 0.5 m, after evaluation of their general levels of fit (see Section 4.2). Although they showed fits better than 0.5 m, this level was chosen to account for the large number of image constraints.
We also use landmark data in our analysis, and they are weighted at 5 pixels, which is much higher than their typical noise level or level of fit. However, we generate image constraints from the landmark data, so these data types are not independent. To account for this, the landmarks, which typically fit at a few tenths of a pixel, are deweighted. We include the direct landmark data because they help in the determination of the GDTM parameters (discussed below), especially the scale factor.  Figure 2) vs. off-nadir angle, for either a smaller GDTM or a higher orbit. A constant instrument bias is indicated with the dashed horizontal line. The central body is assumed to be a sphere with a radius of 290 m, and the orbit is assumed to be circular, at 1000 m from the center of the body. For altimetry, we use only direct altimetry, and we reserve altimetry constraints for evaluation of our results. OLA has two different transmitters-the Low Energy Laser Transmitter (LELT) and the High Energy Laser Transmitter (HELT) (Daly et al. 2017)-resulting in two different sets of data. During Orbital B, only LELT data were collected, whereas during Recon B, HELT and LELT data were collected. HELT fires at a rate of 100 Hz, and LELT at a rate of 10,000 Hz (Daly et al. 2017). The reported range precision for LELT is 1.1 cm and that for HELT is 2.6 cm, whereas the overall accuracy is reported to be 6 cm for LELT and 31 cm for HELT (Daly et al. 2017). Nonetheless, the altimetry data were weighted at a much higher value of 5 m. This is partly to account for the large amount of altimetry data points due to OLA's high firing rates, so that they do not automatically dominate the solutions. Another reason was to account for possible systematic signals in the altimetry data, which we discuss further below.
OLA is a scanning laser altimeter, and during evaluation of the data it was found that a scale factor exists between azimuth and elevation as measured from the center of the scanning mirror. The ratio of azimuth and elevation factors was found to be 1.0073 . After examining the OLA data, we choose to exclude points 1.15°or more away from the center of the mirror. The remaining center data are thus less affected by the azimuth/elevation scale factor. OLA typically obtains linear or raster scans that span ±5°, and our choice thus limits our spatial coverage. Because OLA scans at a high frequency, the OLA data are down-sampled by a factor of 100.
Our processing software GEODYN includes precise force and measurement modeling. For Bennu's gravity field model, we use a spherical harmonics expansion based on the particle ejection data . This model was determined up to degree and order 6 using spacecraft and particle data, and it was extended to degree and order 15 using shape data and a constant density. As explained below, our approach using spacecraft data (we did not include the particles in this analysis) is not sensitive to the higher degrees and orders. We include third-body perturbations for all planets, the Sun, and Pluto, using the DE 424 ephemerides (Folkner 2011). We use the ephemeris of Bennu based on analysis of OSIRIS-REx spacecraft data (Farnocchia et al. 2019), and we do not adjust it in this analysis. This is discussed further in Section 4.2. Solar radiation pressure on the spacecraft is computed externally using a complex model of the spacecraft and ray-tracing (see also Farrés et al. 2017;Geeraert et al. 2019 for details). The contribution of Bennu's albedo is computed within GEODYN, for which a lower-resolution box-wing model is used, where the spacecraft is modeled as a set of plates with surface properties (Marshall & Luthcke 1994). The spacecraft attitude in inertial space is obtained from the onboard star trackers via spacecraft attitude telemetry. We model the force of thermal reradiation on the spacecraft using expected spacecraft component temperatures trended through telemetry data .
We apply the standard corrections in our modeling of the radiometric data (e.g., Lemoine et al. 2013). We use highaccuracy orientation information for Earth and account for the effects of ocean loading and solid-body tides on the DSN sites. The radiometric data are corrected for the media effects of the troposphere and ionosphere. We use the media corrections provided directly by the DSN and use Niell mapping functions (Niell 1996) to map the zenith delays for the troposphere.
Relativistic effects on the radio data are also computed in GEODYN internally, using the ephemeris information. We use information on the center-of-mass (COM) position and the locations of the antenna phase centers for precision modeling of the path of the radio signal from DSN station to spacecraft. Antenna transponder delay information is also taken into account (e.g., Bertone et al. 2018).
The altimetry data are modeled as round-trip light times from the OLA instrument to the surface of Bennu. For direct altimetry, the bounce point of the measurement can be computed for a given spacecraft position and pointing, accounting for relativistic effects. Together with the altimetry range measurement, this gives the local height of the surface in the Bennu-fixed frame, and this can be compared with topography knowledge from a GDTM (e.g., Mazarico et al. 2018;Goossens et al. 2020a). The discrepancy between the two is the direct altimetry residual, which is minimized by updating the spacecraft position, pointing, and other parameters.
We used the OLA v20 GDTM of Bennu based on altimetry data , in particular, the GDTM version that was preprocessed with the free Generic Mapping Tools (GMT) software (Wessel et al. 2013), as discussed in Daly et al. (2020) and Barnouin et al. (2020). This version of the GDTM does not account for overhangs near the edges of boulders. We selected an OLA-based shape model because, at the time of analysis, this model best represented the actual size and variations of Bennu's shape, because of the geodetic control of the altimetry data. Having the size and shape correct is important for the Surface Feature Navigation algorithm (Gnam et al. 2019) that was used for optical navigation that created the landmark data. Subsequent independent analyses used later GDTMs, including the final OLA GDTM that incorporated a refined orientation, spin rate, and other factors; these analyses are discussed elsewhere in this Focus Issue.
For the implementation in GEODYN, we use maps of the shape of Bennu in a digital elevation model representation (not the triangular mesh version described in Daly et al. 2020), which also does not allow overhangs. The maps have a resolution of 64 pixels per degree and were built from a tiled triangular mesh with a resolution of 20 cm.
For the image data, we create landmark observations using landmarks extracted from the same GDTM. We do not adjust the landmark locations in this analysis. We apply a camera model that has information on the focal length, center pixel location, and distortion (Bos et al. 2020). The estimation parameters are divided into global parameters, which are those that have an effect on all the measurements, and local parameters, which have an effect only on parameters within a certain time period (arc). Our global estimation parameters consist of Bennu's GM and the location of its pole (as expressed by the right ascension, declination, and prime meridian in the inertial frame at J2000; linear changes in each of these; and an acceleration term on the prime meridian). We do not adjust spherical harmonic coefficients of Bennu's gravity field in this analysis, as these were already determined using the combination of spacecraft and particle data, and because the spacecraft is not sensitive to higher degrees and orders (Scheeres et al. 2020a). Local parameters include the state (position and velocity) of OSIRIS-REx at the initial epoch (for both the Orbital B and Recon B phases); scale factors on the ray-traced solar radiation pressure accelerations in spacecraft X, Y, and Z directions (with an a priori uncertainty of 5% to account for any mismodeling); and velocity changes (with an a priori uncertainty of 0.2 mm s −1 ) induced by angular momentum desaturation events, where the reaction wheels undergo an unloading that can impart a residual acceleration on the spacecraft. We also include the estimation of empirical accelerations to account for any mismodeling of other small forces, such as solar radiation, thermal rerediation, and albedo, to the extent that the overall scale factors that are estimated are not sufficient to capture variations in these forces. For the Orbital B phase, we estimate a set of accelerations every 4 hr, consisting of constant accelerations in the along-track, crosstrack, and radial directions, as well as accelerations with a frequency of one orbital period in the along-track and crosstrack directions. These are tightly constrained at 5 × 10 −10 m s −2 , and a time-correlation constraint with a correlation length of 4 hr is applied to ensure a smooth set of accelerations. This approach was selected based on analysis of data from the Gravity Recovery and Interior Laboratory mission, where it was found that the estimation of these accelerations improved the orbit determination and gravity field models (e.g., Lemoine et al. 2013). For the Recon B phase, we estimate three constant accelerations (in the along-track, cross-track, and radial directions) every 12 hr, constrained at 3 × 10 −10 m s −2 . We tested using the same amount of accelerations as we applied for Orbital B but did not find a significant difference, whereas orbit evaluations for Orbital B did show improvements.
We also estimate several geometrical parameters related to the measurements. We estimate four parameters related to the GDTM: a scale factor on the entire GDTM, and three Euler rotation angles on the GDTM to account for errors in the pole orientation model obtained during the GDTM creation process. These are also applied to the landmark coordinates for the landmark data. Of these, the scale factor is a parameter of great interest, as it provides information on the shape of Bennu. We also estimate the difference between the COM around which the spacecraft orbits and the center of figure (COF) to which the GDTM is referenced (called the COF−COM offset). Data tied to the body's surface, such as direct altimetry and direct landmarks, are sensitive to such an offset, and this offset is estimated in three orthogonal directions.
In Section 2.2 and Figure 3, we described a remaining bias in the altimetry due to GDTM-or orbit-scale differences. The instrument itself can also have a bias, which is constant as indicated in the figure. We estimate separate biases for the HELT and LELT altimetry data. Because these are instrument measurement biases in range, the LELT bias is constrained to be the same in the two mission phases. No HELT data were collected in Orbital B, so the HELT bias is estimated from Recon B data only. These are global parameters. We also estimate additional local parameters related to pointing. We estimate pointing corrections in roll and pitch for the OLA data. Here, separate angles for LELT are estimated per mission phase. These pointing corrections are given an a priori uncertainty of 5″. We tested the estimation of pointing corrections from the image constraint data but found that this did not improve the orbits. Image data were often collected with simultaneous attitude information available through star imaging, and we deem the pointing knowledge of the camera to be more than sufficient and thus do not perform further corrections. Finally, we estimate pass-by-pass biases for the sequential range data. The a priori uncertainty for these biases was 10 4 range units (RU), which essentially means they are unconstrained, and typical biases are at the level of 4 RU or less.
We then process the data using GEODYN. Our processing software first adjusts only the local parameters until convergence, which is reached when the weighted root-mean-square (rms) of the residuals changes by less than 0.1%. Then, global parameters are adjusted. After this, we iterate again over the local parameters until the rms at the global iteration change also does not change more than 0.1%. In both local and global iterations, GEODYN forms the partial derivative matrices and data residuals. Together with the prescribed constraints and data weights, the system is solved following the square root information filter approach (Bierman 1977). We apply dynamical data editing in our analysis: we choose to delete data with a residual larger than 3.5 times the rms of the previous iteration, where the value of 3.5 is based on experience, with the goal to keep data with residuals within 3.5 times the (expected) standard deviation of the data. On the first iteration a larger factor is applied to account for large residuals in case the initial values are not determined well enough. This only affects the altimetry data and to a small extent the Doppler data. This will be discussed in Section 4.2.
We stress again that an important feature of our approach is the use of a single large solution for orbital and geophysical parameters that spans two mission phases, within one execution of the analysis software GEODYN. Because the interphase image constraint measurement has two parts in different orbital periods, we simultaneously perform a numerical integration of the equations of motion (and the variational equations) for OSIRIS-REx in the Orbital B and Recon B phase. For bookkeeping purposes, OSIRIS-REx is treated as a separate spacecraft in each phase. This means that we estimate separate local parameters for each phase. In a standard multi-arc analysis, Orbital B data and Recon B data would be analyzed separately and data matrices would be added, after convergence of the local iterations. In our case, the image constraint data connect the two periods because they have sensitivity to local parameters in both Orbital B and Recon B. Thus, for each local iteration where the spacecraft states (and the other local parameters listed above) are adjusted, there is a connection between the two parameter sets. This connection is, of course, also there for the global iterations. In this analysis, both local and global iterations are performed simultaneously on the two spacecraft representing Orbital B and Recon B arcs within GEODYN until convergence.

Results
We focus on the contribution of the interphase image constraint data by comparing two separate analyses: one with and one without these interphase image constraints. In this section, we first present the results for the geophysical parameters of interest. We then discuss post-fit data residuals and results for the off-nadir altimetry that were collected during Recon B. Finally, we discuss the evaluation of our results using independent data.

Geophysical Parameter Estimates and Correlations
We present the results for the geophysical parameters of interest in Table 2. When inspecting the results, the differences between the two cases are small, except for GM, where there is a clear effect including the interphase image constraints. At first glance, these results give the impression that the interphase image constraints do not affect the results much. There are, however, two strong points of evidence that they do strongly affect the results. First, we point to the formal errors (one standard deviation) also quoted in Table 2. For several of the important parameters, especially GM, the formal error for the case with interphase image constraints is much lower. This difference cannot be entirely explained by the difference in the amount of data, because the number of interphase image constraints (103,591 data points), compared to the total amount of data points (see also Table 1), is not enough for a factor of 10 change of the formal error for GM. Some of the pole parameters also have greatly improved formal errors, especially the rotation rate and rotation acceleration. This may be expected because image constraints are especially sensitive to changes in orientation, as explained in Section 2.1. However, several caveats need to be taken into account with the formal errors. These errors are not calibrated, and they reflect the combination of geometrical strength, data weights, and amount of data. We include them here because both cases use the same data weights, and thus at least a relative comparison is worthwhile. As discussed later in this section, however, we do not readily take the formal error results to reflect the errors on our parameters.
The second piece of evidence of the contribution of the interphase image constraints is the strongest and comes from the correlations between the parameters in our estimation procedure. Correlations are computed from the scaled covariance matrix of the inverse problem, where correlations close to ±1 may indicate parameters for which independent, individual estimates are difficult to obtain: their value may be strongly influenced by the value of the parameter with which they are correlated, and although their combined effect may be well determined, the individual values may not be. Low correlations close to zero indicate independent parameters in the estimation process. Although we indicated difficulties with the direct use of the formal errors (and thus the covariance), correlations are insensitive to scaling of the covariance and reflect only the dependencies between the parameters. We are thus more confident in directly interpreting the correlations. The correlations between selected parameters are presented in Table 3. We especially highlight the correlations between GM and key GDTM parameters. Correlations between GM and local parameters (not included in the table), especially force parameters such as the empirical acceleration coefficients and scale factors on the ray-traced accelerations, are very small (0.07 or smaller in their absolute value). The correlations between GM and the GDTM scale factor, as well as those between GM and the altimetry biases, are greatly reduced when  we include the interphase image constraints. The correlation between GDTM scale and GM is greatly reduced from 0.705 to 0.059. In addition, correlations between the altimetry biases and GM are also greatly reduced, from −0.918 to −0.168 for LELT data and from −0.766 to −0.085 for HELT data. This indicates that the interphase image constraints indeed help greatly to break the correlations between parameters such as the GDTM scale factor, altimetry biases, and GM. The results from the case using the interphase image constraints are therefore more robust. Sensitivity of the image constraints with respect to changes in orientation is demonstrated with lower correlations between the prime meridian and the rotation rate and acceleration, meaning that these parameters are now more independent. Especially the inclusion of image constraints with a longer baseline between the two images helps improve the estimation of these parameters. We also note that the formal error on the prime meridian is exceptionally large for the case without interphase image constraints. This is closely related to the high correlations between this parameter and the rotation rate and acceleration. If we fix the rotation rate and acceleration in an estimation without interphase image constraints, the formal error on the prime meridian becomes 0°. 071 (an error of ∼0.05%), much smaller than the value of 5°. 730 (an error of ∼4%) that we find when the rotation rate and acceleration are estimated. This also shows the positive influence of the interphase image constraints on the estimation of Bennu's orientation, as there are no such issues with the prime meridian formal error for that case. Finally, for both cases with and without interphase image constraints, correlations between the rotation rate and rotation acceleration are high and effectively the same: −0.94 without and −0.99 with interphase image constraints.
The correlation between the GDTM scale factor and the altimetry bias for LELT increases and switches sign, whereas that between the HELT bias and GDTM scale factor changes sign but does not change its (absolute) value. This is indicative of the trade-off that remains to some extent between the GDTM scale and the altimetry data. The altimetry data strongly determine scale because they constitute a direct measurement tied to the surface. Landmarks do too, but as pixel measurements they are less sensitive to distance. Among the altimetry data, LELT data dominate because of their amount (see Table 4) compared to HELT data. LELT fires 100 times faster than HELT and acquired hundreds of scans in Orbital B. And although the inclusion of interphase image constraints greatly decouples the estimation of GM from the other parameters, the presence of a bias in altimetry makes it more difficult to entirely decouple this bias from the GDTM scale. Hence, there are remaining correlations between the GDTM scale and altimetry biases, in both cases with and without the interphase image constraints. We also point out that the correlation between the two altimetry biases themselves decreases when interphase image constraints are included.
The HELT bias is about 25 cm larger than the LELT bias. This was suspected during initial analysis of both altimetry types, and the results shown here confirm it. To some extent, this makes the use of the altimetry data less straightforward, as also exhibited by the correlations between the GDTM scale and biases. On the other hand, our analysis showed that there is strength in adding the altimetry data from HELT and LELT to the direct landmark data. Initially, we did not include (direct) landmark data in our analysis, and we found a scale factor slightly less than 1 and a negative LELT bias (with there still being a difference of about 25 cm in the LELT and HELT bias values). However, we found that this scale factor did not fit the landmark data well. By combining both data types (with deweighted landmark data because they are not entirely independent of the image constraints, as indicated above), we were able to find consistent altimetry biases and simultaneously fit the altimetry and landmark data, which otherwise was not possible.
With these results, we find the parameter values from the case using interphase image constraint data to be more robust than the values from the case without the interphase image constraints. Additional analysis, by varying the start value for GM and leaving out altimetry data all together (see more about this in Section 6), indicates variations in parameters that are about 10 times the quoted formal errors in Table 2, with the average values being close to those from the case including interphase image constraint data. Hence, for GM and the GDTM scale factor, our results are This translates to a ∼0.007% error in GM and a ∼0.04% error in the scale. With the value for G = (6.67430 ± 0.00015) × 10 −11 m 3 kg −1 s −2 from Tiesinga  (2020), this GM results in a mass for Bennu of (7.3304 ± 0.0005) × 10 10 kg (with a ∼0.007% error).

Data Residuals
Data fit is a way to evaluate to what extent the selected modeling and parameterization agree with the data. Post-fit analysis can be used to verify whether there are any remaining signals in the data not explained by the models, although, of course, the models were determined with the goal of minimizing the data residuals. Here, we discuss residuals and data fit. In Table 4, we list the number of data points, the rms of fit, and the mean in the residuals, per data type and per mission phase. We found that both cases, with and without the interphase image constraints, produce very similar residuals and statistics, so we present only results here from the case that included the interphase image constraints.
All data fit at levels below their weights as given in Table 1. This is especially the case for landmarks and altimetry because we deweighted those data deliberately. Table 4 makes clear that the solution is dominated by altimetry data. The nature and data rate of OLA provided a large amount of data, which we downselected, as discussed in Section 3. There were only a few periods during Recon B where altimetry data were collected, whereas they were collected over a longer span in Orbital B; hence, there is a discrepancy in the amount of altimetry data per phase. We separate the fits for LELT and HELT, and as Table 4 shows, these fits are close together. LELT data residuals have a higher mean in Recon B, which we discuss below.
Our Recon B arc was much shorter than our Orbital B arc, which is also reflected in the amount of data per phase. Nonetheless, we were able to extract 103,591 image constraints between these two phases. Table 4 shows that there is no discrepancy in fit between intraphase and interphase image constraints, as both sets fit at a level of several centimeters. The remaining mean in all data types is also low, indicating noiselevel or better fits with little to no remaining systematic offsets. We estimated a bias for radiometric range data, so these data are expected to have small means. The estimated range biases are shown in Figure 4. They are at the level of 4 RU or less, indicating low levels of unmodeled deviations in Bennu's ephemeris and remaining media effects on the range data. The biases show some remaining ephemeris errors for the Orbital B phase. For the Recon B phase this is also the case, but the time span may be too short to indicate ephemeris errors clearly.
In Figure 5, we show the residuals for each data type, per mission phase. Doppler and range data were routinely collected throughout the mission. The different number of data points in Table 4 thus reflects only the difference in arc length, not any difference in the collection of Doppler data during the mission phase (and the same is true for the other data types apart from altimetry, which we discuss below). Doppler data for Orbital B show very little remaining signal, while Doppler for Recon B shows a slightly higher rms and mean, and some signal in several passes. However, increasing empirical accelerations for Recon B, or even estimating additional biases on the Doppler data in Recon B, did not improve the orbit accuracy, as found from orbit evaluations and comparisons (such as using altimetry constraints, which we address in Section 5, or by computing orbit differences between cases). We thus opted not to increase the number of accelerations, and not to include biases, for Recon B. Range data for both phases show some pass-by-pass remaining signal. As stated above, remaining signal in the range data likely comes from Bennu ephemeris errors (see also Figure 4). Although GEODYN is capable of simultaneously estimating central body ephemeris and other parameters (e.g., Genova et al. 2018), this is left for future analysis. The Bennu ephemeris was updated during the mission, and ephemeris errors affect only the radiometric range and (to a lesser extent) Doppler data, not the body-fixed camera and altimetry data.
The landmark data fit well below the pixel level, despite deweighting them at 5 pixels. They are also well distributed throughout the mission phases and contribute strongly to the overall orbital accuracy. The same is true for the image constraints, which fit very well at a level of several centimeters. As stated before, landmark and image constraint data are not independent, as image constraints are built from landmarks seen in separate images. The levels of fit for these data indicate the quality of the image data and the fidelity of the shape model used (in the case of the landmark data). Because images were routinely collected for optical navigation, we were able to select plenty of image constraints in each phase. Initial tests and orbit evaluations showed that an even distribution with many image constraints provided us with the best orbits. Image constraint time stamps occur often at the same time, as one image can contain many landmarks for which constraint pairs can be constructed.
The altimetry residuals in Figures 5(i) and (j) are clipped because of the dynamical editing that we applied, as discussed in Section 3. This also affected the Doppler residuals and range (none of the other data types had points edited), but to a much lesser extent, and it is thus not as visible in Figures 5(a) and (b). For the altimetry the effect is stronger, but out of about 1 million altimetry data points, only 7600 were edited. When all residuals are plotted as a histogram (not shown), they lie on a very narrow distribution. This is, of course, also reflected in the rms and mean in Table 4. When inspecting residuals with an absolute value larger than 1 m (which includes those that were edited in the processing), we find some instances of limited clustering around particularly large boulders. Daly et al. (2020) indicate that the GDTM that we use here models the overhangs of boulders poorly, and larger residuals there can thus be expected. Figure 5(j) shows both LELT and HELT residuals for Recon B, and both fit at a similar level and show no remaining systematic signal. Figure 5 also shows that the distribution of altimetry in Recon B was much sparser (in the temporal sense, but as a consequence also in the spatial sense), due to operational limitations on the collection of altimetry data. As stated before, our Recon B arc length was chosen such as to include the available altimetry. We discuss altimetry further in Section 4.3.
In Figure 6 we show residuals for the interphase image constraints. They fit at a level very similar to the intraphase data, as can also be seen from the rms and mean in Table 4. Instead of using time stamps, we now plot the residuals versus Figure 5. Data residuals for each data type and per mission phase. The left column has residuals for Orbital B, and the right column has residuals for Recon B. We show (a, b) Doppler, (c, d) range, (e, f) landmark, (g, h) intraphase image constraints, and (i, j) direct altimetry residuals. For each type and phase, the mean is indicated with a dashed-dotted line, and one standard deviation with solid black lines. We show interphase image constraints separately in Figure 6. their time separation. For the interphase constraints, the first image was always taken during the Orbital B phase. These image constraint residuals also show no remaining systematic signal, and there is no trend of residual versus time separation. This is also the case when we plot the intraphase image constraints in this way (not shown).

Off-nadir Altimetry Results
Following the discussion in Section 2.2, we present altimetry residual statistics versus the off-nadir angle. In Figure 7 we show the mean, rms, and number of data points for altimetry residuals versus off-nadir angle. We grouped the residuals per one degree bins for clarity. Figure 7(c) shows the number of data points in each bin. For bins with relatively few points, the results are less conclusive. Figure 7(a) shows the mean in the off-nadir bins, which is the bias that remains because of a GDTM-or orbit-scale issue that we discussed in Section 2.2 (different from the constant HELT and LELT biases that we estimated). Figure 7 shows that only a small signal remains in the Orbital B data. The overall mean of altimetry residuals is very small, at −2.2 mm. There is some signal still visible in the mean for Orbital B, which, following Figure 3, could be a remaining effect of a further GDTM-scale increase or orbit decrease. However, the signal is very small at a few centimeters at most for an off-nadir angle of 4°, but below 1 cm for angles smaller than 2°. For Recon B, the remaining signal is stronger. Interestingly, at nadir, where the full effect of a GDTM scale or orbit scale should be apparent (unless they cancel or work against each other), the mean is only −8 mm. At larger off-nadir angles, the remaining bias is larger, but, despite efforts to point off-nadir, data are few, which may have resulted in Recon B not contributing as much as Orbital B. And overall, the bias is still 10 cm at most. From the results in Table 4, the remaining altimetry bias in Recon B comes mostly from the LELT data. The combination of only a faint signal in the dominating Orbital B data set and the small remaining bias at nadir for Recon B leads us to conclude that if there are any remaining GDTM or orbit effects, they are small, and much below the desired GDTM-scale accuracy level of 30 cm. The rms level of the altimetry residuals per off-nadir bin as shown in Figure 7(b) is flat, indicating that there are no systematic effects remaining in the rms per off-nadir angle.

Evaluation with Independent Data
We also evaluated our results with independent data. Especially because the data fits for the cases with and without the interphase image constraints are not conclusive, an independent evaluation can help test whether the estimated parameters also fit data from a different type or from yet a different orbital period. Here, we test our results in two ways: (1) with a data type not used in our main analysis, called altimetry constraints, from Orbital B altimetry data; and (2) with data from a period not used in our main analysis, Orbital A. Our Orbital A arc was from 2019 January 1, 00:00:00 UTC until 2019 February 28, 00:00:00 UTC.
Altimetry constraints express the discrepancy in total position between two different altimetry tracks and can be used as a measure of orbit accuracy. They are explained in detail in Mazarico et al. (2010). Altimetry constraints are independent of the GDTM because of their differential nature,  like image constraints. We constructed a set of altimetry constraints from OLA data collected during Orbital B. While Orbital B OLA data were used in our main analysis, they were used as direct altimetry. We evaluated our results as follows: after our global iterations to determine the geophysical parameters, we took the resulting orbit for the Orbital B period and evaluated altimetry constraints along this orbit. There was thus no additional estimation of parameters involved when processing these data.
In Table 5 we show the rms and mean of the altimetry constraint residuals for several cases: when using the initial parameter values (such as a GDTM scale factor of 1 and no COF−COM offset), or when using the new values. For both of these cases, we looked at the final orbits from runs with and without the interphase image constraints. Although the values are close, the fit and mean improve when the interphase constraints are used. Moreover, the best results are obtained with the new values and with the interphase constraints. The results deteriorate slightly with the new values for the case without the interphase constraints. This strengthens our conclusion that the results from the case with interphase constraints are more robust, although it should be noted that the results for altimetry constraints are very close among the four cases. Overall, they give an approximate impression of orbital accuracy for OSIRIS-REx at Bennu, which for Orbital B is several tens of centimeters.
We also processed data from the Orbital A phase with our newly estimated parameters. In this analysis, we fixed the global parameters but determined the orbit using radiometric data, image constraints, altimetry, and landmark data. Here, we focus especially on the landmark residuals, which are shown in Figure 8, because they appear to be most affected. The fits for the other data types are always indistinguishable for all runs, i.e., regardless of whether we use the new or old values. We again estimated empirical accelerations as in the Orbital B data, but instead of ray-traced solar radiation pressure, the box-wing model was used here. Data weights were the same as listed in Table 1. Using a value for the GDTM scale factor of 1 and no GDTM rotation angles, etc., we notice strong signals in the landmark residuals. These disappear when we process Orbital A data with our new values. Because these data are independent of our solution, this indicates improvements in the values of the estimated parameters: the new parameters fit the Orbital A landmarks much better, with no systematic periodic signal remaining, whereas there is considerable systematic signal when we process the data without these changes to the parameters. We performed an additional analysis of Orbital A data where we processed the data and updated only one set of global parameters (for example, only the scale factor, or the scale factor and GDTM rotation angles, or only the orientation parameters). This allows us to determine the contribution of each parameter to improvements in the landmark fit. This analysis indicates that the landmark fit improves just slightly when we update only the scale factor. When we update the scale factor and GDTM rotation angles, the landmark fit is close to that of the case with the entire set of new global parameter values. The landmark fit also improves when we only update the COF−COM set of parameter values.

Discussion
The inclusion of interphase image constraints greatly reduces correlations between parameters. Yet some correlations remain, mostly between the GDTM scale and altimetry biases, and to some extent between the HELT altimetry bias and GM, even though the latter correlation is relatively small at −0.168. It also appears that the value for GM does not change much from its start value, which may be counterintuitive given the relation between orbit, GDTM scale, and GM, as outlined in Sections 1 and 2.2. To some extent, despite using data from different orbital periods, the altimetry bias and GM still trade: a lower GM value, for example, 4.8915 m 3 s −2 , results in larger altimetry biases, whereas the difference between the LELT and HELT bias remains about 25 cm. Yet for such a case we also find a remaining 2 cm bias for Recon B altimetry data at the nadir direction, whereas it was closer to zero with our reported GM estimate (see Section 4.3).
As stated in Section 3, altimetry data are deweighted in our analysis by using a weight of 5 m. We also investigated using a weight of 50 cm for the altimetry data, which is closer to the actual level of fit we obtain. This obviously improves the altimetry fit and also results in an even flatter bias for off-nadir results. However, it also affects the fit of the direct landmark data, which now fit at a level more than twice as high as previously and, more importantly, have a substantial mean in the residuals. This indicates a remaining discrepancy between the landmark and OLA data, but as we discussed in Section 4.1, the combination of both with our relative weights allows us to fit both types simultaneously.
As stated in Section 4.1, we also investigated cases without altimetry data, to probe the sensitivity of the solutions with respect to this data type. For such a solution, we find that the scale factor value slightly increases to 1.0009193. Leaving out altimetry data also does not change the start GM value much. We performed additional solutions from only radiometric and direct landmark data where GM was held fixed at different values. The quality of these solutions' post-fit residuals did not vary significantly, as the post-fit rms values were very close.
Thus, it appears that radiometric and landmark data from OSIRIS-REx are not strong contributors to the determination of the semimajor axis of the orbit. As explained in Section 1 and elsewhere, the semimajor axis and GM are related. This makes the results appear to be insensitive with respect to GM. Yet additional analysis, where we leave out both altimetry and interphase image constraints (but still include intraphase image constraints) and where we start from two different GM values, results in a GM of 4.89284 m 3 s −2 , which is within one standard deviation of our result in Equation (5). Hence, OLA altimetry does not appear to strengthen semimajor axis and GM determination as much as we had hoped for.
However, the inclusion of altimetry did have a positive effect on the solution for other parameters, as well as the statistics of those parameters. Although the correlation between GM and the GDTM scale factor when using interphase image constraints but no altimetry is still very low at 0.0589 (which is very close to the result of 0.0587 when also including altimetry), this correlation increases to 0.884 when neither altimetry nor interphase image constraints are used. Also, without altimetry, the formal standard deviations of the three COF−COM parameters indicate that the Z component is not nearly as well determined as the other two components. In these solutions the adjusted value of the Z component (−21 cm) is also much larger than the other two components (2.3 cm for X and 0.9 cm for Y). Such a change in COF−COM Z occurs for all cases without altimetry (independent of whether the interphase image constraints are used). When altimetry is added, the Z component standard deviation improves much more than the standard deviations of the other two components, bringing them all closer together (the Z component's standard deviation improves from 3.1 to 0.9 cm, while the other components' standard deviation is around 0.5 cm for both cases). Also, the adjusted value of the Z component drops from −21 to −3.6 cm, bringing it more in line with the other two components. All indications are that the addition of altimetry improves the determination of COF−COM, because altimetry constitutes a direct link, with distance information, between spacecraft and asteroid surface. In this case, an offset along the direction of the rotation axis, the Z component, is resolved much better with such information. The X and Y components, perpendicular to the rotation axis, may be better sampled from a polar orbit with regular image cadence, or systematically less affected, and thus could be resolved well from image data alone.
As a result, we find that the presence of the biases in the altimetry data can weaken their contribution to our solutions, in contrast to what was anticipated in simulations (Mazarico et al. 2017). In addition, as shown above, the influence of the altimetry data on the solution is not as straightforward as relationships between GM, semimajor axis, and GDTM scale may lead one to believe. But by estimating separate biases for the different lasers, and by combining these data with the interphase image data and direct landmark data, the inclusion of altimetry data in our analysis in general produces the most consistent set of results: reasonable altimeter biases, small remaining means in nadir and off-nadir altimetry data, the smallest COF−COM values, and, most importantly, a consistent GDTM scale factor result that fits both altimetry and landmark data. We take this to indicate that there are fewer systematic signals remaining in the data for the solutions with altimetry.
Using the values for GM and the GDTM scale factor and their associated errors from Equation (5), we can compute Bennu's bulk density. We list the results of previous studies and our study in Table 6. The shape model we used has a volume of 0.061354 ± 0.00006 km 3 . We apply our scale factor (and its uncertainty) simply as a scaling on the volume (to the power of three), and together with the value for G from CODATA 2018 as listed above, we find a bulk density ρ for Bennu of which has a larger uncertainty than our GM value. A slightly higher bulk density of 1194 ± 3 kg m −3 is reported by Daly et al. (2020). Our density is lower owing to our scale factor, which increases the volume. Using the GM and error from Scheeres et al. (2019), we compute a bulk density error of 1.87 kg m −3 using the OLA GDTM volume and error instead of 3 kg m −3 .
Despite the smaller volume error on the OLA GDTM, the error in bulk density in our results is dominated by the reported volume error, because of our smaller error on GM. Furthermore, if we regard the volume errors reported for the shape models to be mainly dominated by scale errors (considering the dense coverage of Bennu with both images and laser altimetry, making it unlikely that coverage is an issue), then our scale Figure 8. Landmark residuals for independent data from phase Orbital A, for a case using a GDTM scale factor of 1, no rotation angles, no COF−COM correction, no altimetry bias, and standard orientation parameters (labeled "Start values"), and for a case using our newly estimated values for GDTM scale and rotation, orientation, altimetry biases, COF−COM, and GM (labeled "Our solution"). We fix the global parameters to these values but determine the orbit by estimating local parameters. The mean of the landmark residuals is indicated with a dashed-dotted line and one standard deviation with a solid line, with gray for the "Start values" case and black for "Our solution." factor error would supplant the reported volume error (in the density error reported above, we account for both volume error and scale factor error). This would mean a volume error of 0.000066 km 3 and, consequently, a bulk density error that is further reduced to 1.29 kg m −3 . This is smaller than the density error of 1.87 kg m −3 that we find for the values from Daly et al. (2020) because of our smaller GM error. If we use our smaller GM error with the Daly et al. (2020) volume and volume error, we find a density error of 1.17 kg m −3 (dominated by volume error). The slightly different values for GM in these analyses do not change the bulk density value itself. Only the contribution of the GM error to the density error is changed.
The bulk density of Bennu is also consistent with that reported for Ryugu by Watanabe et al. (2019): 1.19 ± 0.02 g cm −3 (which is equal to 1190 kg m −3 ; we use the units and significant digits as reported). Additional analysis by Yamamoto et al. (2020) found a slightly lower value, at 1.18 g cm −3 . The latter depends on both laser altimetry and image tracking data, and their density value is lower owing to a slightly smaller value for GM. The density error is not listed in Yamamoto et al. (2020), but it is presumably the same as that for Watanabe et al. (2019), as the density error is also dominated by volume errors. Analyses that do not take correlations between estimated parameters into account can lead to optimistic error estimates for the bulk density. The inclusion of the image constraints results in a natural decoupling of the parameters, as judged from the correlations, and this results in a more realistic error on the density. When the samples of Bennu are returned to Earth, grain density can be determined accurately. Together with an improved error on the bulk density estimate, this will improve knowledge of the porosity of the asteroid, which has implications for models of the internal structure and thermal properties, among others.

Conclusions and Outlook
We have used DSN radiometric data (Doppler and sequential range) in combination with NavCam 1 image landmarks and OLA laser altimetry data from two phases of the OSIRIS-REx mission: Orbital B (2019 July-August) and Recon B (2020 late January-early February). We used direct altimetry, where the tracks are compared with the OLA GDTM of Bennu . Our landmark data are derived from the same shape model. In addition, we analyzed a differenced image data type that are called image constraints. Image constraints are constructed from observations of identical landmarks in images taken at different times. The distance between the two vectors describing the sight lines from camera to landmark is minimized. Because it is a differenced data type, it is independent of the location of the landmarks, and this means that this measurement type is independent of a priori knowledge of the shape of the body. It is more sensitive to changes that occur over time, such as certain orientation parameters (rotation rate and rotation acceleration). We used image constraints between paired images taken during the same mission phase (intraphase image constraints) and between paired images taken during different mission phases (interphase image constraints). The latter tie together the two mission phases. In a standard multi-arc approach, where data from different spans of time are combined, measurements from each span can be processed separately. For the interphase image constraints, we processed the phases simultaneously because each image constraint measurement has information about each span (such as spacecraft position, pointing, and other parameters).
Our goal with this analysis was to accurately and robustly determine the shape and mass of Bennu. Separate analyses using data from the different phases resulted in different, inconsistent estimates for the asteroid's gravitational parameter GM. Here, we used the interphase image constraints to connect the two phases. We estimated the asteroid's GM and orientation; a scale factor on a GDTM, as well as three rotation angles to ensure that the GDTM and pole were aligned; and the coordinates of the difference between the COF and COM. We also estimated biases in the altimetry data, with separate biases for the HELT and LELT data.
Our results show that the estimated values for the parameters of interest are close when we include or exclude the interphase image constraints, except for GM, which is slightly higher when we exclude the interphase image constraints. However, using the interphase data results in greatly reduced correlations between parameters such as GM and the GDTM scale factor (from 0.705 to 0.059), between the altimeter biases (from 0.79 to 0.35), and between orientation parameters (from −0.90 to 0.015 between the rotation rate and prime meridian, and from 0.716 to −0.016 between the rotation acceleration and prime meridian). Notes. Values between brackets indicate the error in percentages of the parameter value. a GM was not estimated. b Used volume, and 7 times volume error, from Barnouin et al. (2019). c GM based on particles or spacecraft data; no density given. d GM was not estimated. e GM and scale (volume) estimated simultaneously.
We went to great effort to include direct altimetry data from OLA in our large solution. Our goal in doing this was to improve the determination of orbit semimajor axis and GM. Although the inclusion of altimetry data improved our solution in some ways, the effect on GM determination was not as straightforward as initially thought. A stronger source of improvement in GM was the use of multiple mission phases tied together by image constraints. In principle, the inclusion of direct altimetry ranges should have significantly improved determination of orbit semimajor axis and GM. However, the need to estimate altimeter range biases diluted the strength of GM determination through correlations. Even so, inclusion of OLA altimetry strengthened our solutions and also helped in solution assessment. Altimeter ranges greatly improved the determination of the Z component of the parameters that describe the difference between the COM and COF. The analysis of altimeter range residuals binned by pointing angle also provides extra confidence in our GDTM scale factor. Only the combination of altimeter range and landmark data (with the interphase image constraints) results in a GDTM scale factor that fits both direct altimetry and direct landmark data. Independent data in the form of altimetry crossover constraints were also useful in assessing our solutions. Orbits determined using parameters from our solutions that included direct altimetry ranges yielded improvements in altimetry constraint discrepancies and in fits to landmark data from a different orbital phase, Orbital A (2019 January-February). Lastly, we obtain a strong determination of the difference between HELT and LELT biases (HELT-LELT = 25 cm).
From our analysis, we find that Bennu's GM = 4.89256 ± 0.00035 m 3 s −2 , and we find a scale factor for the GDTM of 1.000896 ± 0.00036. This translates to ∼25 cm for a radius of 290 m and thus confirms that the knowledge of Bennu's shape was better than the desired level of 30 cm, which was the requirement for TAG. While some of the methods were developed in collaboration with the navigation team to define a forward path for resolving the shape scale discrepancy, the final results presented here were eventually determined independently of the operational time line and did not directly influence TAG operations.
Using these values, we find a bulk density of 1191.57 ± 1.74 kg m −3 . For the density error, we take the conservative approach and account for both reported volume error and our scale factor error. OLA data reduced the volume error. Our density error is, however, dominated by volume error due to our smaller GM error.
The OSIRIS-REx mission consisted of phases with different orbital characteristics: the sensitivities of different orbits to parameters such as GM and scale already can help to decorrelate parameters, yet we showed that this is only to some extent, and not automatically always the case. Specifically, our analysis shows the strength of including the interphase data: they connect two different mission phases, and they decorrelate important parameters, increasing the robustness of their estimates. The data returned by the OSIRIS-REx mission over 2 yr of proximity operations and multiple, distinct orbit configurations have shown the value of such data sets, particularly cross-phase observations. Although not all missions will have different orbital characteristics during their mission span, there are plenty of scenarios, such as multiple flybys, where this is more or less the case. Regardless, our analysis shows that measurements that connect observations separated in time will improve the determination of GM, parameters related to shape, and parameters that describe change over time.