Constraining the interiors of asteroids through close encounters

ABSTRACT

Knowledge of the interior density distribution of an asteroid can reveal its composition and constrain its evolutionary history. However, most asteroid observational techniques are not sensitive to interior properties. We investigate the interior constraints accessible through monitoring variations in angular velocity during a close encounter. We derive the equations of motion for a rigid asteroid’s orientation and angular velocity to arbitrary order and use them to generate synthetic angular velocity data for a representative asteroid on a close Earth encounter. We develop a toolkit AIME (Asteroid Interior Mapping from Encounters) which reconstructs asteroid density distribution from these data, and we perform injection-retrieval tests on these synthetic data to assess AIME’s accuracy and precision. We also perform a sensitivity analysis to asteroid parameters (e.g. asteroid shape and orbital elements), observational setup (e.g. measurement precision and cadence), and the mapping models used. We find that high precision in rotational period estimates (≲0.27 s) is necessary for each cadence, and that low perigees (≲ 18 Earth radii) are necessary to resolve large-scale density non-uniformities with uncertainties of |$\sim 0.1{{\ \rm per\ cent}}$| of the local density under some models.

methods: data analysis, minor planets, asteroids: general, planets and satellites: interiors

1 INTRODUCTION

Over the past 20 yr, the increase in quantity and quality of sensitive all-sky surveys has prompted the discovery of numerous asteroids. Such advances have been made via ground-based surveys such as the Catalina Sky Survey (Larson et al. 1998), Pan-STARRS (Kaiser et al. 2002), and the Lincoln Near-Earth Asteroid Research project (LINEAR; Stokes et al. 2000), as well as space-based instruments such as the Wide-field Infrared Survey Explorer (WISE) mission (Wright et al. 2010). Many of these asteroids are relatively small, but some are kilometre-sized and a few are predicted to closely encounter Earth or other planets in the near future. More encounter candidates are likely to be discovered by new efforts such as the Large-aperture Synoptic Survey Telescope (LSST; Tyson 2002). Their encounters can then be monitored by global ground-based networks such as the Las Cumbres Observatory (LCO; Brown et al. 2013). Such ground-based monitoring is typically used to derive the rotation period of an asteroid and its surface properties (see e.g. Devogèle et al. 2021).

Since the tidal torque acting on an asteroid during an encounter depends on the interior mass distribution, the careful monitoring of angular velocity variations during an encounter also presents a window into the interior properties of asteroids. The gravitational two-body system has been studied in the context of tidal torque to different orders and with several different methods (Paul 1988; Scheeres et al. 2000; Ashenberg 2007; Boué & Laskar 2009; Hou et al. 2017). Further theoretical studies showed that the tidal torque, observed through angular velocity perturbations, is sensitive to asteroid interior density distribution (Scheeres, Marzari & Rossi 2004; Naidu & Margot 2015; Makarov et al. 2022).

Angular velocity perturbations have been observed and used to extract asteroid properties in several cases, including for the 2013 encounter of (367943) Duende with Earth (Benson, Scheeres & Moskovitz 2020; Moskovitz et al. 2020), and asteroid binaries (3905) Doppler and (617) Patroclus (Berthier et al. 2020; Descamps et al. 2020). Orbital and physical properties, including moment of inertia (MOI) ratios, have also been extracted for asteroid 99942 Apophis, discovered on 2004 June 19 by R. A. Tucker, D. J. Tholen, and F. Bernardi (Giorgini et al. 2005, 2008; Smalley et al. 2005), and on target to encounter Earth in 2029 (Yu et al. 2014; Hirabayashi, Kim & Brozović 2021; Lee et al. 2022; Valvano et al. 2022). However, density distribution features beyond the MOI ratios have not yet been extracted for any asteroid encounters. More research is needed to study in what cases these effects are observable, and what factors generally inhibit observation of these new features. It seems pivotal to augment previous work on the affect of tidal torque on Apophis’ angular velocity in particular (Souchay et al. 2014, 2018) so that upcoming observations may constrain these properties and thus improve our predictions.

We address this community need by developing a methodology to translate (1) time series of asteroid angular velocity data into constraints on density moments and (2) constraints on density moments into constraints on an asteroid’s density distribution. Other techniques, such as measurement of tidal distortion (Richardson, Bottke & Love 1998), impact or seismometry experiments (Richardson et al. 2005), or gradiometry (Carroll & Faber 2018), may additionally constrain the density distribution. In Section 2, we introduce the analytical and numerical fundamentals of this methodology. There, we describe a simulation used to integrate the equations of motion and produce synthetic data of angular velocity over time, followed by a Markov chain Monte Carlo (MCMC) fit process which extracts density moments from the fit data. We then describe two methods to generate full density distributions from the density moments. In Section 3, we present the results of a series of injection-retrieval tests demonstrating the extent to which the properties of an asteroid chosen to generate synthetic spin data can be retrieved via our methodology. Finally, in Section 4, we assess the sensitivity of these constraints to various physical, observational, and methodological parameters to provide guidance for monitoring upcoming close encounters.

2 METHODS

In the following mathematical model for the influence of tidal torque on asteroid encounters, we assume for simplicity that (1) the central body is much more massive than the asteroid, (2) both are rigid, (3) there are no distant perturbing objects, and (4) the asteroid is in a short-axis, non-tumbling rotational state before the encounter. All of these assumptions except (2) can be relaxed without drastic changes to the model.

The only properties of an asteroid’s density moments that affect tidal torque interactions are its ‘density moments’, defined here as

$$\begin{eqnarray*} K_{\ell m} = \frac{a_\mathcal {A}^{2-\ell }}{I_\mathcal {A}} \int _\mathcal {A} d^3 r \rho _\mathcal {A}(\boldsymbol r) R_{\ell m}(\boldsymbol r). \end{eqnarray*}$$

(1)

These are complex, unitless quantities. |$\rho _\mathcal {A}(\boldsymbol r)$| is the asteroid density distribution and R_ℓm are the regular solid spherical harmonics (see Supplementary Material A.2 for details). The integral is computed over the entire asteroid mass, denoted |$\mathcal {A}$|⁠, and d³r indicates the three-dimensional volume element throughout the paper. |$I_\mathcal {A}$| denotes a MOI scale defined as

$$\begin{eqnarray*} I_\mathcal {A} = \int _\mathcal {A} d^3 r \rho _\mathcal {A}(\boldsymbol r) r^2 \end{eqnarray*}$$

(2)

while |$a_\mathcal {A}$| is the length-scale

$$\begin{eqnarray*} a_\mathcal {A}^2 = \frac{1}{V_\mathcal {A}} \int _\mathcal {A} d^3 r r^2 \end{eqnarray*}$$

(3)

where |$V_\mathcal {A}$| is the asteroid volume. We call these MOI and length-scales in part because they obey |$I_\mathcal {A} = \mu _\mathcal {A} a_\mathcal {A}^2$| for uniform asteroids where |$\mu _\mathcal {A}$| is the mass of the asteroid. Note that |$a_\mathcal {A}$| is a function only of the surface of the asteroid, so that |$a_\mathcal {A}$| is known if the surface is observed.

The tidal torque experienced by an asteroid is

$$\begin{eqnarray*} \boldsymbol \tau &= & G\frac{I_\mathcal {A}I_\mathcal {B}}{2 a_\mathcal {A}^2a_\mathcal {B}^2}\left[\sum _{\ell m} a_\mathcal {B}^\ell J_{\ell m} \sum _{\ell ^{\prime } m^{\prime }}a_\mathcal {A}^{\ell ^{\prime }}S_{\ell +\ell ^{\prime }, m + m^{\prime }} (\boldsymbol D)^{*} (-1)^{\ell ^{\prime }}\right.\nonumber \\ && \times \left.\sum _{m^{\prime \prime }=-\ell ^{\prime }}^{\ell ^{\prime }} \sqrt{\frac{(\ell ^{\prime }-m^{\prime \prime })!(\ell ^{\prime }+m^{\prime \prime })!}{(\ell ^{\prime }-m^{\prime })!(\ell ^{\prime }+m^{\prime })!}} \mathcal {D}^{\ell ^{\prime }}_{m^{\prime }m^{\prime \prime }}(\alpha , \beta , \gamma)^{*} \right. \nonumber \\ && \times \Big ((i\boldsymbol {\hat{x}} - \boldsymbol {\hat{y}})(\ell ^{\prime }-m^{\prime \prime }+1)K_{\ell ^{\prime },m^{\prime \prime }-1} \nonumber \\ && +(i\boldsymbol {\hat{x}}+\boldsymbol {\hat{y}})(\ell ^{\prime }+m^{\prime \prime }+1)K_{\ell ^{\prime },m^{\prime \prime }+1}+2im^{\prime \prime }\boldsymbol {\hat{z}} K_{\ell ^{\prime }m^{\prime \prime }}\Big) \Bigg ], \end{eqnarray*}$$

