DETECTION OF SEMIMAJOR AXIS DRIFTS IN 54 NEAR-EARTH ASTEROIDS: NEW MEASUREMENTS OF THE YARKOVSKY EFFECT

The Yarkovsky drift manifests itself primarily as a change in mean anomaly (or along-track position), and some observational circumstances are poorly suited to detect such changes. Examples include optical astrometry secured when the line of sight is roughly parallel to the asteroid velocity vector or when the object is at large distances from Earth. In both instances, the differences in astrometric positions can be much smaller than observational uncertainties, resulting in low sensitivity to the Yarkovsky effect. The overall Yarkovsky sensitivity depends on the orbital geometry of the NEA and on the entire set of available observations. This can be quantified rigorously. For each epoch t_i at which optical observations were obtained (1 ⩽ i ⩽ N), we predict the position P⁰_i for the best-fit orbit (da/dt = 0) and the position P*_i for the same orbit modified by a nominal non-zero da/dt. The value of the nominal rate is not important as long as it results in detectable (∼arcsecond) changes in coordinates and as long as it is applied consistently to all objects; we used da/dt = 0.1 AU Myr⁻¹.

We then define the Yarkovsky sensitivity s_Y as

$$\begin{equation} s_Y = \sqrt{\frac{1}{N}\sum _{i=1}^N \frac{\big(P_i^*-P_i^0\big)^2}{\sigma _i^2}}, \end{equation} \tag{ 3 }$$

where σ_i is the positional uncertainty associated with observation i. This root-mean-square quantity provides an excellent metric to assess the relative sensitivity of any given data set to a drift in semimajor axis, including drifts caused by Yarkovsky influences. The metric can be applied to the entire set of available observations, or to the subset of observations that survive the outlier rejection steps described below. We computed both quantities and used the latter for our analysis. We found that data sets with scores s_Y below unity yield unreliable results, including artificially large rates and large error bars. Out of ∼1250 numbered NEAs, only ∼300 have s_Y > 1 and ∼150 have s_Y > 2. In this paper, we focus on a subset of these NEAs.

For this work, we employed orbital fits to optical astrometry to determine semimajor axis drift rates for NEAs. We used the OrbFit software package, which is developed and maintained by the OrbFit Consortium (Milani & Gronchi 2009). OrbFit can fit NEA trajectories to astrometric data by minimizing the root mean square of the weighted residuals to the data, optionally taking into account a given non-zero rate of change in semimajor axis da/dt. We included perturbations from 21 asteroids whose masses were estimated by Konopliv et al. (2011).

We downloaded optical astrometry for all numbered minor planets (NumObs.txt.gz) from the Minor Planet Center (MPC) on 2012 January 31. We have assumed that all the astrometry has been properly converted to the J2000 system. The quality of the astrometry varies greatly, and we applied the data weighting and debiasing techniques implemented in OrbFit, which appear to follow the recommendations of Chesley et al. (2010). Data weights are based on the time the observation was performed, the method of the observation (CCD or plate), the accuracy of the star catalog, and in some cases the accuracy of the observatory. Correction for known star catalog biases was applied when possible. Biases vary depending on the specific star catalog and region of the sky, and can reach 1.5 arcsec in both right ascension (R.A.) and declination (decl.). Correction for these biases can substantially improve the recovery of orbital parameters from observations. However, as discussed in Chesley et al. (2010), not every observation can be debiased. Some observations were reported to the MPC without noting the star catalog used in the data reduction. Although Chesley et al. (2010) deduced the star catalogs used by several major surveys, there remain observations from smaller observatories that do not have associated star catalogs. Accordingly, a fraction of the astrometry used in this paper was not debiased. Based on counts published Chesley et al. (2010), we estimate this fraction to be less than 7.2% of all the observations.

Our procedure for determining the semimajor axis drift rate included three steps: an initial fit to the debiased data, an outlier rejection step, and a search for the best-fit da/dt, with iteration of the last two steps when necessary.

