© ESO, 2017

1. Introduction

The complex motion of near-Earth asteroids (NEAs) is dominated by the gravitational perturbations of the Sun and planets. However, the gravitational interaction with other small bodies and non-gravitational perturbations can affect their behavior and become relevant for the prediction of their future positions (Farnocchia et al. 2015).

In particular, the Yarkovsky effect is a subtle non-gravitational acceleration due to the anisotropic emission of thermal radiation of Solar System objects that causes a secular drift in the semi-major axis (Bottke et al. 2006). This perturbation is important to understand the long-term dynamics of the asteroid population since it is a driving factor for feeding resonances in the main belt and transporting asteroids to the inner Solar System (Farinella et al. 1998; Morbidelli & Vokrouhlickỳ 2003; Bottke et al. 2002b). The Yarkovsky effect is also relevant for impact hazard predictions where high-precision ephemeris predictions are required (Giorgini et al. 2002, 2008; Chesley 2006; Farnocchia et al. 2013b; Farnocchia & Chesley 2014; Chesley et al. 2014; Spoto et al. 2014; Vokrouhlický et al. 2015b).

The diurnal component of the Yarkovsky effect, which is usually the dominant one, is proportional to the cosine of the obliquity of the spin axis (Bottke et al. 2006). Therefore, the spin orientation determines whether an asteroid’s semi-major axis drifts inwards or outwards. More than ten years ago, La Spina et al. (2004) analyzed the distribution of known NEA spin axes, about 21 at the time, and found a ratio of retrograde to direct rotators. The observed ratio was an excellent match to the one expected from the Bottke et al. (2002a) NEA population model and the injection mechanism of asteroids to the inner Solar System through orbital resonances, that is, 2 ± 0.2.

A derivation of the Yarkovsky accelerations from thermophysical modeling is generally impractical as they depend on physical properties such as size, mass, shape, obliquity, and thermal properties (Bottke et al. 2006) and even the surface roughness (Rozitis & Green 2012), which are generally unknown. However, for asteroids with a well-observed astrometric arc, it is possible to directly estimate the Yarkovsky effect by measuring deviations from a gravity-only trajectory (Chesley et al. 2003, 2016; Vokrouhlický et al. 2008, 2015a; Nugent et al. 2012; Farnocchia et al. 2013a, 2014).

We currently have a limited dataset of known NEA spin axes; the EARN¹ and the DAMIT² databases combined list 43 as of Mar 9, 2017. Thus, it is difficult to derive a statistically reliable obliquity distribution. In general, one needs specific observations at multiple observing geometries to constrain the spin axis of an asteroid, for example, light curves (Ďurech et al. 2009, 2011) or radar observations (Benner et al. 2015). Even so, in many cases there are two distinct solutions and it is not always possible to identify the correct one. An interesting example is (29075) 1950 DA, which had two possible pole solutions from radar observations (Busch et al. 2007). The estimate of the Yarkovsky effect on this object (Farnocchia & Chesley 2014) resolves the ambiguity between the two in favor of the retrograde solution.

Along the same lines, in this paper we use a current dataset of Yarkovsky estimates to put constraints on the NEA obliquity distribution. In particular, by using the properties of the NEA population we can derive distributions from most of the parameters on which the Yarkovsky effect depends. By making numerous assumptions, for example, neglecting the dependence of the Yarkovsky effect on bulk density and thermal properties, Farnocchia et al. (2013a) made a previous attempt to infer a four-bin NEA obliquity distribution from a set of 136 Yarkovsky detections. In this paper we use a more sophisticated technique by solving an inverse problem where different obliquity distributions are tested to provide the best match to the Yarkovsky estimate dataset. This technique was introduced, with a preliminary application to a similar dataset, in Cotto-Figueroa (2013).

2. Yarkovsky modeling

