Measurability of the Heliocentric Momentum Enhancement from a Kinetic Impact: The Double Asteroid Redirection Test (DART) Mission

This article presents the Gauss–Radau Small-body Simulator (grss), ¹¹ an open-source small-body orbit propagation and orbit determination tool. The orbit determination algorithm in the grss library is based on the least-squares filter, the algorithm for which comes from Vallado (2022). The goal of a least-squares filter is to minimize the cost function,

$$\begin{eqnarray}&&J={{\boldsymbol{b}}}^{T}{\boldsymbol{W}}{\boldsymbol{b}},\end{eqnarray} \tag{ 2 }$$

where J is the cost and b is the N-dimensional vector of residuals of each observation. The ith element of this vector is

$$\begin{eqnarray}&&{{\boldsymbol{b}}}_{i}={{\boldsymbol{O}}}_{i}-{{\boldsymbol{C}}}_{i}.\end{eqnarray} \tag{ 3 }$$

O is the vector of observables, C is the vector of computed observations, and W is the N × N weight matrix. The observation data used in grss come from two sources: (1) optical astrometry from the International Astronomical Union's Minor Planet Center (MPC), ¹² and (2) radar astrometry from the JPL Small-Body Radar page. ¹³ Additional measurements, such as stellar occultations and spacecraft measurements, can be ingested using the MPC's ADES format. ¹⁴

The computed measurements are acquired by propagating an initial N_par-dimensional vector of solve-for parameters (where N_par is the number of fitted parameters), X ₀ consisting of the cometary orbital elements (Milani & Gronchi 2009), any nongravitational parameters according to the Farnocchia et al. (2013) transverse acceleration model, and, in our case, β_⊙ to the epoch of each observation. This initial guess is updated through a process called differential corrections. The correction at every iteration is calculated using the equation

$$\begin{eqnarray}&&{\boldsymbol{\Delta }}{\boldsymbol{X}}={\boldsymbol{P}}{{\boldsymbol{A}}}^{T}{\boldsymbol{W}}{\boldsymbol{b}},\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{P}}={({{\boldsymbol{A}}}^{T}{\boldsymbol{W}}{\boldsymbol{A}})}^{-1}$ is the covariance matrix of the solution and A is the gradient of the observations with respect to the solve-for parameters. grss calculates this gradient using first-order central differences. This corrected state (calculated as X ₀ + Δ X ) is computed until the change in nominal solution is within a prescribed tolerance. The solution is said to have converged at that point. This classical process of nonlinear weighted least-squares estimation can be traced back to the work of C. F. Gauss (Gauss 1809).

The grss library uses the ias15 algorithm, which is a 15th-order integrator based on Gauss–Radau spacings, which also gives the library its name (Everhart 1985; Rein & Spiegel2014). The propagator submodule within grss uses a force model that incorporates the following:

• 1.  
  Newtonian point-mass interactions (Newton 1687);
• 2.  
  Einstein–Infeld–Hoffman point-mass relativistic effects (Einstein et al. 1938; Moyer 2003);
• 3.  
  Zonal harmonics accelerations from the Sun and Earth (Laplace 1798; Kaula 1966); and
• 4.  
  A₁, A₂, A₃ nongravitational acceleration terms in the heliocentric radial, transverse, and normal directions, respectively (Marsden et al. 1973).

The Newtonian accelerations are calculated for the Sun, planets, the Moon, Pluto, and the Big16 main-belt asteroid perturbers from the DE441 ephemerides. On the other hand, the relativistic accelerations are only calculated for the Sun, planets, the Moon, and Pluto. Marsden et al. (1973) define the nongravitational acceleration terms. However, in this work, only the acceleration in the transverse direction is calculated. This is done using an inverse square expression of the form A_T = A₂(1 au/r)^d , where r is the heliocentric distance and d = 2 is chosen to match the level of absorbed solar radiation (Farnocchia et al. 2013). It is possible to vary this exponent for a better orbit fit as was done by Chesley et al. (2014). However, in the case of Didymos, retaining d = 2 already provides an excellent orbit fit to all the observations.

Before the optical astrometry is passed through the orbit determination code, it needs to be treated to account for star catalog biases and weighted appropriately to ensure statistical reliability of the orbit solution given by the least-squares estimator. We use the astrometric catalog debiasing scheme from Eggl et al. (2020) to account for biases in the star catalogs used to reduce the optical astrometry reported to the MPC.