(4)

where |$\boldsymbol D$| is the position of the asteroid; α, β, and γ are z − y − z Euler angles expressing the orientation of the asteroid; S_ℓm are the irregular solid spherical harmonics; and |$\mu _\mathcal {B}$| and |$a_\mathcal {B}$| are the mass and radius of the central body while J_ℓm are the density moments of the central body. Equation (4) is derived in Supplementary Material A.4 via a novel derivation and is accurate to arbitrary order in D.

Since it is the angular acceleration of the asteroid that is observable, rather than the torque applied, we also compute the MOI of the asteroid around the principal axes:

$$\begin{eqnarray*} I_x &=& \frac{2}{3} I_\mathcal {A} \left(K_{20} - 6 K_{22} + 1 \right)\nonumber \\ I_y &=& \frac{2}{3} I_\mathcal {A} \left(K_{20} + 6 K_{22} + 1 \right)\nonumber \\ I_z &=& \frac{2}{3} I_\mathcal {A} \left(-2K_{20} + 1 \right). \end{eqnarray*}$$

(5)

Note that all moments of inertia are proportional to |$I_\mathcal {A}$|⁠. Equation (4) indicated that tidal torque was also proportional to |$I_\mathcal {A}$|⁠. The Euler equations, giving angular acceleration in terms of moment of inertia and torque (Supplementary Material equation A.14), show that the angular acceleration of the asteroid is therefore independent |$I_\mathcal {A}$|⁠. Thus, the observables do not depend explicitly on |$I_\mathcal {A}$|⁠.

Throughout the paper, we refer to the ‘inertial frame’ (the frame in which the orbit is fixed) and the ‘body-fixed frame’ (the frame in which the asteroid is fixed and in which K_ℓm, |$I_\mathcal {A}$|⁠, and |$a_\mathcal {A}$| are computed), which are also defined in Supplementary Material A.1.

2.1 Simulation design

We built a simulation of an asteroid’s rotational state during a close encounter with a central body. This simulation requires as initial data (1) the orbital parameters of the asteroid r_p (perigee distance) and v_∞ (hyperbolic excess velocity); (2) the cadence of angular velocity observation Δt; (3) the central body moments J_ℓm, mass |$\mu _\mathcal {B}$|⁠, and radius |$a_\mathcal {B}$|⁠; (4) the initial asteroid angular velocity in the inertial frame |$\boldsymbol \Omega _0$|⁠; (5) the asteroid length |$a_\mathcal {A}$|⁠; and (6) the asteroid’s density moments K_ℓm and initial Euler angle γ₀. All parameters except (6) are assumed to be known to high accuracy. For example, |$\boldsymbol \Omega _0$| could be extracted from pre-encounter light-curve data and |$a_\mathcal {A}$| from a model for the asteroid’s surface using radar. The other two initial Euler angles α₀ and β₀ are required to be zero by the assumption of no initial tumbling.

We begin our simulation at D = 10r_p. Since the leading order of equation (4) yields τ ∝ D⁻³, this corresponds roughly to a torque of 10⁻³ times the maximum torque at perigee. Unless otherwise indicated, the simulation is terminated at D = 10r_p as well. With the simulation inputs specified, the equations of motion are integrated via the Runge–Kutta fourth-order method with a variable time-step, manually chosen to limit integration error to only ∼100 times the floating point numerical error.

We extract angular velocity from this simulation assuming that angular velocity is observable. If orientation data are more readily available instead, the same simulation may be used to generate orientation data. Supplementary Material C describes the consequences of this change, which are relatively small.

2.2 Uncertainty model

To add noise to data generated via the above simulation, we use the following uncertainty model. Each asteroid spin vector |$\boldsymbol \Omega$| is assumed to be uncorrelated with other spin vectors, and we model uncertainty in the orientation and in the period as also uncorrelated. Consider a true spin vector |$\boldsymbol \Omega ^{*}$|⁠. For the sake of description, we work in coordinates in which |$\boldsymbol \Omega ^{*} \parallel \boldsymbol {\hat{z}}$|⁠. Then, expressing the observed spin vector |$\boldsymbol \Omega$| in spherical coordinates, we draw the polar angle from a normal distribution with standard deviation σ_θ centred on zero and the azimuthal angle from a uniform distribution. We also draw the ratio Ω/Ω* from a lognormal distribution centred on one, with width σ_P/P_ω, where P_ω = 2π/Ω is the period of the asteroid. Explicitly, the probability density function (PDF) of P_ω is

$$\begin{eqnarray*} P(\rho) = \frac{1}{\rho \sqrt{2\pi (\sigma _P / P_\omega)^2}} \exp \left(-\frac{\ln ^2(\sigma _P / P_\omega)}{2(\sigma _P / P_\omega)^2} \right). \end{eqnarray*}$$

(6)

See Fig. 1 for an illustration of the uncertainty model. A lognormal distribution is chosen such that P_ω > 0, but since σ_P/P_ω ≪ 1 typically in our analysis, the probability distribution is essentially Gaussian.

Figure 1.

Diagram in the inertial frame of the uncertainty model used to define the probability that the true spin vector |$\boldsymbol \Omega ^{*}$| should be observed as |$\boldsymbol \Omega$|⁠. The angular uncertainty on θ and period uncertainty on |$|\boldsymbol \Omega |$| are treated separately.

The log likelihood resulting from this uncertainty model is (excluding additive constants)

$$\begin{eqnarray*} \ln \mathcal {L} &=& -\frac{1}{2}\sum _{i = 0}\left [\frac{\cos ^{-1} (\boldsymbol \Omega _i^{*} \cdot \boldsymbol \Omega _i/(\Omega _i^{*} \Omega _i))^2}{\sigma _\theta ^2}\right. \nonumber \\ && \left. +\frac{\ln \left(\Omega _i /\Omega _i^{*} \right)^2}{(\sigma _P / P)^2} + 2\ln \frac{\Omega _i}{\Omega _i^{*}}\right ]. \end{eqnarray*}$$

(7)

where Ω_i is the ith spin vector in the data set.

This model was chosen because it separates spin pole and period uncertainty. Therefore, if one is more precisely determined by measurement, σ_θ and σ_P/P can be adjusted separately in accordance.

2.3 Extracting density moments from spin data

Given synthetic data, an Affine Invariant MCMC Ensemble sampler was used to generate posterior probability distributions (PPDs) from flat priors. We use the python implementation emcee (Foreman-Mackey et al. 2013). Our parameters were γ₀, K₂₀, K₂₂, and K_3m (10 in total), bounded by |γ₀| < π/4, and the K_2m bounds given in Supplementary Material equation A.5. Computing K_3m for many positive density distributions revealed a |K_3m| ≲ 0.01 typically, so we use more conservative bounds on K_3m of |K_3m| < 1. In general, spin data are most sensitive to γ₀ and K_2m, which we call the ‘first-order parameters’. We call K_3m the ‘second-order parameters’.

The MCMC was determined to converge when the fractional change in autocorrelation time (computed every 100 iterations) was 1 per cent, and the number of iterations already computed was more than 100 times the autocorrelation time. The MCMC fit also was set to terminate if more than 10⁵ iterations were run, but this only occurred for fits computed for asteroid encounters in which degeneracies were present, for example, when the asteroid has rotational symmetry. About 10⁴ iterations was often sufficient.