The Yarkovsky perturbation can be modeled as a transverse acceleration A₂/r² (Farnocchia et al. 2013a), where r is the distance from the Sun in au and A₂ is the sum of two terms, one corresponding to the diurnal effect due to the asteroid’s rotation and one corresponding to the seasonal effect due to the asteroid’s orbital motion: (1)where A is the Bond albedo, Φ(1au) = 3G_S/ (2Dρc) is the standard radiation force factor at 1 au, G_S = 1361 W/m² is the solar constant, D is the asteroid’s diameter, ρ is the bulk density, and γ is the spin obliquity. The thermal parameters θ_rot and θ_rev depend on the rotation and revolution periods, respectively, and also on spin rate, thermal inertia, thermal emissivity, geometric albedo p_v, and r. The function f describes the spin-rate and thermal-inertia dependence of the Yarkovsky acceleration for a Lambertian sphere, it is non-monotonic and f(0) = f(∞) = 0 (Bottke et al. 2006). Finally, α is an enhancement factor that is intended to describe the effect of surface roughness alone, but that effect is itself dependent on thermal inertia and spin rate (Rozitis & Green 2012).

By separating the obliquity γ from all of the other parameters we have (2)where F₁,F₂ are positive functions that do not depend on the obliquity. To simplify we replace the instantaneous heliocentric distance with the solar flux-weighted mean heliocentric distance , where a and e are orbital semi-major axis and eccentricity, respectively. Equation (2)represents the starting point of our inverse process to derive possible obliquity distributions starting from probability distributions for A₂, D, A, ρ, Γ, , P_rot, P_rev, and α.

3. Dataset of Yarkovsky estimates

To obtain a distribution of A₂ on the left-hand side of Eq. (2)we used the Chesley et al. (2016) list of Yarkovsky estimates. This list contains 42 Yarkovsky detections considered “valid”, which means that the signal-to-noise ratio (SNR) of the detection is greater than 3 and its value is compatible with the Yarkovsky mechanism.

Moreover, Chesley et al. (2016) have a second category referred to as “weak” detections where the Yarkovsky estimate uncertainty is small enough that it would permit a clear detection if the Yarkovsky A₂ parameters were scaled from that of Bennu using its 1 /D dependence. However, the astrometric observations are incompatible with such accelerations thus suggesting a lower magnitude of the Yarkovsky effect. Some of Bennu’s physical properties tend to increase the Yarkovsky effect, for example, the extreme obliquity of 178° and the low bulk density of 1.3 g/cm³ (Chesley et al. 2014), thus the category of “weak” detections is likely to include objects that have physical properties (e.g., obliquity, bulk density) that lower the magnitude of the Yarkovsky effect. We included these “weak” detections in our dataset to avoid biasing the sample against non-extreme obliquities.

We updated the Chesley et al. (2016) list by including all of the available optical and radar astrometry as of September 2016, for a final dataset of 125 Yarkovsky estimates (see Table A.1). To limit the spread in A₂ caused by the diverse sizes (from a few meters to a kilometer in diameter) of the objects for which we have a Yarkovsky estimate, we normalized the A₂ values by an absolute magnitude scale factor 1329 km10^{− H/ 5}/ with a constant albedo p_v = 0.154. The resulting A₂ distribution is shown in Fig. 1.

Fig. 1 Histogram of Yarkovsky estimates normalized by an absolute magnitude scale factor of . The average 1σ uncertainty in normalized A2 is 5 × 10-15 au/d2.

4. Probability distribution of physical parameters

To invert Eq. (2)and solve for an obliquity distribution, we need to model the intrinsic distributions of the other parameters needed to compute the Yarkovsky effect, that is, D, A, ρ, Γ, , P_rot, P_rev, and α. The adopted distributions are based on what is known of the NEA population as well as the specific properties of the Yarkovsky estimate dataset.

Diameter. To derive the diameter we use the standard conversion formula from absolute magnitude H and geometric albedo p_v (Pravec & Harris 2007): (3)While the absolute magnitude distribution of NEOs follows a power law (Bottke et al. 2002a), the one of the objects in our dataset of Yarkovsky estimates resembles a normal distribution; see Fig. 2.

Fig. 2 Probability distribution of the absolute magnitude for the objects in Table A.1.

The shape of the distribution can be explained by the contribution of two factors: on one side there are fewer bigger objects in the population, and on the other side there are smaller objects, which are faint and so are harder to observe or even discover. Therefore, small objects are less likely to have long observation arcs, which reduces the chances of obtaining constraints on the Yarkovsky effect. The result is that objects with a magnitude around H = 20 are currently the ones more likely to have a Yarkovsky estimate. To sample H we selected a normal distribution with a mean of 20.12 and standard deviation of 2.44.