We used the orbital elements from the Minor Planet Center's MPCORB database as initial conditions for the first fit for each object (step 1). This first fit, performed with da/dt = 0 and outlier rejection turned off, slightly corrected the orbital elements for our weighted, debiased observations. The orbital elements from each object's first fit became the starting orbital elements for all later fits of that object.

The second fit of each object served to reject outliers and was initially performed with da/dt = 0 (step 2). The residual for each observation was calculated using the usual observed (O) minus computed (C) quantities:

$$\begin{equation} \hbox{\fontsize{9}{11}\selectfont$ \chi _{\rm res}=\sqrt{\left(\dsty\frac{({\rm R.A.}_{\rm O}-{\rm R.A.}_{\rm C}) \times \rm cos({\rm decl.}_{\rm O})}{\sigma _{{\rm R.A.}} }\right)^2+\left(\dsty\frac{{{\rm decl.}}_{\rm O}-{{\rm decl.}}_{\rm C}}{\sigma _{{{\rm decl.}}}}\right)^2},$} \end{equation} \tag{ 4 }$$

where σ_R.A. and σ_decl. are the uncertainties for that observation in R.A. and decl., respectively. We rejected observations when their $\chi _{\rm res} >\sqrt{8}$, and recovered previously rejected observations at $\chi _{\rm res}=\sqrt{7}$, with the rejection step iterated to convergence. Results are fairly robust over a large range of thresholds for rejection (Section 3). If the post-fit residuals were normally distributed, the chosen thresholds would result in <1% of observations being rejected as outliers. Because errors are not normally distributed, our typical rejection rates are 2%–5% of all available astrometry. This second fit produced the set of observations that were used in the third step.

The third step was a series of orbital element fits to the observations over a set of fixed da/dt values. During these fits, we used the set of observations defined by the second fit and did not allow further outlier rejection. The quality of a fit was determined by summing the squares of residuals χ² = ∑χ²_res. To locate the region with the lowest χ², we used a three-point parabolic fit or the golden-section minimization routine (Press et al. 1992). A parabola was then fit to the χ² curve in the vicinity of the minimum, and we used the minimum of the parabola to identify the best-fit da/dt value.

Confidence limits were estimated using χ² statistics. Confidence regions of 68.3% and 95.4% (1σ and 2σ, respectively) were established by the range of da/dt values that yielded χ² values within 1.0 and 4.0 of the best-fit χ² value, respectively (Figure 1).

Zoom In Zoom Out Reset image size

The initial outlier rejection step can in some cases eliminate valid observations simply because the Yarkovsky influences are not captured in a dynamical model with da/dt = 0. To circumvent this difficulty, we iterated the outlier rejection step with the best-fit da/dt value and we repeated the fitting process. In 52 out of 54 cases, the new best-fit value matched the previous best-fit value to within 1σ, and we accepted the new best-fit values as final. For the other objects we repeated the reject and fit processes until successive best-fit values converged within 1σ (which never required more than one additional iteration). Our results report the da/dt values obtained at the end of this iterative process.

We restricted our study to numbered NEAs with the best Yarkovsky sensitivity (Equation (3)), specifically s_Y > 2 (Figure 2).

Zoom In Zoom Out Reset image size

We also chose to focus on objects with non-zero da/dt values by using a signal-to-noise ratio (S/N) metric, defined as the ratio of the best-fit da/dt to its 1σ uncertainty. We accepted all objects with S/N >1 (Figure 2).

Some asteroids have observations that precede the majority of the object's astrometry by several decades and have relatively high uncertainties. In order to test the robustness of our results, we removed these sparse observations, which were defined as 10 or fewer observations over a 10 year period. Fits were then repeated for these objects without the early observations. If the initial best-fit value fell within the 1σ error bars of the new best-fit value, the initial result was accepted, otherwise, the object was rejected.

Superior detections of the Yarkovsky effect are likely favored with longer observational arcs, larger number of observations, and good orbital coverage. For this reason we limited the sample to those NEAs with an observational arc at least 15 years long, with a number of reported observations exceeding 100, and with at least eight observations per orbit on at least five separate orbits.

We report on the 54 objects that met all of these criteria: sensitivity, S/N, sparse test, and orbital coverage.