In addition to the debiasing, the optical astrometry acquired from the MPC also needs to be assigned weights before it is passed through the filter. We use the Vereš et al. (2017) weighting scheme to assign weights to observations depending on the observatory and the star catalog used for each observation. We also adopt the Vereš et al. (2017) rule of inflating the weights by a scaling factor of $\sqrt{{N}_{\mathrm{obs}}/4}$ for each observation if the number of observations from the same observatory on the same night is N_obs > 4. Once these optical observations are debiased and properly weighted, they are ready for the orbit filter. The radar astrometry uncertainty acquired from the JPL radar API is accepted at face value.

During the filtering process, one more algorithm is applied to the residuals. The grss orbit filter also has the capability to perform automatic outlier rejection according to the work of Carpino et al. (2003). This algorithm rejects and readmits optical measurements in the orbit fit based on a rejection threshold χ_rej and a recovery threshold χ_rec. We set χ_rej = 3 and χ_rec = 2.8 for grss. The χ values for each optical measurement pair are calculated as

$$\begin{eqnarray}&&{\chi }_{i}^{2}={{\boldsymbol{b}}}_{{\boldsymbol{i}}}{{\boldsymbol{\Gamma }}}_{{{\boldsymbol{b}}}_{{\boldsymbol{i}}}}^{-1}{{\boldsymbol{b}}}_{{\boldsymbol{i}}}^{T},\end{eqnarray} \tag{ 5 }$$

where b _i is the vector of residuals for the ith optical observation pair (R.A. and decl.) and ${{\boldsymbol{\Gamma }}}_{{{\boldsymbol{b}}}_{{\boldsymbol{i}}}}$ is the residual covariance matrix of that pair of residuals, calculated as

$$\begin{eqnarray}&&{{\boldsymbol{\Gamma }}}_{{{\boldsymbol{b}}}_{{\boldsymbol{i}}}}={{\boldsymbol{\Gamma }}}_{{\boldsymbol{i}}}\pm {{\boldsymbol{A}}}_{{\boldsymbol{i}}}{\boldsymbol{P}}{{\boldsymbol{A}}}_{{\boldsymbol{i}}}^{T}.\end{eqnarray} \tag{ 6 }$$

Here A _i is the 2 × N_par array of the gradient for the ith observation (the first row corresponds to the partial derivative of the R.A. with respect to the fitted parameters, and the second row is for the decl. measurement), and Γ _i is the 2 × 2 covariance matrix of the same optical observation. The second term on the right-hand side of the equation has a negative sign when the observation is included in the fit and a positive sign when it is rejected. In our implementation, the initial guess given to the filter is allowed to converge without any outlier rejection on the first pass. Once a converged solution is found, outlier rejection is turned on, and the filter will successively remove and/or recover any optical measurements with too high of a χ value. The weighted residuals are once again allowed to converge with this outlier rejection scheme turned on, and the converged solution after this stage is the output of the grss orbit filter.

In addition to the data from MPC and the JPL radar astrometry database, grss also has the capability to generate synthetic observations for the purposes of covariance analysis. These synthetic data are generated from a reference trajectory and can be generated with or without any constant biases or Gaussian noise around the reference trajectory. The synthetic observation feature is heavily used in this work to investigate the number and type of observations needed to measure β_⊙.

Before the grss orbit propagator and orbit filter can be used to estimate β_⊙, they must be validated against tools that have been demonstrably successful in the past. We individually tested the propagator and orbit determination code with the heliocentric trajectory of the (65803) Didymos system barycenter. This was a natural choice, since the Didymos system was also the target of the DART mission. Additional testing with different asteroids was also done and is shown in the examples available on the grss documentation. ¹⁵

To validate the propagator, we took the latest Didymos orbital solution from the JPL Small-Body Database ¹⁶ (Solution 204 at the time of writing) and propagated it for 3000 days into the future. The initial solution is shown in Table 1. The difference in the final propagated position when compared to JPL's internal small-body propagation code was 2.4 m. Another longer-term propagation of 250 yr exhibited differences on the order of 1 km, which can be attributed to differences in the numerical scheme used to integrate the equations of motion.

Testing the orbit determination code is inherently more comprehensive, as it accumulates errors from both the integration and the observation model. The observation data set for Solution 204 includes a diverse set of observation types. There are more than 3000 optical measurements from ground-based observatories that include optical astrometry from telescopes around the world, radar delay measurements from the Goldstone Solar System Radar facility and Arecibo Observatory, optical measurements from the DART spacecraft before impact, and ground-based stellar occultation measurements. These stellar occultation measurements provide high-quality R.A.–decl. astrometry, and therefore very accurate observation models must be used when fitting these data.