Geometric albedo. For the geometric albedo we consider three major taxonomic classes: C, S, and X. The split and frequencies are from Chesley et al. (2002) with the difference that we combine the M-class and the E-class into X. For each of the classes, we use lognormal distributions with median and standard deviation shown in Table 1, which are informed by the statical analysis presented in Thomas et al. (2011). The usage of lognormal distributions prevents the occurrence of nonphysical negative values of p_v. From the drawn values of H and p_v we sample the diameter using Eq. (3).

Among the objects of Table 1, only 39 objects have known taxonomy in the C, S, and X classes. The split is consistent with that of Table 1, in fact (7.7 ± 4.1)% are of type C, (74 ± 7.0)% are of type S, and (17.7 ± 6.2)% are of type X.

Bond albedo. The Bond albedo A is a function of the geometric albedo p_v and the slope parameter G: A = (0.29 + 0.684G)p_v (Bowell et al. 1989). We already described the distribution for p_v. Following Mommert et al. (2014a) and Mommert et al. (2014b), we analyzed the current statistics from the JPL Small-Body Database³ and obtained normal distribution for G with a mean 0.18 and a standard deviation 0.13. The geometric albedo was derived from the distributions described earlier.

Bulk density. Similar to what was done for the geometric albedo, for the bulk density ρ we considered lognormal distributions depending on the taxonomic class; see Table 1. The distribution parameters are based on the census of asteroid densities performed by Carry (2012). We note that, once the taxonomic class is drawn, the distributions of both p_v and ρ are chosen consistently with that taxonomic class, thus accounting for the fact that p_v and ρ are not independent.

Thermal inertia. To account for thermal inertia, we computed the thermal parameter for each NEA with a measured thermal inertia value from Table 2 of Delbò et al. (2015). We excluded (54509) YORP from this list because of the very large uncertainty on its derived thermal inertia value. A lognormal distribution was then fit to the 13 measured thermal parameter values to give the mean NEA thermal parameter as θ_rot = 10^{0.6 ± 0.3}. For our synthetic population, thermal parameters were then randomly drawn from this log-normal distribution.

Rotation period. The rotation period P_rot is size dependent. In particular, there is the so-called spin barrier of 2.2 h (Warner et al. 2009): few objects with H < 20 spin faster than this limit. Therefore, we divide the absolute magnitude in bins of 1 mag and for each bin we sample the rotation period according to the mean and standard deviation of the rotation periods available from the JPL Small-Body Database (most of which are from the Warner et al. 2009, asteroid light curve database) in that bin (see Fig. 3). To avoid nonphysical values of the rotation period, we removed samples with a period greater than 1000 h. Moreover, we removed the samples with a rotation period smaller than 2 h for H ≤ 20, or smaller than 0.01 h for H > 20.

Fig. 3 Median (crosses) and 1σ standard deviation of the rotation periods from the JPL Small-Body Database for bins of 1 mag from H = 10 to H = 30.

Orbital period and solar-flux-weighted heliocentric distance. Some orbital configurations favor a Yarkovsky estimate, for example, potentially hazardous asteroids (PHAs) are more likely to come close to Earth and so are easier to observe, especially using radar. In particular, none of the objects in the dataset has a perihelion q > 1.15 au or an aphelion Q < 1 au. Since the orbital geometry can introduce a selection effect, we took the distribution in semi-major axis a and e corresponding to our set of Yarkovsky detections: the semi-major axis distribution is approximated with a lognormal distribution whose associated normal distribution has mean 0.13 au and standard deviation 0.33 au, while the eccentricity distribution is approximated with a normal distribution with mean 0.4 and standard deviation 0.2 truncated at 0 and 1. Finally, we filtered out the objects with q > 1.15 au and Q < 1 au, and converted the semi-major axis to the orbital period P_rev (see Fig. 4).

The distribution in solar-flux-weighted heliocentric distance is derived from the a and e distributions described above.

