Constraints on the nearEarth asteroid obliquity distribution from the Yarkovsky effect
^{1} Department of Mechanical & Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK
email: Chiara.Tardioli@strath.ac.uk
^{2} Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
^{3} School of Physical Sciences, The Open University, Milton Keynes, MK7 6AA, UK
^{4} Department of Physics and Electronics, University of Puerto Rico at Humacao, 00792 Humacao, Puerto Rico
^{5} Department of Astronomy, University of Maryland, College Park, MD 20742, USA
^{6} Planetary Science Division, Science Mission Directorate, NASA Headquarters, Washington, DC 20546, USA
Received: 9 June 2017
Accepted: 3 September 2017
Aims. From light curve and radar data we know the spin axis of only 43 nearEarth asteroids. In this paper we attempt to constrain the spin axis obliquity distribution of nearEarth asteroids by leveraging the Yarkovsky effect and its dependence on an asteroid’s obliquity.
Methods. By modeling the physical parameters driving the Yarkovsky effect, we solve an inverse problem where we test different simple parametric obliquity distributions. Each distribution results in a predicted Yarkovsky effect distribution that we compare with a χ^{2} test to a dataset of 125 Yarkovsky estimates.
Results. We find different obliquity distributions that are statistically satisfactory. In particular, among the considered models, the bestfit solution is a quadratic function, which only depends on two parameters, favors extreme obliquities consistent with the expected outcomes from the YORP effect, has a 2:1 ratio between retrograde and direct rotators, which is in agreement with theoretical predictions, and is statistically consistent with the distribution of known spin axes of nearEarth asteroids.
Key words: methods: statistical / celestial mechanics / minor planets, asteroids: general
© ESO, 2017
1. Introduction
The complex motion of nearEarth asteroids (NEAs) is dominated by the gravitational perturbations of the Sun and planets. However, the gravitational interaction with other small bodies and nongravitational 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 nongravitational acceleration due to the anisotropic emission of thermal radiation of Solar System objects that causes a secular drift in the semimajor axis (Bottke et al. 2006). This perturbation is important to understand the longterm 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 highprecision 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 semimajor 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 wellobserved astrometric arc, it is possible to directly estimate the Yarkovsky effect by measuring deviations from a gravityonly 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^{1} and the DAMIT^{2} 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 fourbin 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 CottoFigueroa (2013).
2. Yarkovsky modeling
The Yarkovsky perturbation can be modeled as a transverse acceleration A_{2}/r^{2} (Farnocchia et al. 2013a), where r is the distance from the Sun in au and A_{2} 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^{2} 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 spinrate and thermalinertia dependence of the Yarkovsky acceleration for a Lambertian sphere, it is nonmonotonic 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_{1},F_{2} are positive functions that do not depend on the obliquity. To simplify we replace the instantaneous heliocentric distance with the solar fluxweighted mean heliocentric distance , where a and e are orbital semimajor 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_{2}, D, A, ρ, Γ, , P_{rot}, P_{rev}, and α.
3. Dataset of Yarkovsky estimates
To obtain a distribution of A_{2} on the lefthand 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 signaltonoise 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_{2} 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^{3} (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 nonextreme 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_{2} 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_{2} values by an absolute magnitude scale factor 1329 km10^{− H/ 5}/ with a constant albedo p_{v} = 0.154. The resulting A_{2} 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 A_{2} is 5 × 10^{15} au/d^{2}. 

Open with DEXTER 
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. 

Open with DEXTER 
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 Mclass and the Eclass 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).
Geometric albedo, bulk density and frequency for different taxonomic classes in an Hlimited sample.
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 SmallBody Database^{3} 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 lognormal distribution.
Rotation period. The rotation period P_{rot} is size dependent. In particular, there is the socalled 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 SmallBody 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 SmallBody Database for bins of 1 mag from H = 10 to H = 30. 

Open with DEXTER 
Orbital period and solarfluxweighted 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 semimajor axis a and e corresponding to our set of Yarkovsky detections: the semimajor 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 semimajor axis to the orbital period P_{rev} (see Fig. 4).
The distribution in solarfluxweighted 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. 