Before the MCMC was run, local minima in the likelihood were found via the Nelder–Mead algorithm implemented in scipy (Gao & Han 2012). Generally, only one local minimum existed, except when K₂₂ = 0 in which case rotational symmetry caused all values of γ₀ to be degenerate. Ensemble walkers were initialized near the local minimum, distributed by a Gaussian approximation of the likelihood as determined via the Hessian of the likelihood at the minimum. Due to the high sensitivity of the angular velocity data to density moments, the minimization procedure sometimes failed to isolate the minimum likelihood. Therefore, a simpler simulation without the K_3m terms of equation (4) was first used to minimize likelihood as a function of the first-order parameters γ₀ and K_2m, and then the full simulation was used to find the second-order parameters K_3m with the first-order parameters fixed.

To further ensure convergence, we first minimized with respect to data truncated soon after perigee. After convergence, we refined the minimum by minimizing based on the full data, with the previous minimum as the initial estimate.

2.4 Density distribution constraints

The asteroid density distribution |$\rho _\mathcal {A}(\boldsymbol r)$| is not uniquely determined via the density moments. For example, the mass of the asteroid is unconstrained so |$\rho _\mathcal {A}(\boldsymbol r)$| cannot be determined on an absolute scale. However, by making sufficient assumptions about the density distribution, we can nevertheless measure fluctuations in |$\rho _\mathcal {A}(\boldsymbol r)$| across the asteroid from K_ℓm. To best understand the density distribution of an asteroid, an ensemble of models with differing assumptions is desired so that common traits across the models can be identified. To this end, we outline two possible models here and discuss two more in Supplementary Material C.

We assume that the asteroid’s surface is known in the inertial frame from radar data. Since the asteroid tumbles during the encounter, we also assume that the centre of mass of the asteroid is known in the inertial frame. The density moments are extracted from flyby data, but these are known instead in the body-fixed frame.

To compare the asteroid surface and extracted density moments in the same coordinate system, we define a new frame called the ‘hybrid frame’. The hybrid frame is co-located with the body-fixed frame, but its orientation is known with respect to the inertial frame; it has |$\boldsymbol \Omega _0 \parallel \boldsymbol {\hat{z}}_\text{hybrid}$| with third Euler angle γ = 0. The hybrid frame initially differs from the body-fixed frame only by a rotation around |$\boldsymbol {\hat{z}}_\text{hybrid}=\boldsymbol {\hat{z}}_\text{body-fixed}$| of γ₀. Such rotations affect density moments via

$$\begin{eqnarray*} K_{\ell m}^\mathrm{hybrid} = e^{-im\gamma _0}K_{\ell m}^\mathrm{body\text{-}fixed}. \end{eqnarray*}$$

(8)

Thus, values and uncertainties on |$K_{\ell m}^\mathrm{body-fixed}$| and γ₀ (obtained from the encounter data) can be translated into values and uncertainties on |$K_{\ell m}^\mathrm{hybrid}$|⁠. The surface model can be rotated from the inertial frame to the hybrid frame so that both the surface and K_ℓm are now in the same frame and are comparable. Henceforth, we will operate only in the hybrid frame and suppress the label.

When K_3m are extracted, 12 moments are constrained (excluding the trivial K₀₀ and K_1m, which govern the centre of mass). In cases where γ₀ is known very accurately such that the uncertainty increase imposed by equation (8) is very small, it is numerically favourable to treat ℑK₂₂ and K₂₁ as fixed, just as |$\Im K_{22}^\text{body-fixed}$| and |$K_{21}^\text{body-fixed}$| are fixed. In this case, there are nine constrained moments. Since γ₀ is precisely known for all cases studied in this paper, we will always study this case.

Since the asteroid mass cannot be determined by this analysis, it is convenient to additionally set the mass equal to its volume, so that the average density is ρ_avg = 1 and thus the extracted densities can be interpreted as ratios of ρ/ρ_avg. We then place additional constraints that 0.25 < ρ < 3 to ensure realistic densities. Using this constraint as a prior, we design another MCMC given one of the two density distribution models discussed below. The likelihood used is the multivariate-Gaussian approximation of the posterior distribution for K_ℓm, extracted from flyby data.

2.4.1 Finite element model

To define the ‘finite element’ model, we divide the asteroid into N finite elements of uniform density and use the density ρ_i of each as parameters. The four constraints imposed by the known mass and centre of mass of the asteroid (seven with K₂₁ and ℑK₂₂ are also fixed) are used to fix some of the ρ_i. These constrained densities are easily computed since the asteroid mass |$\mu _\mathcal {A}$| and the product |$I_\mathcal {A}K_{\ell m}$| are both linear functions of ρ_i. Because |$I_\mathcal {A}K_{00} = a_\mathcal {A}^2\mu _\mathcal {A}$| by definition (which is known), and the other fixed moments are zero, |$I_\mathcal {A} K_{\ell m}$| is known for all the fixed moments, and computing the corresponding constrained densities is simply a matrix inversion. The size and location of each finite element must be chosen before extracting its density.

The value of N must be carefully chosen, since it embodies a balance between the precision and accuracy of the resulting distribution. If N is set equal to the number of data points, an accurate solution is guaranteed, but uncertainties are inflated. If N is chosen lower, the choice of element layout might exclude a distribution that exactly matches K_ℓm, but uncertainties are diminished due to less redundancy in the model. For the rest of this paper, we use N = 12 finite elements, which corresponds to seven constrained elements and five degrees of freedom (d.o.f).

To arrive at this choice of five d.o.f, we extracted five density distributions with uncertainties from five random grids for the asymmetric asteroid. We calculated the mean density uncertainty over the asteroid for all five grids σ_ρ, as well as the average deviation from the true density distribution Δρ and the significance of that deviation Δρ/σ_ρ. We did this for different levels of observational precision and with either 9, 7, 5, 3, or 2 d.o.f. Our data revealed that five d.o.f appeared to show the lowest density uncertainty while still producing low significance of deviations from the true density distribution.

2.4.2 Lumpy model

A drawback of the finite element model is that the generated density distribution might not be representative of the asteroid if the elements are not optimally placed. We therefore describe an alternate model which includes the positions of the elements as parameters at the cost of resolution similarly to de Wit et al. (2012). We call this the ‘lumpy’ model.

Suppose the asteroid is formed of N constant-density, possibly overlapping ‘lumps’, enclosed within a constant-density substrate whose surface is visible to observers. The substrate mass and added mass of the lumps are denoted by μ_i, each with position |$\boldsymbol r_i$| (relative to the asteroid’s centre of mass), density moments |$K_{\ell m}^{(i)}$|⁠, and length a_i. Here, i denotes the index of the lump where the substrate is i = 0. We do not need to include I_i as a free parameter because these lumps have constant density, so that |$I_i = \mu _i a_i^2$|⁠. Furthermore, by requiring that a lump’s density moments be computed relative to the centre of mass of the lump, we have |$K_{1m}^{(i)} = 0$|⁠.

The translation rules of spherical harmonics (van Gelderen 1998) give that the asteroid density moments in the hybrid frame are

$$\begin{eqnarray*} K_{\ell m} &= &\left[ \sum _{i=0}^N \frac{a_i^2 + r_i^2}{a_\mathcal {A}^2}\mu _i \right]^{-1} \nonumber \\ &&\times \left[ \sum _{i=0}^N \sum _{\ell ^{\prime } m^{\prime }}\mu _i \frac{a_i^{\ell ^{\prime }}}{a_\mathcal {A}^\ell } R_{\ell - \ell ^{\prime }, m - m^{\prime }}(\boldsymbol r_i)K_{\ell ^{\prime } m^{\prime }}^{(i)} \right] \end{eqnarray*}$$

(9)

where the unmarked sum limits are 0 ≤ l′ ≤ l and −ℓ′ ≤ m′ ≤ ℓ′. We also have total mass and centre of mass constraints:

$$\begin{eqnarray*} \mu _\mathcal {A} = \sum _{i=0}^N \mu _i; \qquad 0 = \sum _{i=0}^N \mu _i \boldsymbol r_i. \end{eqnarray*}$$

(10)