Table 2 shows the orbit solution produced using grss. The grss solution and the JPL Solution 204 show excellent agreement with each other, with all orbital elements within 1σ of each other. Additionally, the Bhattacharya distance, which is a measure of the distance between two statistical distributions, is only 0.060. The complement to this distance is the Bhattacharya coefficient, which measures the similarity between two statistical distributions, which is 0.942 in this case. These two metrics indicate that the two solutions are exceedingly close. Given a similar input, orbit solutions computed by the grss library are providing results of comparable quality to those of the JPL solar system dynamics group. Figure 2 contains further support for this statement. Here, the postfit residuals of the optical measurements used by grss are shown. The red circles around the scatter points indicate that those residuals were automatically removed by the outlier rejection scheme. We see that most observations with high residuals are automatically rejected. The few that remain are weighted such that, even with the high residuals, they do not have a meaningful effect on the final solution.

Zoom In Zoom Out Reset image size

The remaining differences in the solutions can be explained by the fact that the JPL solutions sometimes incorporate manual weighting of the observations on top of the Vereš et al. (2017) rule. Additionally, the JPL solution assumes a 1 s uncertainty on the timing of all observations. This is not modeled by grss, thus resulting in a lower uncertainty than if time uncertainty were taken into account. There are also some differences introduced by the fact that JPL propagates analytical partials along with the state and grss uses numerical derivatives. This leads to small differences in the covariance matrix between the two solutions. However, Tables 1 and 2 show that the uncertainties are quite close to each other. Additionally, Figure 3 shows the correlation factors for each of the two orbit solutions. The two sets of correlation factors are in good agreement with each other. This indicates that it is not just the 1σ uncertainties on the diagonals of the covariance matrices that agree well, but it is both covariance matrices that are very similar to each other.

Zoom In Zoom Out Reset image size

grss does not just perform well at fitting optical data. Figure 4 shows the postfit χ values of each observable in the fit. We see that the radar delay values (shown in green) are closely packed around 0, indicating that the grss observation model does not have trouble fitting radar data either. The radar data from 2003 have lower rms values than the ones in 2022 because they were regenerated using a shape model and therefore have lower uncertainties on the measurement (Naidu et al. 2020). Due to the lower uncertainties, the filter assigns a tighter weight to the 2003 delay points, which leads to better unnormalized residuals.

Zoom In Zoom Out Reset image size

A consequence of the DART impact on the Didymos system barycenter is a momentum transfer of magnitude β_⊙ times the momentum of the DART spacecraft at impact. Figure 5 depicts the direction of this impulse and its associated formal 1σ uncertainty using the observation scenario from Section 3.5, as well as the directions of the velocity vector of the DART spacecraft and the Didymos system barycenter at the time of impact. The blue vector represents the direction of the velocity of the DART spacecraft (relative to the Didymos system barycenter) at the time of impact, and the red arrow represents the direction of the velocity of the system barycenter relative to the Sun at the same epoch. The uncertainty of the delta-v imparted by DART is shown by the green ellipsoid. More details on how this ellipsoid was calculated are presented at the end of Section 3.5.

Zoom In Zoom Out Reset image size

The minor axis of the ellipsoid is significantly smaller than the major and intermediate axes. Additionally, the minor axis is oriented along nearly the same direction as the DART spacecraft velocity at impact. In this work, we therefore make the assumption that the resulting impulse is instantaneously applied in the same direction as the DART spacecraft's velocity relative to the Didymos system barycenter at impact. Fixing the direction of the impulse in a direction that cuts through the smallest axis of the ellipsoid helps the filter when estimating β_⊙.

Since this work focuses on measurability of β_⊙, rather than on an actual estimate of β_⊙, this assumption of alignment of the DART velocity vector and the impulse to the barycenter is justified. Measurement of β_⊙ from the DART impact would need a better estimate of this impulse direction. This direction would then be informed by high-fidelity ejecta dynamics simulations and observations of the ejecta (Fahnestock et al. 2022; Li et al. 2023). Conversely, accurately estimating the three-dimensional β_⊙ properties can help better constrain the ejecta plume as well.

The first question that needs to be answered when looking at β_⊙ measurability is whether it can be estimated with the already-existing observation data. Therefore, we took all of the Didymos astrometry up until the time of writing (2023 August) and attempted to estimate β_⊙. This data set included all the observations in the JPL Solution 204 data set and more than 2000 additional ground-based optical astrometry points from 2022 August to 2023 June from the MPC.