Fig. 4 Orbital period distribution as derived from the objects in Table A.1.

Enhancement factor. Small-scale surface roughness enhances the diurnal component of the Yarkovsky effect through thermal-infrared beaming, that is, re-radiation of absorbed sunlight back towards the Sun (Rozitis & Green 2012). The degree of enhancement is a non-linear function of the asteroid thermal parameter, albedo, and heliocentric distance. In particular, it has been previously shown that the enhancement factor increases for decreasing thermal parameter and decreasing heliocentric distance, and it also increases for increasing albedo. The enhancement factor for a set of properties can be calculated for a spherically shaped asteroid covered with hemispherical craters (i.e., the craters represent the surface roughness) using the thermophysical model described in Rozitis & Green (2012).

Using this model, we generated a lookup table to obtain the enhancement factor as a function of A, θ_rot, and . The top panel of Fig. 5 shows the enhancement factor corresponding to a 100% roughness as a function of θ_rot.

Fig. 5 Above: enhancement factor α corresponding to 100% roughness as a function of the thermal parameter. The level curves correspond to a distance of 1 au. Below: histogram of the enhancement factor α.

Finally, we obtain α by scaling by the asteroid’s surface roughness, which is uniformly drawn between 0% and 100%. On average we obtain a 15% enhancement of the diurnal component of the Yarkovsky effect. The bottom panel of Fig. 5 shows the resulting distribution of α, obtained from the drawn A, θ_rot, and and interpolating the lookup table.

5. Models for the obliquity distribution

We considered three different parametric models for the distribution of the cosine of the obliquity. These parametric formulations enable us to generate synthetic A₂ distributions to be compared to the Yarkovsky dataset of Sect. 3 and in turn find the ones providing the best match. Because of the YORP effect, we expect local maxima of the distribution for extreme obliquities (0° and 180°) and a minimum close to γ = 90° (Čapek & Vokrouhlický 2004).

Bin model. The simplest model considers a number of bins N and tests different frequencies p_i,i = 1,N for each bin. Each p_i has to be positive and their sum has to be 1. The number of degrees of freedom is N−1. Figure 6 shows an example with N = 3.

Fig. 6 Example of a three-bin distribution.

Piecewise linear model. The second model we considered is that of a continuous piecewise linear (PL) function, as shown in Fig. 7. The three parameters defining this distribution are the ordinates Q₁, Q₂, and Q₃ in −1, 0, and 1, respectively. Since the integral has to be one, that is, Q₁ + 2Q₂ + Q₃ = 2, the number of independent parameters is 2. This model can be generalized by having a variable abscissa x₂ for the middle point, thus leading to three independent parameters. We refer to this generalization as PLMP.

Fig. 7 Example of a piecewise linear distribution.

Quadratic model. The final model we consider is that of a quadratic function f(γ) = acos²γ + bcosγ + c (see Fig. 8). We allow only concave-up parabolas, that is, a> 0, as we know that the YORP effect favors extreme obliquities (Čapek & Vokrouhlický 2004). The parabola’s minimum must be positive and its abscissa between −1 and 1, and the integral has to be 1. Therefore, the number of independent parameters is 2.

Fig. 8 Example of a quadratic distribution.

6. Solution of the inverse problem

Starting from the distributions described in Sect. 4 and a given parametric distribution in the obliquity, we can draw samples and use Eq. (2)to obtain samples in A₂ (see Fig. 1). Therefore, for each parametric obliquity distribution we obtain a predicted distribution in A₂ to be compared with that coming from the set of Yarkovsky estimates. As already described in Sect. 3, the A₂ values are normalized by absolute magnitude to reduce the spread caused by the range of different sizes considered.

To measure how well the predicted distribution matches the one from the Yarkovsky estimates we perform a χ² test. The range of normalized A₂ values is divided in m bins so that each bin contains the same probability mass from the predicted distribution, that is, the integral of the predicted distribution over each bin is 1 /m. Then, from the list of Yarkovsky estimates we compute the probability p_i,i = 1,m of falling within each bin. Finally, χ² is computed as: (4)We look for the obliquity distributions providing a lower χ², which from standard χ² statistics should be close to m−1−n_p, where n_p is the number of estimated parameters. To avoid small number statistics and after checking the stability of χ², we based our predicted distribution on 10⁵ samples and the number of bins is m = 11.