Additional assumptions can be imposed on |$K_{\ell m}^{(i)}$|⁠. For example, we can require that the lumps be ellipsoids, so that |$K_{3 m}^{(i)} = 0$|⁠. The most extreme case is spherical lumps, which have K_ℓm = 0 for ℓ > 0. K₀₀ = 1 is also guaranteed by definition, meaning that each spherical lump has only five d.o.f (a_i, μ_i, and |$\boldsymbol r_i$|⁠). The substrate has one degree of freedom, since only μ₀ is unknown. Thus, this spherical lumpy model has 5N − 3 total d.o.f. Again, the choice of N affects the accuracy and uncertainty of the model results. For the rest of this paper, we use N = 1 spherical lump for simplicity, corresponding to two d.o.f. In Section 4.1, we also use N = 2 spherical lumps which possesses seven d.o.f.

2.4.3 Density distribution uncertainties

Once the parameters of the density distribution models have been extracted, each MCMC sample corresponds to its own density distribution. To generate an average distribution, we randomly choose 1000 of these samples and define the density distribution ρ at each point to be the mean of the densities at that point across the 1000 samples and take the standard deviation of the sample densities as the (uncorrelated) density uncertainty σ_ρ. In the case of the finite element model, we run the MCMC 20 times and select 1000 parameter samples from each to reduce dependence on the initial choice of finite element locations.

3 RESULTS

The techniques described above were implemented in a publicly available toolset called AIME (Asteroid Interior Mapping from Encounters). To demonstrate, we provide a full density distribution retrieval applied to synthetic data for two ‘reference asteroids’, which correspond to the following encounter parameter choices.

• An orbit around a spherical, non-rotating, Moonless Earth with 6 km s⁻¹ excess velocity and perigee at 5 Earth radii. This orbit was chosen to roughly match that of Apophis and corresponds to an eccentricity of 3.88.

• An initial roll of γ₀ = π/8.

• A cadence of 2 min and observational uncertainty of σ_θ = 10⁻² and σ_P/P = 10⁻⁷.

• A rotational period of 9 h, with the initial angular velocity vector distributed between the |$\boldsymbol {\hat{X}}$|⁠, |$\boldsymbol {\hat{Y}}$|⁠, and |$\boldsymbol {\hat{Z}}$| axes in a 1: 2: −2 ratio.

• An asteroid with radius |$a_\mathcal {A} = 1$| km and K_3m = 0. For K₂₂ and K₂₀, we use two standard values: one with (K₂₂, K₂₀) = (0, −0.097) and one with (0.052, −0.202). Including the third point obtained by reflection K₂₂ → −K₂₂ (a 90° rotation), these are the three points that minimize the mean distance between an arbitrary point in the allowed parameter space (Supplementary Material equation A.5) and these reference values. The first point is called the symmetric case because the corresponding uniform-density-ellipsoid model is rotationally symmetric around |$\boldsymbol {\hat{z}}$|⁠. The second case is called the asymmetric case. The asymmetric case has a = 1140 m, b = 1839 m, and c = 565 m, while the symmetric case has a = b = 1411 m and c = 1008 m.

In the following sections, we first introduce the retrieval capabilities regarding density moments, then we turn to an example of a retrieval of the density distribution from the moments. For both stages of information retrieval, we find that the results are consistent with the synthetic data and the true asteroid interior.

3.1 Density moment retrieval

Fig. 2 shows our synthetic spin data for the asymmetric reference asteroid. The best-fitting model is overlaid in the top panel and residuals are shown the bottom panel. Uncertainties are plotted on the residuals with correlations between the vector components. The fit results are consistent with the data. This figure also reveals which data points are most informative. The at-perigee data are irregular and reveal information about the density moments, and the post-perigee data show torque-free tumbling behaviour which constrains K_2m via the MOI ratios as in Moskovitz et al. (2020). The post-encounter periods and phase also indirectly sheds light on the at-perigee data by constraining the rotational velocity the asteroid must have had when leaving the perigee.

Figure 2.

Data, best-fitting results, and residuals for a fit to synthetic data simulated for the asymmetric reference asteroid. Uncertainty bands are also shown. The best-fitting results are consistent with the data.

Fig. 3 shows a corner plot of the PPDs of the ten parameters (namely γ₀ and K_ℓm for ℓ ≤ 3), marginalized to functions of one (histograms) or two (contours) variables. The true parameters are shown and usually lie within 1σ or 2σ of the ΔK_ℓm = 0, where ΔK_ℓm is the difference between the mean posterior K_ℓm and the true K_ℓm. The PPDs are generally Gaussian and sometimes show strong correlation between parameters, but no continuous degeneracy occurs. We performed 48 independent minimizations of the likelihood before the MCMC fit began, each with an initial point chosen randomly in the parameter space. All converged to the same minimum, demonstrating that the model lacks discrete degeneracy as well.

Figure 3.

PPDs extracted from synthetic encounter data for the asymmetric reference asteroid. Samples from the MCMC fit are shown as individual points, and the contours enclose 1σ, 2σ, and 3σ confidence regions. True values are shown as blue lines. PPDs are Gaussian and show no degeneracies.

3.2 Density distribution retrieval

To provide an example of a density distribution extraction, we consider an asteroid with a core and ask whether AIME can resolve the location, mass, and size of the core. Specifically, we use a spherical core of radius 300 m, placed 500 m from the centre and with density 50 per cent greater than that of the surrounding asteroid. We use the lumpy model to extract a distribution in this section because it is designed to look for cores.

Synthetic data were generated for the new asteroid and moments were extracted via the process described above. The true density moments were within the confidence intervals of the extracted moments. A density distribution was then extracted and shown in Fig. 4 next to the true distribution. Visually, the extracted distribution is virtually identical to the true distribution, which is emphasized by the low uncertainties of |$\langle \sigma _\rho / \rho \rangle = \mathcal {O}(10^{-3})$| throughout the asteroid. Again, the extracted density distribution is consistent with the truth despite the low density uncertainties. Both the density moments and the extracted density distribution are inconsistent with a uniform distribution; i.e. the deviation from uniform shown in Fig. 4 is statistically significant.

Figure 4.

Cross-sectional slices of the extracted distribution (left) and true distribution (right) of the asymmetric reference asteroid. The density distribution was created via the lumpy model. The difference between the extracted and the true distributions and density uncertainty are also shown (bottom left / right). The extracted distribution is nearly identical to the true distribution and is statistically consistent with the truth. These figures are available in animated form in the Supplementary Material.

4 DISCUSSION

In Section 4.1, we discuss the differences between the two density distribution models for many asteroids and find that the distributions and uncertainties they generate are strongly model-dependent.

Given this understanding of our density extraction models, we answer the important questions of which encounters can be successfully studied using AIME, and how observational campaigns can be designed to best take advantage of an encounter. To quantify AIME’s success, we use the median of a density uncertainty distribution 〈σ _ρ/ρ〉 in Section 4.2. But since this median is so model-dependent, we also analyse the success of the method via the density moment uncertainty σ(K_ℓm) (which is not moment dependent) in Section 4.3. In Section 4.4, we study AIME’s dependence on the properties of the central body.

4.1 Density distribution model caveats

To study the output of the density distribution models, we begin by executing them on synthetic data generated for the two reference asteroids. The resulting density distributions and uncertainties are shown in Fig. 5. These distributions are consistent with the moments extracted by the MCMC in Section 2.3 – indeed all the distributions we show have moments consistent with the data, even for non-ellipsoidal asteroids (not shown).

Figure 5.

Cross-sectional slices of the density distributions extracted via the finite element model for the asymmetric (top two rows) and symmetric (bottom two rows) uniform-density reference asteroids. The finite element model (first and third rows) and the lumpy model (second and fourth rows) are employed. From left to right, the densities (divided by the average density), deviations from the true density, uncertainties, and significance of the deviations are plotted. These figures are available in animated form in the Supplementary Material. Extracted densities are generally within 10 per cent of the truth.