Results in this case showed a β_⊙ estimate of −4.8 ± 16.5, indicating a signal-to-noise ratio (S/N) of about 0.3. This shows that the existing data are not enough to give a statistically significant estimate. Since β_⊙ is the link between the pre- and post-impact orbits of the Didymos system, both of them need to be well constrained to obtain a significant β_⊙ estimate. Since the post-impact observation arc is less than 9 months long at the time of writing this article, a β_⊙ estimate is not possible at this time. Therefore, we look at the future observation scenarios needed to get a reliable β_⊙ estimate.

In the first observation scenario, we consider the possibility of a handful of new stellar occultation measurements at a monthly cadence between 2024 June and October. The 2024 June start time was chosen since there is a good chance of getting successful occultation measurements in the second half of 2024. This prediction is based on the density of star background that would enable an occultation measurement, regardless of the duration and magnitude of the occulted star. Additionally, the uncertainty of the current orbit of Didymos propagated to the beginning of its observability period in 2024 is only ≈1 Didymos body radius. The 2024 October cutoff was chosen because that is the planned launch date of ESA's Hera mission (Michel et al. 2022), and we wanted to explore the possibility of estimating β_⊙ before the launch.

These synthetic future stellar occultation measurement uncertainties were assumed to be at 2.5 mas. This is slightly higher than the average uncertainty of the existing Didymos occultation measurements (Chesley et al. 2023), to ensure a conservative estimate. These observations were generated from the JPL Solution 204 for Didymos as a reference trajectory. A reference value of β_⊙ = 3 was set as the "true" value to assess the accuracy of the estimates returned by the orbit filter. This assumption is justified because the nominal β value (when not accounting for Dimorphos mass uncertainties) for the DART impact is about 3.6. β_⊙ must be lower than the β value since only a subset of ejecta that escaped Dimorphos could have escaped the system's Hill sphere.

Additionally, even though it is not statistically significant, the β_⊙ estimate from Section 3.1 is less than 0.5σ away from this assumption. This implies that β_⊙ = 3 is not an unreasonable assumption to make. It is, however, important to stress that this is not the β_⊙ value from the DART impact, but a reference value to assess how well we could measure β_⊙ under different scenarios. The β_⊙ estimate in this case was 3.2 ± 0.7, at an S/N ∼ 4.9. Therefore, occultations in the second half of 2024 already give a long-enough observation arc for a reliable β_⊙ estimate before the Hera launch.

After analyzing the effects of the 2024 occultations on the β_⊙ S/N, we decided to add more occultations in 2027, when the Hera spacecraft will be at the Didymos system. Therefore, five additional monthly occultations from 2027 February to June were added to the observation data set for Case 1. The β_⊙ estimate in this case was 3.0 ± 0.4, at an S/N ∼ 8. Therefore, the 2027 occultations increased the S/N by a factor >1.5. However, all five of these additional occultations might not be possible in early 2027, so we decided to explore another case.

After the Hera mission's planned launch in 2024 and a 2 yr cruise, the spacecraft would arrive at the Didymos system in late 2026 and would begin operations in early 2027. As part of these operations, the spacecraft would regularly communicate with radar ground stations on Earth. Some of these tracking passes could then be converted into pseudorange measurements of the asteroid barycenter, as was done with the OSIRIS-REx mission at Bennu (Farnocchia et al. 2021). Case 3 considers converting some of these Hera tracking measurements into pseudorange measurements of the Didymos system barycenter at a monthly cadence from 2027 February to June.

These pseudorange measurements are added to existing observations without any new occultations. The expected uncertainties are ∼5 m, based on detailed numerical simulations of the Hera Radio Science experiment (Zannoni et al. 2018). However, all synthetic pseudorange measurements of the barycenter of the Didymos system in this study are weighted at an uncertainty of 15 m to ensure a conservative estimate. This is also conservative compared to the 2 m uncertainties of the Farnocchia et al. (2021) OSIRIS-REx pseudorange points and should not produce overly optimistic S/N values. In this case, the β_⊙ estimate is 3.0 ± 0.2, with S/N ∼ 17. This increase in S/N is due to the fact that the Hera pseudorange measurements provide extremely accurate astrometry, and these measurements are taken at a point where the post-impact arc is almost 5 yr long.

The last observation scenario considered for estimating the β_⊙ for DART is the superset of the three future observation scenarios considered up to this point. These data include all existing observations and assume that new occultation measurements in 2024 and 2027 will be successful and that the Hera mission will provide pseudorange measurements for the Didymos system barycenter. In this scenario, the estimate of β_⊙ is 3.0 ± 0.2, with an S/N ∼ 17. This is the best estimate of β_⊙ (negligibly better S/N than Case 3), as it is based on the largest set of high-precision observations.