7. Results

We first test the obliquity distribution parametric models with two independent parameters, that is, a three-bin distribution, a piecewise linear function with middle point in 0, and a quadratic function. To find the best fitting parameters, we scan a grid and compute χ² for each grid point. Figures 9−11 show the level curves of χ² as a function of the parameters.

Fig. 9 Discrete three-bin distribution. Thin lines indicate level curves of χ2 in the two-dimensional (2D) search space. The minimum χ2 is marked with a diamond. Thick lines mark loci of constant retrograde to direct rotator ratios. The dotted area corresponds to inadmissible values of the parameters.

Fig. 10 Level curves of χ2 and ratio for the PL distribution model with two parameters. Lines and markers are the same as in Fig. 9.

Fig. 11 Quadratic distribution: level curves of χ2 and ratio in the 2D search space. Lines and markers are the same as in Fig. 9.

Table 2 gives the best-fit χ², the p-value, that is, the probability of randomly obtaining a larger χ², and the corresponding ratio between retrograde and direct rotators. All the models provide statistically acceptable p-values, with a quadratic model giving the best fit with a χ² of 7.4.

All of the models give a retrograde to direct rotators ratio (R_R/D) that is statistically consistent with the ratio found by La Spina et al. (2004) and also with the theoretical 2 ± 0.2 ratio derived from NEA population models (Bottke et al. 2002a), which suggests that NEAs generally maintain their rotation direction. However, if the timescales required to complete a YORP cycle (Rubincam 2000) are much shorter than the typical NEA dynamical lifetime, the YORP effect should have randomized the distribution of prograde versus retrograde rotators. YORP self-limitation (Cotto-Figueroa et al. 2015) may provide a means to reconcile the high present-day retrograde fraction, where the YORP-driven deformation of aggregate objects confines their rotation rates to far narrower ranges than would be expected in the classical YORP-cycle picture, therefore greatly prolonging the times over which objects can preserve their sense of rotation.

We now try to increase the number of independent parameters to three for the bin and piecewise linear models. Therefore, the expected χ² decreases to 7. To find an absolute minimum value of χ² on a three-dimensional (3D) space we use the IDEA global optimizer (Vasile et al. 2011). For a four-bin distribution we find a best fit solution (p₁,p₂,p₃,p₄) = (0.56,0.15,0.06,0.28) with minimum χ² = 9.55 (p-value of 22%) and R_R/D = 2.4. An F-test shows that the Δχ² = 13.60−9.55 improvement due to the addition of a third parameter is only significant at the 13% level. For the generalized piecewise linear function with variable mid-point abscissa the best fit solution is (Q₁,Q₂,Q₃,x₂) = (1.16,0.04,0.73,0.074) with minimum χ² = 8.14 (p-value of 32%) and R_R/D = 1.8. Again, the statistical significance of Δχ² = 10.43−8.14 as measured by an F-test is only 20%.

Therefore, the addition of a third parameter is only marginally significant and the quadratic model with two independent parameters provides the lowest χ². This suggests that our set of Yarkovsky estimates does not have enough signal to solve for more than two parameters.

Figure 12 shows the best-fit distribution obtained with the different models. Interestingly, all of the models favor extreme obliquities, which is consistent with the expected consequences of the YORP mechanism (Čapek & Vokrouhlický 2004). Figure 13 shows the corresponding distributions in normalized A₂ with that coming from the Yarkovsky estimates. All the models capture the bimodal behavior of the dataset of Yarkovsky estimates and find a larger fraction of negative A₂ values, which is consistent with the preliminary results from Cotto-Figueroa (2013). The distributions from our models overestimate the heights of the peaks of the A₂ distribution. Besides statistical noise due to Poisson statistics in the Yarkovsky sample, this effect could be caused by inaccurate assumptions in the distributions adopted in Sect. 4. For instance, the peaks drop by 10% when doubling the standard deviation of the density probability distributions in Table 1. However, the qualitative results on the obliquity distribution still stand; for example, the ratio between retrograde and direct rotators does not change.