For the finite element model with the reference observational setup, the uncertainty on observations is such that the density distribution is generally within 10 per cent of the true uniform density (second column) while the density uncertainty is generally less than 40–50 per cent of the density value at any point in the asteroid (third column). In no place is the significance of these deviations from the true distribution greater than 1σ (last column). The lumpy model yields distributions with much lower uncertainty than the finite element distributions (maximum uncertainty on the order of 1 per cent of the local density or less) due to its few degrees of freedom and its particular design; with the one-lump-model and with asteroids whose surface’s centroid is the asteroid’s centre of mass (such as these uniform asteroids), the single lump must lie at the centroid with mass close to zero and with unconstrained radius. The uncertainty of regions far from the asteroid centre, where the lump is unlikely to extend, is typically very small, while regions close to the centre are more likely to be contained inside a lump and hence have greater density uncertainty. This uncertainty is entirely model-driven and can be discarded.

We also explore model behaviour in non-uniform density asteroids. First, we consider again the off-centre core asteroid used in Section 3.2; distributions are shown in Fig. 6. As remarked before, the lumpy model fit is successful. However, the finite element model does not reproduce the core. The deviation from the true density extends to as much as 17 per cent in some locations, leading to a maximum significance of 0.52σ. Visually the finite element model has spread the high-density core into the rest of the asteroid.

Figure 6.

Cross-sectional slices of the density distributions extracted via the finite-element (top) and lumpy (bottom) models for an asteroid with an off-centre core. From left to right, the densities, deviations from the true density, uncertainties, and significance of the deviations are plotted. These figures are available in animated form in the Supplementary Material. The lumpy model successfully extracts the core.

This ‘spreading out’ of the density distribution is not necessarily a failure of the finite element model; the last column of Fig. 6 shows that the deviation from the true density distribution does not have high significance. In effect, the finite element model is acknowledging many possible distributions that could all have the same moments as the true distribution, while the lumpy model picks one. In our case, the lumpy model was correct, but for asteroids without a discrete, spherical core like this, it may not be.

To further highlight the model-dependence of the extracted density distributions, we consider two final asteroids which lead to inflated density distribution uncertainties independent of the data quality. Specifically, we will highlight the lumpy model’s degeneracy for a centred core and the two-core lumpy model.

Fig. 7 shows density distributions extracted via the finite element model and the lumpy model for a centred core of density three times the surrounding density. Results are similar to the off-centre core shown in Fig. 6 in that the finite element model does not isolate the lump, instead spreading the excess mass over the asteroid. Unlike the off-centre core example, the lumpy model is not able to recover the true distribution either. It produces deviations from the true density distribution of roughly the same size as the finite element model, with large significance.

Figure 7.

Cross-sectional slices of the density distributions extracted via the finite-element (top) and lumpy (bottom) models for an asteroid with a centred core. From left to right, the densities, deviations from the true density, uncertainties, and significance of the deviations are plotted. These figures are available in animated form in the Supplementary Material. The resulting density distribution is consistent with the density moments but does not represent the true distribution due to degeneracy.

The success of the model in the off-centre case was due to the fact that the shape of the asteroid was offset from the centre of mass by a corresponding amount, assumed to be known precisely. The position of the lump was therefore observed up to one free parameter: the lump’s mass. In the centred core case, the core mass does not affect the asteroid centre of mass so the mass is unconstrained and uncertain. The underlying assumption that the asteroid’s centre of mass is so precisely known stems from the fact that the shape of the asteroid is observed to rotate around its centre of mass. If observations do not allow the centre of mass to be determined in this way, then the lump’s position will be more uncertain for the off-centre case.

We also consider an asteroid with two lumps of radii 300 m and density three times the surrounding density. Each lump is located 500 m from the centre of the asteroid, so that they counterbalance and the asteroid’s observed centre of mass is its surface’s centroid. The corresponding two-lump lumpy model has seven d.o.f, in contrast to the five d.o.f of the finite element model or the two d.o.f of the one-lump lumpy model. Both models are run on this asteroid and the resulting distributions are shown in Fig. 8.

Figure 8.

Cross-sectional slices of the density distributions extracted via the finite-element (top) and the two-lump lumpy (bottom) models for an asteroid with two counterbalancing cores. From left to right, the densities, deviations from the true density, uncertainties, and significance of the deviations are plotted. These figures are available in animated form in the Supplementary Material. The additional lump greatly increases uncertainty, but the resulting distribution is close to accurate.

Fig. 8 again shows that the finite element model is unable to isolate the lumps except to predict generally increased density near the centre. On the other hand, the lumpy model detects that the two lumps are opposite each other and of roughly the same radius and mass. These radius and mass values are also close to the true values. The model places the lumps correctly near the xy-plane of the asteroid but does not perfectly align them with the y-axis, resulting in high Δρ where the true lumps and predicted lumps do not intersect. Uncertainty is also high, due to the extra d.o.f. The difference between the predicted distribution and the true distribution is due entirely to the uncertainty in the moments computed by the fit to flyby data. The density distribution model reproduces those moments perfectly, misplacing the lumps by a few hundred metres because that offset is determined only by the K_3m density moments, which are not as well constrained as K_2m.

4.2 Probing density sensitivity to encounter properties

Here, we investigate how density uncertainty 〈σ_ρ/ρ〉, extracted by both the finite element and lumpy models, depends on physical asteroid properties (the encounter’s orbit, the asteroid’s true shape, and its initial rotational velocity) and observational properties (the data uncertainty, the cadence of observations, and gaps in data coverage).

Fig. 9 depicts 〈σ _ρ/ρ〉 as a function of these other encounter properties, except for the dependence of 〈σ _ρ/ρ〉 on initial spin direction which is better shown as a map in Fig. 10. We immediately see that the two density distribution models produce very different density uncertainties, due to the strong model dependence of their output. We focus on the lumpy model’s output because it is more sensitive to the encounter properties than the finite element model.

Figure 9.

Density uncertainty as a function of physical / observational encounter properties (top / bottom). Both the lumpy and finite element models are used, though the finite element model uncertainty is model-dominated and not very sensitive to encounter properties. The uncertainty threshold (red line) is shown and the properties that exceed this threshold are shaded. The vertical black dotted lines are the property values of the reference asteroid. Perigee, rotational period, and observational uncertainty most strongly affect density uncertainty.

Figure 10.

Density uncertainty computed via the lumpy model as a function of initial spin pole direction. The weak uncertainty cut-off (red line) is shown, as is the reference spin pole direction (orange star). The strong cut-off is never exceeded and does not appear. Initial angular velocities perpendicular to the orbital plane lead to greatest uncertainty.

Fig. 9 reveals that the most limiting physical properties of the asteroid are its perigee and period. Excess velocity does not greatly affect the uncertainty of the density distribution, and asteroid length is only important for |$a_\mathcal {A}\sim$| tens of metres. By contrast, the radius of Apophis is ∼300 m and its perigee and rotational period also obey the strong constraints (Giorgini et al. 2008). The comparison to Apophis is complicated by the fact that Apophis is smaller than our |$a_\mathcal {A}$| value, is tumbling (Pravec et al. 2014), and may change slightly in physical properties due to tidal interaction during the encounter (Yu et al. 2014; Hirabayashi et al. 2021). Further work must therefore be done to apply this analysis to Apophis.

In addition to the period constraint shown in this table, Fig. 10 shows an uncertainty dependence on initial spin direction. The poles (i.e. angular velocities perpendicular to the orbital plane) lead to greater uncertainty, potentially by as much as a factor of 10 over to the orbital plane. Compared to other physical parameters, initial spin direction is not very constraining.

Fig. 9 also demonstrates the strong affect of properties of the observational campaign on the density uncertainty. Most vital are the uncertainties on the asteroid’s instantaneous rotational period σ_P/P and direction σ_θ, which require precision of the order of tens to hundreds of milliseconds and degrees every cadence, respectively. This could be accomplished by multiple precise angular velocity measurements from multiple telescopes or by increasing the time between observations to maximize the change in period between observations. Acknowledging and correcting for correlations in uncertainty between data points could also reduce uncertainty without requiring such high period uncertainty, as would increasing the data set size to include more post-flyby tumbling data.

On the other hand, the cadence of observations Δt and the presence of gaps in the data T_gap do not affect results as strongly. Short cadence appears preferable, though as much as 20–40 min between observations is still not as detrimental to data quality as other properties. Likewise, gaps appearing with size T_gap, even hours in length, do not greatly increase density uncertainty.