We validated our optical-only technique whenever radar ranging observations were available on at least two apparitions. This could only be done for a fraction of the objects in our sample. In the remainder of this paper, optical-only results are clearly distinguished from radar+optical results. For the radar+optical fits, we included all available radar astrometry and disallowed rejection of potential radar outliers. The internal consistency of radar astrometry is so high that outliers are normally detected before measurements are reported.

We also verified that a fitting procedure that holds successive da/dt values constant is equivalent to performing seven-parameter fits (six orbital parameters and da/dt simultaneously). The da/dt values obtained with both procedures are consistent with one another.

In addition to the measurements described above, we produced numerical estimates of the diurnal Yarkovsky drift for each of the objects in our sample. Comparing the measured and estimated rates provides a way to test Yarkovsky models. In some instances, e.g., robust observations irreconcilable with accurate Yarkovsky modeling, it could also lead to the detection of other non-gravitational forces, such as cometary activity. Our numerical estimates were generated as follows. At each timestep, we computed the diurnal Yarkovsky acceleration according to Equation (1) of Vokrouhlický et al. (2000), which assumes a spherical body, with the physical parameters (Opeil et al. 2010) listed in Table 1 and an assumption of 0° or 180° obliquity. We assumed that the thermal conductivity did not have a temperature dependence, but found that adding a temperature-dependent term according to the prescription of Hütter & Kömle (2008; K = K₀ + K₁T³, with K₁ = 0.0076) did not change our predictions by more than 1%. We then resolved the acceleration along orthogonal directions, and used Gauss' form of Lagrange's planetary equations (Danby 1992) to evaluate an orbit-averaged da/dt.

The physical parameters chosen for these predictions mimic two extremes of rocky asteroids; one is intended to simulate a rubble pile with low bulk density, the other a regolith-free chunk of rock (Table 1). These parameters correspond to a thermal inertia range of 77–707 J m⁻² s^−0.5 K⁻¹, enveloping the results of Delbó et al. (2007), who found an average NEA thermal inertia to be 200 J m⁻² s^−0.5 K⁻¹. In most cases, the drift rates produced by these two extreme cases encompass the drift produced by a rubble-pile object that has a regolith-free surface, or the drift produced by a solid object with regolith.

There is no simple relationship between these physical parameters and predicted drift rates, but for most cases the rubble pile exhibits the larger da/dt values due to its low bulk density (Equation (2)). The smaller values of density of the surface and thermal conductivity for rubble piles produce a smaller thermal inertia, and therefore a longer thermal lag. Generally, but not always, this longer thermal lag, combined with the rotation of the asteroid, allows for a larger fraction of departing thermal emission to be aligned with the asteroid's velocity, resulting in a larger drift.

When available, measured values of the geometric albedo, diameter, and spin rate from the JPL Small-Body Database (Chamberlin 2008) were incorporated into our predictions for Yarkovsky drifts. When not available, the diameter D in km was estimated from the absolute magnitude H using (Pravec & Harris 2007)

$$\begin{equation} D=\frac{1329}{\sqrt{p_{V}}} * 10^{-0.2H}, \end{equation} \tag{ 5 }$$

where we used two values of the V-band geometric albedo p_V (0.05 and 0.45), a range that captures observed albedos for the majority of NEAs. When spin rate was unknown, we assumed a value of 5 revolutions day⁻¹, based on the average spin rate values for asteroids 1–10 km in diameter shown in Figure 1 of Pravec & Harris (2000). Emissivity was assumed to be 0.9. Bond albedo was estimated with a uniform value of the phase integral (q = 0.39) on the basis of the IAU two-parameter magnitude system for asteroids Bowell et al. (1989) and an assumed slope parameter G = 0.15.

We have assumed p_V = 0.14 for the purpose of quantifying the Yarkovsky efficiency when the asteroid size was unknown.