Fig. 12 Best-fit obliquity distributions for the three-bin, piecewise linear, quadratic, and four-bin models.

Fig. 13 Comparison between the probability distributions in normalized A2 obtained from the Yarkovsky estimates and that from the different obliquity distribution models.

Figure 14 compares the Farnocchia et al. (2013a) result with the distributions found in this paper. For this purpose we convert our best-fit distribution to a four-bin one by computing their integral over each of the four bins. Interestingly, Farnocchia et al. (2013a) found a ratio of 2.7 between retrograde and direct rotators, while our results are closer to the theoretical expectation, which further suggests that the technique used in this paper has higher fidelity.

To validate our results, we consider the known spin axes of 41 NEOs from the EARN database⁴ and of 17 objects (15 already in the EARN database) from the DAMIT database⁵. Table 3 reports the fractions of obliquities with cosγ< −1/3, −1/3 < cosγ< 1/3, or cosγ> 1/3. Since the objects in our dataset of Yarkovsky estimates had H in the range 14–29 with a peak around H = 20, Table 3 also shows the split for different H ranges. Larger objects (H< 18) have more mass in the middle bin than in that with cosγ> 1/3. On the other hand, smaller objects have a smaller fraction of mid-range obliquities, which is again consistent with the YORP mechanism. It is interesting to point out that all of our best-fit solutions give a split that is well within 1σ of the observed obliquity distribution for objects with H> 18, which is also a condition that most objects of the Yarkovsky estimate dataset meet.

Fig. 14 Comparison between the obliquity distribution presented by Farnocchia et al. (2013a) and the four-bin split of our the best-fit obliquity distributions converted to four bins. In the legend we recall the ratio between retrograde and direct rotators RR/D.

Finally, we check the significance of the Rozitis & Green (2012) enhancement factor α of the diurnal component of the Yarkovsky effect, discussed in Sect. 4. Table 4 gives the best-fit χ² and p-values for the different obliquity distributions obtained by setting α = 1. For all of the models, the low p-values further support the presence of a Yarkovsky enhancement factor as predicted by Rozitis & Green (2012).

8. Conclusions

We used the Chesley et al. (2016) list of Yarkovsky estimates to derive constraints on the obliquity distribution of near-Earth asteroids. We solved an inverse problem where we adopted probability distributions on the physical parameters other than obliquity (e.g., albedo, thermal inertia, bulk density) that determine the Yarkovsky effect. Then, we considered different parametric models for the probability density function of the cosine of the obliquity: piecewise constant (i.e., bins), piecewise linear, and quadratic functions. Finally, we performed χ² tests to quantify the goodness of the match to the distribution of the Yarkovsky estimates.

Using models with only two independent parameters seems to be the best compromise between model complexity and goodness of the fit. F-tests show that a third additional parameter is only marginally significant, thus suggesting that the dataset of Yarkovsky detections does not have enough resolution to constrain a third parameter. Among the analyzed obliquity distributions with two parameters, the one that produces the best-fit is obtained with a quadratic model: f(cosγ) = 1.12cos²γ−0.32cosγ + 0.13. This solution favors extreme obliquities, which is consistent with the action of the YORP effect. Moreover, the corresponding ratio between retrograde and direct rotators, , provides an excellent match to the La Spina et al. (2004) result, that is, , and the Bottke et al. (2002a) theoretical expectation for feeding asteroids into the inner solar system, that is, 2 ± 0.2. Finally, this distribution is very much consistent with the set of known spin axis orientations.

It is possible that the results we obtained are affected by selection effects. In fact, the starting dataset of Yarkovsky detections is comprised of objects with a strong observational dataset, possibly including radar observations. Therefore, objects that have a more favorable observing geometry such as PHAs, dominate the dataset. In particular, all the objects had an aphelion larger than 1 au and so it is possible that Interior-Earth-Objects have a different obliquity distribution than what we derived. Future work will include assessing how this sort of orbital selection bias affects the extrapolation to the entire near-Earth asteroid population.

We also find that the current sweet spot for Yarkovsky detections is around H = 20, with most of the Yarkovsky estimates having an absolute magnitude between H = 17 and H = 23. Brighter objects are likely larger, making the Yarkovsky effect smaller and harder to detect. Fainter objects are less likely to be observed over a long time span, making them more difficult to discern.