The S/N value in this case has another point of validation. In Figure 5, the ellipsoid of the delta-v uncertainty is calculated using this observation scenario. A separate capability of grss that can be used to estimate unconstrained ΔV vectors was used to obtain the 3 × 3 subcovariance matrix corresponding to the DART ΔV. The eigenvalues of this covariance matrix yield the uncertainty of this ΔV (and therefore the dimensions of the ellipsoid), and the eigenvectors give us the orientation of the uncertainty ellipsoid.

Using the estimate of this unconstrained ΔV vector and its uncertainty, we find that the theoretical maximum S/N value of the ΔV magnitude for this observation scenario is 20.3. We consider the S/N of the magnitude of the ΔV vector instead of the components since this is similar to estimating the scalar β_⊙ value. (The S/N for the components of the ΔV vector in such a scenario would be much smaller owing to a lack of directional constraint.) Additionally, there is a small 16° offset between the eigenvector corresponding to the maximum S/N and the assumed direction of the escaping ejecta momentum in this work. Due to this offset, the maximum achievable S/N for this observation scenario reduces to 19.5. This maximum possible S/N is only 15% higher than the one we achieve with the assumption about the direction of the escaping ejecta momentum, and therefore it serves as a sanity check for the β_⊙ S/N values we get from the least-squares estimator. Therefore, the assumption about the direction of the escaping ejecta is not significantly affecting the measurability of β_⊙.

Table 4 shows the summary of the simulated estimates for β_⊙ and the A₂ transverse nongravitational acceleration parameter for all four cases described in this section. We see the expected progression of the S/N as more observations of the Didymos system are added from Case 1 to Case 4. Additionally, we present the fitted values and S/Ns for the fit of the A₂ parameter along with β_⊙ for each of the cases. This is presented because the two parameters are highly correlated (as shown by the last column of Table 4). The Yarkovsky effect causes a secular semimajor-axis drift that has been detected in many other asteroids, with Bennu having the best detection owing to the OSIRIS-REx mission (Farnocchia et al. 2021). On the other hand, the DART impact impulsively changed the orbit of the Didymos system barycenter. So when estimating both parameters simultaneously, the filter has to disentangle the semimajor-axis drift from the Yarkovsky effect, $\dot{a}$, from the semimajor-axis change Δa due to the DART impact.

Table 4. DART β_⊙ and A₂ Estimation Summary for Various Simulated Future Observation Scenarios

Scenario	β⊙	${\sigma }_{{\beta }_{\odot }}$	S/N${}_{{\beta }_{\odot }}$	A2 (au day−2)	${\sigma }_{{A}_{2}}$ (au day−2)	S/N${}_{{A}_{2}}$	Corr (β⊙, A2)
Case 1	3.21	0.66	4.85	−1.03E−14	7.21E−16	14.21	0.57
Case 2	3.03	0.39	7.84	−1.03E−14	6.99E−16	14.75	0.67
Case 3	3.01	0.18	16.75	−1.04E−14	4.19E−16	24.78	1.00
Case 4	3.03	0.18	17.01	−1.03E−14	4.23E−16	24.43	1.00

Note. The target value to be retrieved in these simulations is β_⊙ = 3, A₂ = −1.04E−14 au day^−2.

Download table as:  ASCIITypeset image

To further illustrate this correlation, Figure 6 shows the 3σ covariance ellipses generated using the 2 × 2 subcovariance matrix for the four cases. We see the uncertainty ellipses get smaller as the S/N on both β_⊙ and A₂ increases when more observations are added in each scenario. The correlation between the two parameters in Cases 1 and 2 is 0.568 and 0.670, respectively. Once the Hera pseudorange measurements are added to the estimate in Cases 3 and 4 (highlighted in Figure 7), the correlation between β_⊙ and A₂ jumps to 0.996 in both scenarios.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Additionally, Figure 8 shows the complete correlation matrix of the orbit fit for Case 4. The strong correlation between β_⊙ and A₂ is even more apparent in this depiction. This highlights the importance of having an accurate picture of the acceleration due to the Yarkovsky effect when it comes to estimating β_⊙ for kinetic impact-based missions. In the case of Didymos, as well as many other asteroids, such as Apophis, an accurate estimate of the Yarkovsky acceleration has been enabled by either radar delay measurements or stellar occultations, or both. It is important to note that all near-Earth asteroids that have had successful occultation campaigns at the time of writing have been observed first by radar facilities. For this reason, when considering β_⊙ estimation scenarios for a future KI mission, we consider these two types of measurements both before and after impact.

Zoom In Zoom Out Reset image size