Because a clear connection exists between asteroid physical properties and Yarkovsky drifts, we explored the constraints that can be placed on bulk density and surface thermal conductivity for seven objects with well-known diameters and (excepting one case) spin periods: (1620) Geographos, (1862) Apollo, (2100) Ra-Shalom, (2062) Aten, (2340) Hathor, (1566) Icarus, and (3361) Orpheus. We compared the measured Yarkovsky rates to numerical estimates obtained with a range of physical parameters. For these estimates, we assumed a constant heat capacity C = 500 J kg⁻¹ K⁻¹ (Table 1) and a single value of the bulk density of the surface ρ_s = 1700 kg m⁻³, but we explore a wide range of bulk density and surface thermal conductivity values. Because the obliquities are uncertain or ambiguous in many cases, we chose to illustrate outcomes for two obliquity values, typically 180° and 135°.

Our results are shown in Figures 7 and 8, which are similar to Figure 4 in Chesley et al. (2003). The shaded range consistent with the 1σ confidence limits on da/dt delineates the space of acceptable bulk densities and thermal conductivities, assuming that the Yarkovsky effect is being modeled correctly. By acceptable, we mean consistent with observed da/dt values, even though some of the K–ρ values may not be appropriate for asteroids.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Infrared observations indicate that (2100) Ra-Shalom has a thermal conductivity between 0.1 and 1 W m⁻¹ K⁻¹ (Delbó et al. 2003; Shepard et al. 2008). If we assume a minimum bulk density of 1500 kg m⁻³, this conductivity value is consistent with the range suggested by our Yarkovsky rate determination.

If we make the same minimum density assumption for (1620) Geographos, our measurements suggest that its surface thermal conductivity is greater than 0.002 W m⁻¹ K⁻¹.

For (1862) Apollo, we show the range of physical properties that are consistent with both the optical-only fits and the radar+optical fits. The precision of the radar measurements dramatically shrinks the size of the measured error bars, with correspondingly tighter constraints on density and surface thermal conductivity. This example illustrates that reliable obliquity determinations will be important to extract physical properties from Yarkovsky rate determinations.

Our measurement of (2062) Aten's drift provides some useful insights. If we assume that its bulk density exceeds 1500 kg m⁻³, then its surface thermal conductivity K must exceed 0.3 W m⁻¹ K⁻¹. Furthermore, if we assume that its bulk density exceeds 1600 kg m⁻³, the 1σ confidence region on the measured Yarkovsky drift suggests that its obliquity is between 180° and 135°.

The Yarkovsky simulations for (2340) Hathor were computed with an assumed spin period of 4.5 hr. If the actual period is longer, the curves shown would shift to the left, and if the period is shorter, the curves would shift to the right. Consequently, we cannot make inferences about the K value for this object until its spin period is measured. However, looking at the height of the curve, and with an assumption that the object's bulk density is greater than 1500 kg m⁻³, we can conclude that (2340) Hathor likely has an obliquity lower than 180°.

The assumption of 135° or 180° obliquity for (1566) Icarus restricts this object to low surface conductivity values and low bulk density values, or high surface conductivity values and high bulk density values. Although these obliquities do produce physically plausible parameter combinations, it seems likely that the obliquity for this object is ⩽135°.

The curves for (3361) Orpheus were calculated with an assumed geometric albedo of 0.15. As (3361) Orpheus has a positive da/dt value, obliquities were assumed to be 0° and 45°. The curve representing an obliquity equal to 0° for this object requires very low (<0.002 W m⁻¹ K⁻¹) or very high (>0.7 W m⁻¹ K⁻¹) surface thermal conductivity values for most densities. A more likely scenario is that this object has an obliquity >0°, or perhaps even >45°. An independent measurement of the obliquity could be used to validate obliquity constraints derived from Yarkovsky measurements.

La Spina et al. (2004) and Chesley et al. (2008) examined the predominance of retrograde spins and negative Yarkovsky drift rates and concluded that they were consistent with the presumed delivery method of NEAs from the main belt of asteroids. The ν₆ and 3:1 resonance regions deliver NEAs to near-Earth space (Bottke et al. 2002a). A main belt asteroid can arrive at the 3:1 resonance at 2.5 AU via a positive (if it originates in the inner main belt) or negative (if it originates in the outer main belt) Yarkovsky drift. However, a main belt asteroid can only arrive at the ν₆ resonance (at the inner edge of the main belt) by the way of a negative drift. According to Bottke et al. (2002a) and Morbidelli & Vokrouhlický (2003), 30%–37% of NEAs are transported via the ν₆ resonance, with the rest from other resonances. The net result is a preference for retrograde spins.