To summarize the information contained in Fig. 9, we define a ‘threshold’ on each encounter property beyond which we say the data quality is too low to extract meaningful information concerning the density distribution of the asteroid. The value of this threshold is set at the point when the lumpy model produces |$\langle \sigma _\rho / \rho \rangle = 0.1{{\ \rm per\ cent}}$|⁠, since this is just under the maximum uncertainty that the lumpy model produces. This threshold is marked as a horizontal red line in Fig. 9, and the values at which the encounter parameters exceed the threshold are shown in Table 1.

Table 1.

Thresholds on physical / observational properties (top / bottom) necessary to obtain density distributions with useable uncertainty. Perigee and observational uncertainty are the most constraining properties.

Encounter property	Threshold
Perigee (rp)	<18 Earth radii
Spin period uncertainty (σP)	<0.27 s
Spin pole uncertainty (σθ)	<35°

Table 1 is consistent with the conclusions we drew for Fig. 9. The physical properties of Apophis’ 2029 encounter are consistent with the thresholds, with the perigee and observational uncertainties posing the greatest challenge to successful density distribution extraction. The values of these thresholds are interdependent and may change if the encounter conditions are changed, so they should not be taken to be precise values.

4.3 Moment sensitivity to encounter properties

To address the dependence of 〈σ_ρ/ρ〉 on the choice of model, we also investigate the model-independent moment uncertainty σ(K_ℓm) as well, defined as the range of K_ℓm values that contains 68.27 per cent of the marginal PPD. Since there is no degeneracy between moments and the actual encounter data, σ(K_ℓm) is well-defined and less noise-prone than 〈σ_ρ/ρ〉, though its physical relevance is not as obvious.

Figs 11 and 12 display moment uncertainty as a function of physical and observational encounter properties, respectively. Fig. 13 additionally depicts moment uncertainty as a function of initial spin pole. The thresholds of table 1 are depicted as red lines. Each panel of each figure is generated by creating about 50 synthetic data sets with different encounter parameters and extracting the density moments from each.

Figure 11.

1σ and 2σ confidence intervals for the first-order parameter PPDs (top) and second-order parameters (bottom) as a function of (left to right) perigee, excess velocity, asteroid length, and rotational period. The vertical dotted line indicates the reference asteroid values. The red line indicates the uncertainty threshold.

Figure 12.

1σ and 2σ confidence intervals for the first-order parameter PPDs (top) and second-order parameters (bottom) as a function of (left to right) period and spin pole uncertainty, observational cadence, and length of gaps in the data. The vertical dotted line indicates the reference asteroid values. The red line indicates the uncertainty threshold.

Figure 13.

1σ uncertainties for the first-order parameters (top) and second-order (bottom) as a function of the initial direction of spin in the inertial frame. All maps are made in the Mollweide projection. The orange star indicates the reference spin pole. Green dots are the sampled spin pole directions. The red contours enclose regions above the strong uncertainty threshold on 〈σ_ρ/ρ〉. The weak threshold is never exceeded. Beyond the |$\pm \boldsymbol {\hat{Z}}$| increase in uncertainty, there is also increased moment uncertainty for |$\pm \boldsymbol {\hat{Y}}$|⁠.

Fig. 11 reveals that 〈σ _ρ/ρ〉 from Fig. 9 is more sensitive to σ(K_2m) than to σ(K_3m). For instance, σ(K_2m) is constant as |$a_\mathcal {A}$| is varied despite a dramatic increase in σ(K_3m) for low |$a_\mathcal {A}$|⁠. The resulting 〈σ _ρ/ρ〉 is mostly constant. The opposite is true for P_ω, where σ(K_3m) (except for m = 0) are mostly constant, and 〈σ _ρ/ρ〉 follows the trend of σ( K_2m) and increases for low rotational period. A consequence is that if more post-flyby tumbling data are collected, placing stronger constraints on K_2m rather than K_3m, then K_2m at some point might have essentially no uncertainty. In this case, uncertainty on K_3m will be dominant and the most constraining parameters will change. Rotational period P_ω will cease to be a vital parameter but asteroid length |$a_\mathcal {A}$| will be because σ(K_3m) are much more dependent on |$a_\mathcal {A}$| than P_ω.

Since Figs 11 and 12 show such clear uncertainty dependence on encounter properties, we discuss the implication of each panel individually in the following sections.

4.3.1 Orbital elements

A Keplerian orbit is completely described by five parameters, but three describe the orbit’s orientation with respect to the central body. They are therefore redundant with the orientation of the inertial frame and we do not investigate them here. We parametrize the remaining two parameters by the perigee distance r_p and excess velocity v_∞.

Fig. 11 shows that moment uncertainty depends so strongly on perigee that for r_p > 10 Earth radii, σ(K_3m) is constrained by the prior boundaries of ±1 for m ≤ 2. Fig. 9 also demonstrates that low perigee yield more certain density distributions as extracted by both the lumpy and the finite element models. Such a strong strong dependence is expected from equation (4); it is caused by the |$D^{-\ell ^{\prime }}$| factor contained in |$S_{\ell ^{\prime } m}(\boldsymbol D)$|⁠.

By contrast, density moment uncertainty shows only a slight increase with v_∞. This is likely due to the fact that larger v_∞ leads to a faster and flatter orbit with less time spent close to the planet, where tidal torque is strongest. The change in encounter timing also adjusts the orientation of the asteroid at perigee, which has a separate effect on moment uncertainty. We correct this undesired orientation dependence by adjusting γ₀ so that γ is roughly the same value at perigee (where torque is highest) for all points in the data set.

4.3.2 Initial angular velocity

In Fig. 11, we show σ(K_ℓm) as a function of initial rotational period P_ω. As in the above section, the value of γ₀ was corrected to ensure roughly constant orientation at perigee. K₂₀ and K₂₂ show very large uncertainty for P_ω ≲ 4 h because these fast rotators tumble very little after perigee. This increases uncertainty on the K_2m parameters, which are largely constrained by post-encounter tumbling.

We expect that fast rotators would not tumble post-encounter for the following reason, For small P_ω, all the dynamical variables vary much more slowly than the orientation γ. Approximating each variable as constant over one full rotation of γ, the integral of the first-order contribution of |$\boldsymbol \tau$| over γ ∈ (0, 2π) gives no secular first-order torque to force the asteroid to tumble. However, this effect does not apply to the second-order parameters, since the integral over the second-order term of |$\boldsymbol \tau$| does not vanish, as seen in the figure. An asteroid with large K_3m moments might therefore yield better uncertainties at these low rotational periods (e.g. a non-uniform or non-elliptical asteroid).

The tidal torque experienced by the asteroid is also affected by the initial direction of asteroid spin |$\boldsymbol \Omega _0$| both because spin sets the initial asteroid orientation up to γ₀ and because of the spin-dependence of the rotational equations of motion (Supplementary Material equation A.14).

Fig. 13 shows increased moment uncertainty for initial spin pole |$\boldsymbol \Omega _0 \parallel \boldsymbol {\hat{Z}}$|⁠, just as Fig. 10 shows increased density distribution uncertainty. Many other moments also exhibit increased uncertainty for |$\boldsymbol \Omega _0 \parallel \boldsymbol {\hat{Y}}$|⁠. This pattern is explained by the tidal torque equation (equation 4). By plugging in values for the Euler angles, |$\boldsymbol z \parallel \boldsymbol {\hat{Z}}$| and |$\boldsymbol z \parallel \boldsymbol {\hat{Y}}$| at perigee lead to |$\boldsymbol \tau \propto K_{22} \boldsymbol {\hat{z}}$| to first-order, and |$\boldsymbol \tau \parallel \boldsymbol {\hat{X}}$| at perigee leads to |$\boldsymbol \tau = 0$| to first-order. |$\boldsymbol \tau \parallel \boldsymbol {\hat{z}}$| implies that only the period of the asteroid is changed, not its spin pole direction, implying that the asteroid tumbles less after the flyby. As discussed above, tumbling allows precise constraints on K_2m, so that reduced tumbling results in greater uncertainty.