The obliquity distribution outside of this magnitude range could be different. In particular, for near-Earth asteroids with known spin axes there is a larger fraction of mid-range obliquities for H< 18 than for H> 18. This difference could be caused by the 1 /D² dependence of the YORP effect (Vokrouhlický et al. 2015a), which makes it less effective at driving the spin to extreme obliquities on larger objects. On the other hand, the timescales of YORP cycles or stochastic YORP (Statler 2009) for objects smaller than the ones in our dataset can be shorter. Therefore, these objects would have a more frequent reconfiguration of the rotation state as the spin-up limit is reached.

Another limitation of our model is that we neglect non-principal-axis rotation. As discussed by Vokrouhlický et al. (2015b), for a tumbling asteroid, the Yarkovsky effect essentially depends on the rotational angular momentum rather than the spin axis. Therefore, if the fraction of non-principal-axis rotators were significant, our obliquity distribution would be more reflective of the orientation of the angular momentum vector than that of the spin axis.

The obliquity distributions presented here can be useful when an a priori distribution on the Yarkovsky perturbations is desired to model the trajectory of a target asteroid and no signal is yet visible from the orbital fit to the astrometry. This was the case for Apophis (Farnocchia et al. 2013b) and 2009 FD (Spoto et al. 2014) where the impact predictions are sensitive to the Yarkovsky effect though a direct estimate was not available.

Acknowledgments

The work of C. Tardioli was supported by the Marie Curie FP7-PEOPLE-2012-ITN Stardust, grant agreement 317185.D. Farnocchia and S. Chesley conducted this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. B. Rozitis acknowledges support from the Royal Astronomical Society (RAS) in the form of a research fellowship. T. Statler conducted some of this work while on an Intergovernmental Personnel Act appointment at NASA Headquarters; the results presented in this paper do not in any way constitute a statement of opinion, policy, or practice from NASA.

References

Appendix A: Additional table

All Tables

All Figures

	Fig. 1 Histogram of Yarkovsky estimates normalized by an absolute magnitude scale factor of . The average 1σ uncertainty in normalized A2 is 5 × 10-15 au/d2.
In the text

	Fig. 2 Probability distribution of the absolute magnitude for the objects in Table A.1.
In the text

	Fig. 3 Median (crosses) and 1σ standard deviation of the rotation periods from the JPL Small-Body Database for bins of 1 mag from H = 10 to H = 30.
In the text

	Fig. 4 Orbital period distribution as derived from the objects in Table A.1.
In the text

	Fig. 5 Above: enhancement factor α corresponding to 100% roughness as a function of the thermal parameter. The level curves correspond to a distance of 1 au. Below: histogram of the enhancement factor α.
In the text

	Fig. 6 Example of a three-bin distribution.
In the text

	Fig. 7 Example of a piecewise linear distribution.
In the text

	Fig. 8 Example of a quadratic distribution.
In the text

	Fig. 9 Discrete three-bin distribution. Thin lines indicate level curves of χ2 in the two-dimensional (2D) search space. The minimum χ2 is marked with a diamond. Thick lines mark loci of constant retrograde to direct rotator ratios. The dotted area corresponds to inadmissible values of the parameters.
In the text

	Fig. 10 Level curves of χ2 and ratio for the PL distribution model with two parameters. Lines and markers are the same as in Fig. 9.
In the text

	Fig. 11 Quadratic distribution: level curves of χ2 and ratio in the 2D search space. Lines and markers are the same as in Fig. 9.
In the text

	Fig. 12 Best-fit obliquity distributions for the three-bin, piecewise linear, quadratic, and four-bin models.
In the text

	Fig. 13 Comparison between the probability distributions in normalized A2 obtained from the Yarkovsky estimates and that from the different obliquity distribution models.
In the text

	Fig. 14 Comparison between the obliquity distribution presented by Farnocchia et al. (2013a) and the four-bin split of our the best-fit obliquity distributions converted to four bins. In the legend we recall the ratio between retrograde and direct rotators RR/D.
In the text