An observational consequence of this process would be an excess of retrograde rotators in the NEA population. La Spina et al. (2004) conducted a survey of 21 NEAs and found the ratio of retrograde/prograde rotators to be 2.0⁺¹_{− 0.7}.

We note that out of the 42 Yarkovsky-dominated NEAs, 12 have a positive da/dt value. For this sample, our ratio of retrograde/prograde rotators is 2.5 ± 0.1, similar to the value found by La Spina et al. (2004).

The semimajor axis drifts described in this paper affect NEA trajectory predictions. An order of magnitude estimate for the along-track displacement due to a non-zero da/dt is given in Vokrouhlický et al. (2000):

$$\begin{equation} \Delta \rho \simeq 7 \dot{a}_4 (\Delta _{10} t)^2 a_{{\rm AU}}^{-3/2}, \end{equation} \tag{ 6 }$$

where Δρ is in units of km, $\dot{a}_4$ is da/dt in 10⁻⁴ AU Myr⁻¹, Δ₁₀t is the time difference between observations in tens of years, and a_AU is the semimajor axis of the object in AU. For instance, the estimated along-track displacement due to the observed da/dt for (1862) Apollo is 9 km after 10 years. Similarly, the estimated along-track displacement for faster-moving (1864) Daedalus is 67 km after 10 years.

Our data indicate that (101955) 1999 RQ36, the target of the OSIRIS-REx mission, has a measurable Yarkovsky drift of (− 18.9 ± 0.2) × 10⁻⁴ AU Myr⁻¹. Although it has a relatively short arc (12 years) it has been observed three times by radar, allowing for an accurate da/dt measurement. We estimated the along-track displacement of (101955) 1999 RQ36 over the 6 month duration of the OSIRIS-REx mission to be 0.3 km, which will be easily detectable by a radio science instrument.

(1862) Apollo is a binary asteroid (Ostro et al. 2005). Binary asteroids present a unique opportunity for the determination of physical parameters. If mass and density can be measured from the binary orbit and component sizes, the Yarkovsky constraint on thermal conductivity can become much more meaningful. If the orientation of the plane of the mutual orbit can be measured, a plausible obliquity can be assumed, which makes the constraints on thermal properties tighter still. In some cases, actual obliquity measurements can be obtained from shape modeling efforts.

Yeomans (1991, 1992) identified a non-gravitational perturbation acting on the orbit of (1862) Apollo, but was not able to determine a drift magnitude. To 1σ our observed da/dt value for (1862) Apollo agrees with our Yarkovsky predictions.

(1036) Ganymed has by far the largest Yarkovsky efficiency value (f_Y ∼ 15 × 10⁻⁵) among the objects presented in Table 3. With a nominal value of ∼ − 7 × 10⁻⁴ AU Myr⁻¹, the measured da/dt value is comparable to that of other NEAs. Combined with Ganymed's large diameter estimate (∼32 km based on IRAF measurements), this Yarkovsky rate results in an unusually high f_Y value.

How can this anomaly be explained? One possibility is that some of the early astrometry, dating back to 1924, is erroneous. This could be due to measurement errors, timing errors, bias errors, or reference frame conversion errors. We evaluated the semimajor axis drift with various subsets of the available astrometry and found values ranging between −3 × 10⁻⁴ and −8 × 10⁻⁴ AU Myr⁻¹. On that basis we modified the adopted uncertainties for this object, and our preferred value is (− 6.62^+3.6_{− 1.4}) × 10⁻⁴ AU Myr⁻¹. Doing so does not eliminate the possibility of systematic bias in the astrometry, and we are still left with anomalously high f_Y values.

Another possibility is that the diameter of Ganymed, an S-type asteroid, is much smaller than reported. This seems unlikely considering the more recent WISE albedo measurement of p_V = 0.212 (Masiero et al. 2012), which suggests a diameter of ∼36 km.