Since |$\Omega _0 \parallel \boldsymbol {\hat{z}}$| at the start of the simulation, high uncertainty for |$\boldsymbol {\hat{z}} \parallel \boldsymbol {\hat{Y}}$| at perigee means high uncertainty for |$\boldsymbol \Omega _0 \parallel \boldsymbol {\hat{Y}}$| assuming that the perigee torque effects are dominant. Since |$\boldsymbol \tau = 0$| at perigee for |$\boldsymbol {\hat{z}} \parallel \boldsymbol {\hat{X}}$|⁠, the encounter may be dominated by non-perigee effects in that case. This may explain the increased uncertainty for |$\boldsymbol \Omega _0 \parallel \boldsymbol {\hat{Z}}$| and |$\boldsymbol {\hat{Y}}$| but not |$\boldsymbol {\hat{X}}$|⁠.

4.3.3 Observational uncertainty

Two parameters, σ_θ and σ_P, govern the observational uncertainty of the data set. These parameters are defined in Section 2.2; σ_θ represents the standard deviation of the angle between the true spin pole and the observed spin pole, while σ_P represents the standard deviation of the rotational period. Rather than explore the full space spanned by these two values, we fix one and allow the other to vary to better assess whether uncertainty in spin pole or uncertainty in period more strongly affects uncertainty. This dependence is displayed in Fig. 12. Moment uncertainty σ(K_ℓm) grows linearly with observational uncertainty (σ_θ or σ_P).

In particular, we might ask if some error |$\boldsymbol {\delta }_\omega$| is added to angular velocity |$\boldsymbol \omega$|⁠, does it affect results more strongly if it is parallel to |$\boldsymbol \omega$| (affects the period) or perpendicular (affects the spin pole direction)?

Let |$\delta = |\boldsymbol {\delta }_\omega | / |\boldsymbol \omega |$| and δ ≪ 1. Then if |$\boldsymbol {\delta }_\omega \parallel \boldsymbol \omega$|⁠, it decreases the period P_ω by P_ωδ. This is a fractional change in period of δ. If |$\boldsymbol {\delta }_\omega \perp \boldsymbol \omega$|⁠, then the spin pole angle changes by δ radians. Comparing the σ_P/P (fractional change in period) and σ_θ (spin pole angle) columns of Fig. 12, one can see that a given value of σ_θ contributes a much smaller moment uncertainty than the same value of σ_P/P. This is also visible in Fig. 9 for the lumpy model. In other words, if |$\boldsymbol {\delta }_\omega \perp \boldsymbol \omega$| using the symbols defined above, then |$|\boldsymbol {\delta }_\omega |$| can be large. But if |$\boldsymbol {\delta }_\omega \parallel \boldsymbol \omega$|⁠, then |$|\boldsymbol {\delta }_\omega |$| must be very small. Period precision is therefore more vital than spin pole direction precision when it comes to decreasing uncertainties.

4.3.4 Cadence and data gaps

The time between observations of asteroid angular velocity, (cadence, Δt), may vary depending on the observational schedule of the observing telescopes and the path of the asteroid through the sky. We measure how the moment uncertainty σ(K_ℓm) varies with cadence ranging from 2 min to 1 h in Fig. 12.

Fig. 12 displays little dependence of uncertainty on cadence Δt for Δt ≲ 40 min. We also see flaring of uncertainty for very large cadence, largely driven by the paucity of data points. However, uncertainty dramatically increases for many parameters at about Δt = 30–40 min, a time-scale which depends both on the asteroid rotational period P_ω and the time-scales of its orbit.

Fig. 12 shows that as long as Δt is less than this threshold, the influence of cadence on σ is small, but shorter cadence leads to lower uncertainties.

In certain circumstances, spin data might not be able to be captured for a close encounter at perigee. The asteroid might dip below the horizon, or it might pass too close to the sun to be observed. The resulting gap in data is intended to be captured by the T_gap parameter of Fig. 12, which deserves to be more fully defined.

We mask the perigee of the counter by removing a duration T_gap of data centred on the perigee, where T_gap ranges from 0 to 3 h. To prevent lack of precision induced by lower amounts of data when T_gap is large, we always cut 3 h−T_gap from the data set, half from the beginning and half from the end, so that each data set produced for all T_gap has the same size. We cut around the perigee because tidal torque is the greatest at perigee, and we expect that part of the data set to be most valuable. Indeed, Fig. 12 shows that K_3m especially are more uncertain for T_gap ≳ 1.5 h. However, Fig. 9 shows that T_gap never increases density uncertainty above the threshold, indicating that large amounts of data can be cut without compromising AIME. As with the threshold for cadence discussed above, this 1.5 h cut-off may depend on the asteroid rotational period or the orbital time-scales. It likely also depends on the observational cadence used.

4.3.5 Other parameters

We also study the dependence of moment uncertainty on the asteroid MOI and on the asteroid’s initial orientation, but these relationships are simple enough that they are not included in Figs 11 and 12. Moment uncertainty is generally unrelated to the MOI ratios, except when the asteroid is rotationally symmetric (e.g. for the symmetric reference asteroid). In this case, the initial orientation of the asteroid γ₀ is undefined, creating degeneracy and inflating density moment uncertainty. For symmetric or near-symmetric asteroids, this issue could be resolved by re-parametrizing the MCMC to remove this degeneracy.

Moment uncertainty does depend on the asteroid length, which sets the MOI itself, as shown in Fig. 11. This dependence is relatively simple; moments are damped by factors of |$a_{\mathcal {A}} / D$|⁠, so that reducing asteroid size without reducing perigee yields poor constraints on K_3m

Moment uncertainty is also affected by γ₀. We measured the moment uncertainty as a function of γ₀ for the asymmetric reference asteroid, keeping all other parameters constant. Moment uncertainties generally varied by factors of two or less. The details of this dependence are strongly dependent on the initial spin pole and the asteroid shape, so the data are not shown.

4.4 Sensitivity to central body properties

In all the above studies, we assumed a spherical planet (J_ℓm = 0 for ℓ ≥ 1). J_1m = 0 is enforced by the coordinate definitions, so the effect of central body non-sphericity is limited to the J_2m terms and damped by a factor of |$(a_\mathcal {B} / D)^2$|⁠. We expect these parameters to have small effect on the asteroid behaviour.

Here, we define oblateness as |$\epsilon = (I_z - I_x)/(\mu _\mathcal {B} R_\mathcal {B}^2)$|⁠, where I_x, y, z are the central body moments of inertia along the principal axes, and I_x = I_y. |$R_\mathcal {B}$| is the true radius of the body (not |$a_\mathcal {B}$| from equation 3). For an equatorial orbit, ϵ and the central body density moments (Supplementary Material equation A.3) are related by J_ℓm as ϵ = −10J₂₀/3 and J₂₂ = 0. Since an oblate ellipsoid is mirror-symmetric around all three axes, J_3m are all zero. The next order of tidal torque is therefore J_4m, damped by an additional |$(a_\mathcal {B}/D)^2$| factor, and non-ellipsoid corrections to the central body shape. We do not consider these extra terms.

Fig. 14 displays moment uncertainty σ(K_2m) of the first-order parameters as a function of ϵ across a reasonable range of central body oblatenesses based on those of Solar system planets (Pater & Lissauer 2015). It also shows linear best-fitting curves for moment uncertainty as a function of oblateness. All other parameters, include |$I_\mathcal {B}$| which parametrizes the central body radius, are kept constant. Almost no dependence of σ(K_ℓm) on oblateness ϵ is apparent, although moment uncertainty does measurably decrease for oblate central bodies.

Figure 14.

Top: 1σ and 2σ confidence intervals for the first-order parameter PPDs as a function of oblateness ϵ. All other parameters, including the central body radius, are kept constant. Linear best-fitting lines to σ(K_2m) (black, dotted) are plotted. Bottom: The difference between PPD means extracted from a zero-oblateness model and the true parameters given data with true oblateness ϵ_true ≠ 0. Also shown in both figures are the oblatenesses of reference Solar system bodies. Moment uncertainty depends little on oblateness, but the best-fitting parameter estimates are affected enough by oblateness that oblateness must still be modelled.