Open with DEXTER 
Enhancement factor. Smallscale surface roughness enhances the diurnal component of the Yarkovsky effect through thermalinfrared beaming, that is, reradiation of absorbed sunlight back towards the Sun (Rozitis & Green 2012). The degree of enhancement is a nonlinear 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 α. 

Open with DEXTER 
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_{2} 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 threebin distribution. 

Open with DEXTER 
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_{1}, Q_{2}, and Q_{3} in −1, 0, and 1, respectively. Since the integral has to be one, that is, Q_{1} + 2Q_{2} + Q_{3} = 2, the number of independent parameters is 2. This model can be generalized by having a variable abscissa x_{2} 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. 

Open with DEXTER 
Quadratic model. The final model we consider is that of a quadratic function f(γ) = acos^{2}γ + bcosγ + c (see Fig. 8). We allow only concaveup 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. 

Open with DEXTER 
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_{2} (see Fig. 1). Therefore, for each parametric obliquity distribution we obtain a predicted distribution in A_{2} to be compared with that coming from the set of Yarkovsky estimates. As already described in Sect. 3, the A_{2} 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 χ^{2} test. The range of normalized A_{2} 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, χ^{2} is computed as: (4)We look for the obliquity distributions providing a lower χ^{2}, which from standard χ^{2} 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 χ^{2}, we based our predicted distribution on 10^{5} 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 threebin 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 χ^{2} for each grid point. Figures 9−11 show the level curves of χ^{2} as a function of the parameters.
Fig. 9 Discrete threebin distribution. Thin lines indicate level curves of χ^{2} in the twodimensional (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. 

Open with DEXTER 
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. 

Open with DEXTER 
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. 

Open with DEXTER 
Table 2 gives the bestfit χ^{2}, the pvalue, that is, the probability of randomly obtaining a larger χ^{2}, and the corresponding ratio between retrograde and direct rotators. All the models provide statistically acceptable pvalues, with a quadratic model giving the best fit with a χ^{2} of 7.4.
Bestfit parameters, χ^{2}, pvalue, and retrograde to direct rotators ratio (R_{R/D}) for the models with two independent parameters.
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 selflimitation (CottoFigueroa et al. 2015) may provide a means to reconcile the high presentday retrograde fraction, where the YORPdriven deformation of aggregate objects confines their rotation rates to far narrower ranges than would be expected in the classical YORPcycle 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 χ^{2} decreases to 7. To find an absolute minimum value of χ^{2} on a threedimensional (3D) space we use the IDEA global optimizer (Vasile et al. 2011). For a fourbin distribution we find a best fit solution (p_{1},p_{2},p_{3},p_{4}) = (0.56,0.15,0.06,0.28) with minimum χ^{2} = 9.55 (pvalue of 22%) and R_{R/D} = 2.4. An Ftest shows that the Δχ^{2} = 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 midpoint abscissa the best fit solution is (Q_{1},Q_{2},Q_{3},x_{2}) = (1.16,0.04,0.73,0.074) with minimum χ^{2} = 8.14 (pvalue of 32%) and R_{R/D} = 1.8. Again, the statistical significance of Δχ^{2} = 10.43−8.14 as measured by an Ftest 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 χ^{2}. This suggests that our set of Yarkovsky estimates does not have enough signal to solve for more than two parameters.
Figure 12 shows the bestfit 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_{2} 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_{2} values, which is consistent with the preliminary results from CottoFigueroa (2013). The distributions from our models overestimate the heights of the peaks of the A_{2} 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.
Comparison of the bestfit obliquity distributions with the known spin axes from the EARN and DAMIT databases.
Fig. 12 Bestfit obliquity distributions for the threebin, piecewise linear, quadratic, and fourbin models. 

Open with DEXTER 
Fig. 13 Comparison between the probability distributions in normalized A_{2} obtained from the Yarkovsky estimates and that from the different obliquity distribution models. 

Open with DEXTER 
Figure 14 compares the Farnocchia et al. (2013a) result with the distributions found in this paper. For this purpose we convert our bestfit distribution to a fourbin 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^{4} and of 17 objects (15 already in the EARN database) from the DAMIT database^{5}. 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 midrange obliquities, which is again consistent with the YORP mechanism. It is interesting to point out that all of our bestfit 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 fourbin split of our the bestfit obliquity distributions converted to four bins. In the legend we recall the ratio between retrograde and direct rotators R_{R/D}. 

Open with DEXTER 
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 bestfit χ^{2} and pvalues for the different obliquity distributions obtained by setting α = 1. For all of the models, the low pvalues further support the presence of a Yarkovsky enhancement factor as predicted by Rozitis & Green (2012).
Bestfit χ^{2} and pvalues with enhancement factor α = 1.
8. Conclusions
We used the Chesley et al. (2016) list of Yarkovsky estimates to derive constraints on the obliquity distribution of nearEarth 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 χ^{2} 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. Ftests 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 bestfit is obtained with a quadratic model: f(cosγ) = 1.12cos^{2}γ−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 InteriorEarthObjects 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 nearEarth 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 nearEarth asteroids with known spin axes there is a larger fraction of midrange obliquities for H< 18 than for H> 18. This difference could be caused by the 1 /D^{2} 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 spinup limit is reached.
Another limitation of our model is that we neglect nonprincipalaxis 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 nonprincipalaxis 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 FP7PEOPLE2012ITN 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
 Benner, L. A. M., Busch, M. W., Giorgini, J. D., Taylor, P. A., & Margot, J.L. 2015, Radar Observations of NearEarth and MainBelt Asteroids, eds. P. Michel, F. E. DeMeo, & W. F. Bottke, 165 [Google Scholar]
 Bottke, W. F., Morbidelli, A., Jedicke, R., et al. 2002a, Icarus, 156, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Bottke, W. F., Vokrouhlickỳ, D., Rubincam, D. P., & Broz, M. 2002b, Asteroids III (Tucson: University of Arizona Press), 395 [Google Scholar]
 Bottke, W. F., Vokrouhlickỳ, D., Rubincam, D. P., & Nesvornỳ, D. 2006, Ann. Rev. Earth Planet. Sci., 34, 157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bowell, E., Hapke, B., Domingue, D., et al. 1989, in Asteroids II, eds. R. P. Binzel, T. Gehrels, & M. S. Matthews, 524 [Google Scholar]
 Busch, M. W., Giorgini, J. D., Ostro, S. J., et al. 2007, Icarus, 190, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Čapek, D., & Vokrouhlický, D. 2004, Icarus, 172, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Carry, B. 2012, Planet. Space Sci., 73, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Chesley, S. R. 2006, in Asteroids, Comets, Meteors, eds. L. Daniela, M. Sylvio Ferraz, & F. J. Angel, IAU Symp., 229, 215 [Google Scholar]
 Chesley, S. R., Chodas, P. W., Milani, A., Valsecchi, G. B., & Yeomans, D. K. 2002, Icarus, 159, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Chesley, S. R., Ostro, S. J., Vokrouhlický, D., et al. 2003, Science, 302, 1739 [NASA ADS] [CrossRef] [Google Scholar]
 Chesley, S. R., Farnocchia, D., Nolan, M. C., et al. 2014, Icarus, 235, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Chesley, S. R., Farnocchia, D., Pravec, P., & Vokrouhlický, D. 2016, in Asteroids: New Observations, New Models, eds. S. R. Chesley, A. Morbidelli, R. Jedicke, & D. Farnocchia, IAU Symp., 318, 250 [Google Scholar]
 CottoFigueroa, D. 2013, Ph.D. Thesis, Ohio University [Google Scholar]
 CottoFigueroa, D., Statler, T. S., Richardson, D. C., & Tanga, P. 2015, ApJ, 803, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Delbò, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, Asteroid Thermophysical Modeling, eds. P. Michel, F. E. DeMeo, & W. F. Bottke, 107 [Google Scholar]
 Ďurech, J., Kaasalainen, M., Warner, B. D., et al. 2009, A&A, 493, 291 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ďurech, J., Kaasalainen, M., Herald, D., et al. 2011, Icarus, 214, 652 [NASA ADS] [CrossRef] [Google Scholar]
 Farinella, P., Vokrouhlickỳ, D., & Hartmann, W. K. 1998, Icarus, 132, 378 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., & Chesley, S. 2014, Icarus, 229, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S., Vokrouhlickỳ, D., et al. 2013a, Icarus, 224, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S. R., Chodas, P. W., et al. 2013b, Icarus, 224, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S. R., Tholen, D. J., & Micheli, M. 2014, Celes. Mech. Dyn. Astron., 119, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Farnocchia, D., Chesley, S. R., Milani, A., Gronchi, G. F., & Chodas, P. W. 2015, Orbits, LongTerm Predictions, Impact Monitoring, eds. P. Michel, F. E. DeMeo, & W. F. Bottke, 815 [Google Scholar]
 Giorgini, J., Ostro, S., Benner, L., et al. 2002, Science, 296, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Giorgini, J. D., Benner, L. A., Ostro, S. J., Nolan, M. C., & Busch, M. W. 2008, Icarus, 193, 1 [NASA ADS] [CrossRef] [Google Scholar]
 La Spina, A., Paolicchi, P., Kryszczyńska, A., & Pravec, P. 2004, Nature, 428, 400 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Mommert, M., Farnocchia, D., Hora, J. L., et al. 2014a, ApJ, 789, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Mommert, M., Hora, J. L., Farnocchia, D., et al. 2014b, ApJ, 786, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., & Vokrouhlickỳ, D. 2003, Icarus, 163, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Nugent, C. R., Margot, J. L., Chesley, S. R., & Vokrouhlický, D. 2012, AJ, 144, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Pravec, P., & Harris, A. W. 2007, Icarus, 190, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Rozitis, B., & Green, S. F. 2012, MNRAS, 423, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Rubincam, D. P. 2000, Icarus, 148, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Spoto, F., Milani, A., Farnocchia, D., et al. 2014, A&A, 572, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Statler, T. S. 2009, Icarus, 202, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, C. A., Trilling, D. E., Emery, J. P., et al. 2011, AJ, 142, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Vasile, M., Minisci, E., & Locatelli, M. 2011, IEEE Trans. Evolutionary Computation, 15, 267 [CrossRef] [Google Scholar]
 Vokrouhlický, D., Chesley, S. R., & Matson, R. D. 2008, AJ, 135, 2336 [NASA ADS] [CrossRef] [Google Scholar]
 Vokrouhlický, D., Bottke, W. F., Chesley, S. R., Scheeres, D. J., & Statler, T. S. 2015a, The Yarkovsky and YORP Effects, eds. P. Michel, F. E. DeMeo, & W. F. Bottke, 509 [Google Scholar]
 Vokrouhlický, D., Farnocchia, D., Čapek, D., et al. 2015b, Icarus, 252, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Warner, B. D., Harris, A. W., & Pravec, P. 2009, Icarus, 202, 134 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Additional table
List of Yarkovsky estimates as of September 2016.
All Tables
Geometric albedo, bulk density and frequency for different taxonomic classes in an Hlimited sample.
Bestfit parameters, χ^{2}, pvalue, and retrograde to direct rotators ratio (R_{R/D}) for the models with two independent parameters.
Comparison of the bestfit obliquity distributions with the known spin axes from the EARN and DAMIT databases.
All Figures
Fig. 1 Histogram of Yarkovsky estimates normalized by an absolute magnitude scale factor of . The average 1σ uncertainty in normalized A_{2} is 5 × 10^{15} au/d^{2}. 

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

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

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

Open with DEXTER  
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 α. 

Open with DEXTER  
In the text 
Fig. 6 Example of a threebin distribution. 

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

Open with DEXTER  
In the text 
Fig. 8 Example of a quadratic distribution. 

Open with DEXTER  
In the text 
Fig. 9 Discrete threebin distribution. Thin lines indicate level curves of χ^{2} in the twodimensional (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. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 12 Bestfit obliquity distributions for the threebin, piecewise linear, quadratic, and fourbin models. 

Open with DEXTER  
In the text 
Fig. 13 Comparison between the probability distributions in normalized A_{2} obtained from the Yarkovsky estimates and that from the different obliquity distribution models. 

Open with DEXTER  
In the text 
Fig. 14 Comparison between the obliquity distribution presented by Farnocchia et al. (2013a) and the fourbin split of our the bestfit obliquity distributions converted to four bins. In the legend we recall the ratio between retrograde and direct rotators R_{R/D}. 

Open with DEXTER  
In the text 