If Ganymed's bulk density was especially low, a higher than usual f_Y value would be expected, but this would likely explain a factor of two or three at most, and would not explain the anomalous value.

Perhaps Ganymed departs significantly from a spherical shape, with an effective diameter and mass that are much smaller than those implied by the diameter values reported in the literature. The relatively low light-curve amplitudes do not seem to support such an argument, unless the asteroid is particularly oblate. In that case one could plausibly arrive at volume and mass estimates that are off by a factor of 5–10.

If we can rule out these possibilities (i.e., Ganymed is roughly spherical with no substantial concavities, its diameter estimate is reasonably accurate, and the early astrometry can be trusted), and if no other modeling error can be identified, then we would be compelled to accept an anomalously high Yarkovsky efficiency for this object.

In the course of our study we observed drift values that cannot be accounted for easily by Yarkovsky drift, because they considerably exceed the predicted Yarkovsky rates. In most cases, these can be attributed to poor sensitivity to Yarkovsky influences (Figure 2). Therefore, the high rates can generally be safely discarded. In other cases, the high rates may be due to erroneous optical astrometry or mismodeling of asteroid–asteroid perturbations. However, we cannot entirely rule out the possibility that some of the high drift rates are secure and will be confirmed by further observation and analysis. If the high rates cannot be ascribed to poor Yarkovsky sensitivity or faulty astrometry, one would need to invoke other non-gravitational forces.

One possibility is that orbits are perturbed when NEAs are losing gas or dust in an anisotropic manner. To estimate a rough rate of mass loss that would be needed to account for the drifts measured, we used the basic thrust equation

$$\begin{equation} F=qV_{e}, \end{equation} \tag{ 7 }$$

where F is the force, q is the rate at which the mass departs the asteroid, and V_e is the ejection speed. For an asteroid of mass m this yields

$$\begin{equation} a_{\rm mass\ loss}=\frac{qV_{e}}{m} \end{equation} \tag{ 8 }$$

which can be incorporated into Gauss' form of Lagrange's planetary equations (Danby 1992) as an acceleration aligned with the velocity of the object. The dependence of the force on heliocentric distance r is not known precisely; we assumed F∝r⁻², similar to the Yarkovsky dependence, for simplicity, and because the amount of outgassing likely scales with the amount of incident radiation (as in Figure 4 of Delsemme 1982). We assumed V_e = 1.5 m s⁻¹, the value derived by Hsieh et al. (2004) for 133P/Elst-Pizarro, and we assumed that the mass is departing in the optimal thrust direction.

We quantified the mass-loss rates needed to produce the observed drifts of NEAs with the highest Yarkovsky efficiencies. We estimated a rate of 0.16 kg s⁻¹ for (154330) 2002 VX94 and 2.3 kg s⁻¹ for (7889) 1994 LX. Although these estimates represent the minimum amount of mass loss necessary to account for the observed drifts (if due to mass loss), they are smaller than typical levels from comets. Comets have mass-loss rates that span a wide range of values. On the high side a rate of 2 × 10⁶ kg s⁻¹ was estimated for Hale–Bopp (Jewitt & Matthews 1999). On the low side Ishiguro et al. (2007) measured mass-loss rates for three comets, averaged over their orbits: 2P/Encke (48 ± 20 kg s⁻¹), 22P/Kopff (17 ± 3 kg s⁻¹), and 65P/Gunn (27 ± 9 kg s⁻¹). Mass-loss rates of active asteroids have been estimated to be in the range from ⩽0.04 kg s⁻¹ (113P/Elst-Pizarro) to ⩽150 kg s⁻¹ (107P/Wilson-Harrington) (Jewitt 2012).

Mass loss does not seem to be a viable mechanism to explain the semimajor axis drift rate of (1036) Ganymed, as it would require a minimum mass-loss rate of ∼2500 kg s⁻¹. This would presumably have left detectable observational signatures, which have not been reported to date.

We explore a couple of possibilities for mass-loss mechanisms that could cause semimajor axis drifts.

4.6.1. Associations with Meteoroid Streams

4.6.2. Rock Comet Phenomenon