Given the small effect of ϵ on K_ℓm, it might be tempting to neglect the planetary oblateness when fitting K_ℓm to data. However, the bottom panel of Fig. 14 demonstrates that doing so is invalid. This figure displays K_ℓm as extracted by a fit assuming ϵ = 0, but run on data generated with non-zero ϵ. The difference between the PPD means and true parameters are shown. Moment uncertainties are also shown as bands. The figure shows that even for low (Earth-scale) oblateness, the fit results are inconsistent with the true K_ℓm values, since ΔK_ℓm = 0 is not contained in the 2σ band. This effect is much worse for large oblateness, growing to a difference on the order of |$\mathcal {O}(100)\sigma$| for Jupiter’s oblateness. Therefore, accurately modelling central body oblateness to high precision is essential for the accurate estimation of fit parameters. For non-equatorial orbits, with J₂₂ ≠ 0, we also expect J₂₂ to affect the accuracy of the fit results to a similar degree, with the additional requirement of using the correct asteroid orbital plane.

J₂₀ has a slightly more general definition than oblateness. If the planet has a moon, the integral defining J₂₀ (equation 1) can be extended to include this extra mass, though this can only be done when the asteroid never passes inside the moon’s orbit. As an order-of-magnitude estimate for this effect, two spherical objects with masses and radii of Earth and the Moon, separated by one Lunar distance, and both lying in the orbital plane has a combined oblateness of ϵ = 0.82. Extrapolating moment uncertainties via the slopes of the best-fitting lines given earlier yields a reduction in σ(K_2m) by about 25 per cent. Furthermore, J₂₂ is non-zero for this case, which likely decreases moment uncertainty even more.

This analysis suggests that large moons such as ours can improve fit quality, but further study of this effect (e.g. investigating an encounter that approaches both the Earth and the Moon closely) is beyond the scope of this paper.

Aside from oblateness, central body mass may also affect the success of AIME. To address this possibility, we run our reference asteroid through a Jupiter encounter to analyse the differences in moment uncertainty σ(K_ℓm) compared to an Earth encounter. The physical parameters of the asteroid body are kept the same as the Earth encounter case, as are the observational uncertainty and cadence. The orbit is adjusted for the Jupiter case by setting a perijove distance of r_p = 5 Jupiter radii (compared to perigee radius r_p = 5 Earth radii for the Earth encounter). The excess velocity does not strongly affect σ(K_ℓm) as shown in Fig. 11, so we keep it at the reference value.

The ratio between the moment uncertainties in the Jupiter and the Earth encounters are shown in Table 2. In all cases, the Jupiter-encounter moments are more uncertain than Earth-encounter moments. These uncertainty ratios can be understood as follows. The leading order of tidal torque is proportional to |$\mu _\mathcal {B} / D^3$|⁠. If |$D/a_\mathcal {B}$| (the ratio of the encounter distance to the central body radius) is roughly constant as in this case, then |$\mu _\mathcal {B} / D^3 \propto \mu _\mathcal {B} / a_\mathcal {B}^3 \propto \rho _\mathcal {B}$| where |$\rho _\mathcal {B}$| is the density of the central body. Therefore, tidal torque is not increased around massive planets when we keep |$r_p \propto a_\mathcal {B}$|⁠. The second-order terms are damped by an additional factor of |$a_\mathcal {A}/D$|⁠, which decreases if a massive central body is used. Since Jupiter is about 10 times larger in radius than Earth, we expect that the K_3m terms are about ten times more uncertain than the K_2m components, which is the case. Furthermore, since the orbit size is increased without a decrease in asteroid rotational velocity, the asteroid tends to tumble less for the same reasons as described in Section 4.3.2. This also increases moment uncertainty.

Table 2.

Ratio of moment uncertainty for all density moments K_ℓm between an Earth encounter and a Jupiter encounter with identical properties except for an increased perigee. Observational uncertainty and cadence are assumed to be equivalent for the Jupiter and Earth encounters. Without taking the frequencies of close encounters into account, massive planets such as Jupiter yield less precise density moment estimates.

Kℓm	σ(Kℓm)Jupiter/σ(Kℓm)Earth
γ0	1.6
K22	2.3
K20	11
ℜK33	18
ℑK33	18
ℜK32	18
ℑK32	18
ℜK31	25
ℑK31	10
K30	53

There are additional effects of central body mass which are not captured in this analysis. For example, encounters with massive planets are more plentiful, so that observation for a fixed period of time will lead to a larger number of observed encounters conducive to low-uncertainty moment extraction (small r_p, large |$a_\mathcal {A}$|⁠, etc.). This can be seen via the following equation for the area of the keyhole through which the asteroid must fly to have a perijove r_p or lower:

$$\begin{eqnarray*} A = 2 \pi r_p^2 \left(1+2\frac{G\mu _\mathcal {B}}{r_p v_\infty ^2} \right). \end{eqnarray*}$$

(11)

It is true that σ(K_ℓm) decreases with central body radius, but A increases so fast that the number of encounters that meet the uncertainty threshold will still grow, assuming that the flux of asteroids through keyholes of equal area is the same for Earth and Jupiter. Other effects, such as a change in the physical properties of the encountering asteroids, changes in asteroid rarity, and decreased observational uncertainty due to the distance between Jupiter and Earth-based telescopes, may also affect the fit uncertainties. These effects contradict, and which dominates depends on the asteroid population near Jupiter and the observation method.

5 CONCLUSIONS

We develop and demonstrate a methodology, AIME, which constrains density fluctuations in an asteroid from angular velocity changes occurring during a close encounter. We find that this inversion process is most sensitive to asteroid perigee and period, specifically requiring that the encounter perigee be ≲18 Earth radii, though this threshold depends on a number of properties of the asteroid and the observational campaign. Asteroids that tumble strongly after the encounter are also better constrained by AIME. Highly precise data, especially in the instantaneous rotational asteroid period, are also required in order to extract precise constraints. Nevertheless, for the reference asteroid and observational campaign used in this paper, we were able to extract large-scale density non-uniformities accurately and precisely, achieving density uncertainties |$\sim 0.1{{\ \rm per\ cent}}$| of the density and excluding uniform distributions for non-uniform asteroids.

We also find that the density distributions inferred are model-dependent, and their uncertainties can be dominated by model-driven uncertainties when the degrees of freedom available to the model exceeds the number of density moments that can be precisely extracted from the data. Models which are specialized to fit for certain features (such as the lumpy model) produce less uncertain results than generalized models (such as the finite element model). To efficiently use encounter data, it is therefore important to investigate multiple models and compare results.

SUPPORTING INFORMATION

Figure 4. Cross-sectional slices of the extracted distribution (left) and true distribution (right) of the asymmetric reference asteroid.

Figure 5. Cross-sectional slices of the density distributions extracted via the finite element model for the asymmetric (top two rows) and symmetric (bottom two rows) uniform-density reference asteroids.

Figure 6. Cross-sectional slices of the density distributions extracted via the finite-element (top) and lumpy (bottom) models for an asteroid with an off-centre core.

Figure 7. Cross-sectional slices of the density distributions extracted via the finite-element (top) and lumpy (bottom) models for an asteroid with a centred core.

Figure 8. Cross-sectional slices of the density distributions extracted via the finite-element (top) and the two-lump lumpy (bottom) models for an asteroid with two counterbalancing cores.

Appendix A. Tidal torque and equations of motion

Appendix B. Comparing orientation and angular velocity data

Appendix C. Additional density distribution models

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We warmly thank Emmanuel Jehin, Maxime Devogele, and Marin Ferrais for meeting with the authors to discuss this work and pointing out possible future initiatives. We also thank the anonymous reviewer for their careful reading of our paper and their comments. JTD thanks the Massachusetts Institute of Technology Undergraduate Research Opportunities Program (MIT UROP) office for funding his work. This paper made substantial use of MIT Supercloud’s facilities.

DATA AVAILABILITY

All the code used in this paper – namely, AIME – are available on GitHub.¹ The data showcased here was generated from that code. Please contact JTD with questions.

Footnotes

1

https://github.com/disruptiveplanets/AIME