© ESO, 2016

1. Introduction

The amount of photometric data of asteroids has been growing rapidly in recent years. These data are a source of information about shapes and rotational states of asteroids. Currently, the main method for determination of spin states and shapes of asteroids from photometry is the inversion of light curves, which was developed by Kaasalainen & Torppa (2001) and Kaasalainen et al. (2001). Models obtained with this method are stored in the Database of Asteroid Models from Inversion Techniques (DAMIT; Ďurech et al. 2010), which now contains models for 907 asteroids. The photometric data can be formally divided into two groups: (i) data dense in time that sample the rotational period well and that are typically used for the light curve inversion method, and (ii) data sparse in time (few measurements per night) that are produced by all-sky surveys, such as Pan-STARRS, Catalina, or LONEOS. Kaasalainen (2004) and Ďurech et al. (2005, 2007) showed that it is possible to get the solution of the inverse problem from sparse photometry if the data are of good quality (noise ≲5%). New asteroid models were also derived with a combination of dense and sparse photometry (Ďurech et al. 2009; Hanuš et al. 2011, 2013, 2016).

In the first statistical study of pole orientation of asteroids (based on 20 bodies), Magnusson (1986) revealed the lack of poles close to the ecliptic plane. That was later confirmed in analyses by Pravec et al. (2002), Skoglöv & Erikson (2002), and Kryszczyńska et al. (2007) for slightly less than 100 asteroids. Hanuš et al. (2011), using a sample of 206 main belt asteroids, found the dependence of the distribution of ecliptical latitudes β on the diameter D. They found basically isotropic distribution of β value with only a slight excess of prograde rotators for D ≳ 60 km, while the distribution of β value for D ≲ 30 km asteroids was found to have a strong preference for either low or high values indicating pole orientation near the pole of the ecliptic. The lack of poles near the ecliptic is most probably due to the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect, which can alter the direction of the spin axes of asteroids smaller than ≈40 km on a timescale shorter than their collisional lifetime (e.g., Pravec & Harris 2000; Rubincam 2000). The distribution of ecliptical longitudes λ of spin axes was, however, supposed to be rather uniform. For instance, Davis et al. (1989) came to this conclusion from the simulations of the collisional evolution of asteroids. With the growing number of asteroids for which pole orientation have been determined, a reliable statistics could be achieved and this hypothesis could be tested. However, even for a sample of 206 asteroids, Hanuš et al. (2011) did not reveal any nonuniformity in distribution of λ, but at the same time the data sample was too small to indicate meaningful nonuniformities. On the contrary, Slivan (2002) and Slivan et al. (2003) revealed a nonuniform pole distribution for 20−35 km size members in the Koronis family. In particular, the prograde-rotating asteroids all had ecliptic longitude between 24° and 73°. This conundrum was resolved by Vokrouhlický et al. (2003), who showed that these objects underwent a 2−3 Gyr long dynamical evolution during which the YORP effect tilted their spin axis near the ecliptic pole. Since YORP also continued to decrease the rotation frequency in their model, the spin state was captured in the Cassini resonance between the pole precession owing to solar torque and orbit precession from Jupiter-Saturn perturbations. The stationary point of this particular secular, spin-orbit resonance is currently at ≃35° ecliptic longitude. Thus all bodies whose spin axes librate about this point must have λ near this value. More recently, Bowell et al. (2014) estimated the ecliptical longitudes λ for more than 350 000 asteroids of the main belt using the magnitude method (Magnusson 1986), based on the variation of brightness with the ecliptical longitude: the maximum of brightness corresponds with the spin axis pointing either toward or opposite from the Earth. Surprisingly, the resulting distribution is clearly nonuniform with an excess of asteroids with λ from 30° to 110° and with minimum for 120° to 160°.

The success of getting a unique solution of the inverse problem with currently available sparse photometric data (which are not accurate enough) is low. Nevertheless, using the distributed computing project Asteroids@home (Ďurech et al. 2015), which significantly reduces the computational time of the period search, Ďurech et al. (2016) derived 328 new models from the analysis of Lowell photometric data. This is an impressive, but still small increase in number to enable a population-wide study. For this reason, we describe a new method for the determination of the orientations of spin axes and shapes of asteroids to utilize photometric data sparse in time. The uncertainties of spin vectors are large for individual bodies, therefore we work with groups of asteroids and construct distributions of tested parameters because working with large samples of bodies should smear uncertainties of individual solutions and, if uncorrelated, the results should hold in a statistical sense.

The structure of this paper is as follows: in Sect. 2, we describe our model and test its reliability on synthetic data; in Sect. 3, we apply the model to the photometric data from the Lowell Observatory database and construct the distributions of ecliptical longitudes for main-belt asteroids and for several groups of asteroids; Sect. 4 deals with distributions of the ratio a/b of axes of asteroids and, in Sect. 5, we summarize the main results.

2. Model

In the light curve inversion method, all parameters describing the rotational state (i.e., the rotational period and orientation of the spin axis), shape, and light scattering on the surface are fitted, and the unique sidereal rotational period P has to be determined. In the case of dense photometric data, we can substantially reduce the computational time necessary for the determination of P by only searching the interval around the value estimated from dense light curves. For sparse data, we usually do not have any estimate of P and we have to search the interval of all possible values, which is time consuming. Moreover, for the majority of asteroids we currently do not have sparse data that is accurate enough to derive a unique rotational period. Therefore, to fully utilize sparse photometry, we developed a new model, which does not allow us to determine the rotational period, but provides an approximate solution for the orientation of the spin axis and the shape parameters of the asteroid.

We model asteroids as geometrically scattering triaxial ellipsoids (a ≥ b ≥ c = 1) rotating about the shortest axis of the inertia tensor. The parameters of the model are the ecliptic longitude λ and latitude β of the pole and the ratios of axes a/b and b/c of the ellipsoid, alternatively axes a and b. The advantage of this model is that the brightness L, which is proportional to the projected area of the illuminated and visible part of the surface, can be computed analytically (Connelly & Ostro 1984) as follows: (1)where e, s are unit vectors defining the position of the Earth and the Sun in the asteroid coordinate system of principal axes of the inertia tensor, and (2)In a special case of opposition e = s, the Eq. (1) simplifies to (3)The direction toward Earth can be described by the rotational angle φ and aspect angle θ (i.e., angle between e and the direction of the spin axis), (4)Having set c = 1, the squared brightness L² normalized by the maximal possible value πab is (5)The mean quadratic brightness over one rotational period is then (6)and the normalized dispersion of squared brightness is (7)We used Eqs. (6) and (7), to compute and η_model for each asteroid and for each of its apparition; we defined apparitions as sets of observations with the gap between these sets of at least 100 days.

For the observational data, we used the following procedure:

• 1.

  We remove the dependence on solar phase angle. The changes in brightness in the light curve of an asteroid are not only due to the rotation but also the geometry of observation. In the model, we assume the case of opposition, which means the solar phase angle α = 0. For the observational data, we fitted the dependence of the brightness on the solar phase angle α by a linear-exponential dependence similar to Hanuš et al. (2011), i.e., (8)where g, h, d, k are parameters fitted for each asteroid, and we divided the observed brightness by that function. As an example, the corrected data for asteroid (511) Davida are shown in Fig. 1.

• 2.

  Then, we required that there were enough data for each asteroid: at least 20 points in one apparition and at least five apparitions for one asteroid (in Fig. 1 there are data from ten apparitions that can be used).

Computed theoretical values of the mean brightness ⟨L²⟩ and of the dispersion of the brightness η can be compared with observation by χ² that we define as: (9)where σ denotes the standard deviation and w denotes the weight for χ² of brightness. To normalize values of calculated and observed mean quadratic brightness, we divided them by and , respectively, which are mean values calculated over all apparitions. The value of w is not set in advance and has to be found by testing on known data. Since ⟨L²⟩ and η are not Gaussian random variables, the χ² in relation (9) is not χ² distributed. Nevertheless, we use this χ² formalism to define the best solution, which has the minimum χ².

Fig. 1 Photometric data of the asteroids (511) Davida corrected by the influence of the solar phase angle (black points). Red points with vertical lines denote the observed mean brightness and its dispersion in individual apparitions, green color denotes the same but calculated quantities for the best-fit model. To normalize L, we divided each value by mean value of L calculated over all apparitions.

In passing we note that in combining Eqs. (6) and (7) we obtain, for a given asteroid, the relation between ⟨L²⟩ and η, (10)This implies that for larger ⟨L²⟩ the model predicts smaller dispersion η. This is in accord with the intuition that a larger brightness corresponds to the pole-on geometry of view (i.e., smaller aspect angle θ).

To find a model with the best agreement (the lowest χ²) between the calculated values and the observation, we computed model values on a grid in parameter space: the ecliptical longitude of the pole from 0° to 360°; the latitude from 0° to 90°, both with a 5° step, and the axes a and b, from 1.1 to 4 and from 1 to a, respectively, both with a 0.1 step (an elongation larger than 4:1 would be unrealistic). As mentioned above, we corrected the observed brightness to the solar phase angle α = 0, however, the geometry remained unchanged and the aspect angle θ, which appears in Eqs. (6) and (7), was calculated for each apparition as a mean value as follows: (11)where v = [cosβcosλ,cosβsinλ,sinβ] ^T is the vector defining the direction of the spin axis and e_mean is the mean vector defining the position of the Earth during one apparition. From the relation (11) we can see that we obtain the same aspect angle for λ, β and λ ± 180° , −β, which is the reason why we test β only in the interval from 0° to 90°. Relation (11) also indicates that, for most asteroids, there is only a slightly worse second minimum of χ² for λ ± 180°. For zero inclination of orbit (e_z = 0), the aspect angle would be the same for λ and λ ± 180°. Owing to this ambiguity in λ, we constructed distributions of λ only in the interval 0°–180° and for λ> 180° we used modulo 180°.

2.1. Testing of the model on synthetic data

To test our model and confirm its reliability, we created synthetic data. We computed the brightness of asteroids using the models from DAMIT database and the Hapke scattering model (Hapke 1981, 1993) with randomly chosen parameters, and we assigned these new (synthetic) values to asteroids contained in the Lowell database (to the time of observation and the appropriate geometry). The distribution of poles for this synthetic data was isotropic.

We added the Gaussian noise (we tested noise σ_L = 0.15 and 0.2), which was then subtracted according the relation (12)if η ≥ 2σ_L, else η_obs = 0. For the real data, we only have an estimate of the noise level and we attempt to subtract different values from the data to find the best results. We also tested synthetic data without any noise (σ_L = 0).

After applying our model on these data, we should obtain uniform distributions of the ecliptical longitudes λ and latitudes sinβ. This was satisfied for the resulting distribution of λ, however, the distribution of latitudes showed a preference for high β. The possible explanation is that we did not include the uncertainties from the Hapke model and from the assumption that asteroids are triaxial ellipsoids. That means, for example, that for synthetic data without any noise and for an asteroid with β = 0, there are still some changes in brightness that our model interprets as nonzero β. To improve the model we added a new parameter that we called model noise σ_model. Then Eq. (12) had to be changed to (13)if , else η_obs = 0.

We tested values σ_model = 0.05, 0.06, 0.07 and 0.1. The resulting distributions of λ were uniform independently on σ_model. This is probably because λ is principally determined from the mean brightness ⟨L²⟩, which is comparatively more stable than the dispersion of brightness η from which β is determined. In the left panel of Fig. 2, there are shown distributions of sinβ for the two best values of σ_model and for the data noise σ_L = 0. The distributions are clearly nonuniform, nevertheless this is the best result we obtained. When we added noise to the synthetic data, we found that there is no significant difference between distributions of sinβ for σ_model = 0.06 and 0.07 (see Fig. 2 on the right) and, therefore, we decided to use the value 0.06 for the real data.

Fig. 2 Distributions of ecliptical latitudes β calculated from synthetic data with noise σL = 0 (left) and σL = 0.15 (right) for two best values of σmodel.

The takeaway message from our tests is that (i) determination of λ is reliable in a statistical sense, while (ii) determination of β is subject to systematic bias that needs to be corrected before interpreting the results.

3. The distribution of ecliptical longitudes

Having tested our approach and calibrated its parameters, we now construct the distribution of ecliptical longitudes for the real data from the Lowell Observatory photometric database (Bowell et al. 2014). This database contains data from 11 observatories, which are stored in the Minor Planet Center. The data were calibrated using the broadband accurate photometry of the Sloan Digital Sky Survey; the accuracy is ~0.1−0.2 mag. For more information about the data reduction and calibration, see Oszkiewicz et al. (2011).

First, we applied our model to 765 asteroids included in DAMIT database (from the first 10 000 numbered asteroids, which are included in the Lowell Observatory database and satisfy the conditions on the number of apparitions and the number of measurements in one apparition) and tried different values of noise σ_L (0.08, 0.1, 0.12, 0.15) and weight w (1, 5, 25); the value of model noise was 0.06. To decide on the best noise level and weight, we compared the calculated λ and β with λ_DAMIT (values from DAMIT derived with the light curve inversion) and β_DAMIT, respectively. From the distributions of Δλ = | λ − λ_DAMIT |, we found the best value of weight as w = 5 and from the distributions of Δβ we found the best value of noise level as σ_L = 0.08. However, we revealed that with this assumed data noise, the model produces hardly any spheroidal asteroids a/b ~ 1. This is because the photometric data for less bright asteroids have higher noise level than for brighter asteroids. In DAMIT, there are preferentially brighter asteroids, hence the noise level 0.08 works for them, but for less bright asteroids, such noise level is underestimated. To estimate the dependence of σ_L on L we used the amplitudes A_mag of light curves stored in the Asteroid Lightcurve Database (LCDB)¹ compiled by Warner et al. (2009). For A_mag we can write (14)where L is given by Eq. (5). The normalized dispersion of brightness η, defined by Eq. (7), is then related with amplitude as (15)where A = L_min/L_max = 10^{− 0.4A_mag}. For 9698 asteroids included in LCDB, we calculated η_A according Eq. (15) and then the appropriate noise level in data for each asteroid is written as: (16)if , else σ_L = 0. We calculated the running mean of σ_L for the sample of 500 bodies to obtain the dependence of σ_L on the mean brightness over all apparitions L_mean. The resulting dependence, with dispersion of σ_L among corresponding 500 bodies, is shown in Fig. 3. We applied this dependence in our model as follows: we assumed the noise level σ_L = 0.07 for asteroids with L_mean> 80; the brightness here is a dimensionless quantity calculated from magnitude M as L = 10^{− 0.4(M − 15)}. For asteroids less bright than 80, we calculated the noise level according to the equation of parabola, (17)The appropriate curve is shown in Fig. 3 (red line). We can see it does not fit the data perfectly, nevertheless, considering the dispersion of values of σ_L (gray lines), such deviation is insignificant.

Fig. 3 Dependence of the mean value of noise level σL on the mean value of Lmean; Lmean is mean brightness over all apparitions. Each mean value was calculated from a sample of 500 asteroids (as described in the main text). Gray lines denote dispersions of σL among corresponding 500 bodies. The red line denotes the fit that was applied in the model.

The Lowell Observatory database contains, in total, data for 326 266 asteroids. For 69 053 asteroids, there were enough apparitions and data points to calculate ecliptical longitude λ and latitude β; the vast majority of these asteroids belong to the first 100 000 numbered asteroids. For this sample, we used our model with weight w = 5, model noise 0.06, and data noise calculated for each asteroid according to the rule described above. The resulting distribution of the ecliptical longitude λ of asteroid pole orientation, shown in Fig. 4, is clearly nonuniform. As we can see, there is an excess of asteroids with λ from 40° to 100° and a minimum for λ ~ 150°. We calculated the Kolmogorov-Smirnov (KS) test of this distribution with a uniform one. The probability that they belong to the same parent distribution Q_KS is almost zero. A similar result was obtained by Bowell et al. (2014), who determined λ from the maximum of a sinusoid curve fitting the variation of brightness.

Fig. 4 Distribution of λ derived for 69 053 asteroids from the Lowell Observatory photometric database with model noise σmodel = 0.06 and weight w = 5. The values of λ larger than 180° map to the values λ − 180°, because of the symmetry of the model.

The distribution of ecliptical latitudes β shows strong preference for sinβ ≳ 0.9, however, since the distribution of sinβ for the synthetic data was not uniform (Fig. 2), the determined ecliptical latitudes are affected by biases and uncertainties that are not properly modeled here; therefore, in the following text we mainly study the distribution of ecliptical longitudes λ.

3.1. Searching for an explanation

Up to now, there is no satisfactory explanation of such nonuniformity in the distribution of ecliptical longitudes λ. We considered the observational and method biases described in Marciniak et al. (2015) and Santana-Ros et al. (2015), nevertheless, we found these do not influence our results; therefore, we searched for some other observational biases and geometrical and dynamical effects as well.

3.1.1. Galactic plane bias

First, we tested the influence of the measurements near Galactic plane, where the stellar background is more dense and thus the measurements may have higher uncertainties. We eliminated the observations with Galactic latitude | b | < 10° and repeated the analysis; for one asteroid there were on average about 6% less points. The differences between computed λ and λ from the DAMIT database were comparable with values for the model with the Galactic plane, however, the nonuniformity in λ was even larger. This result could suggest that, on the contrary, the shortage of observations near the Galactic plane could cause the nonuniformity of λ. However, if such a bias could influence our results, it would have also been seen in our test with synthetic data, since the geometry of observations was kept unchanged. Nevertheless, the resulting distribution of ecliptical longitudes was uniform, therefore, we believe our results are not influenced by such bias and we had to look for another explanation.

3.1.2. Correlation with longitude of ascending node

Next, we studied the role of the orbital longitude of node Ω by examining a possible correlation between asteroid’s pole longitude λ and Ω. The orbital data were taken from the AstOrb catalog². Figure 5 shows distribution of Ω values for 566 089 multi-opposition orbits of main-belt asteroids.

Focusing first on the data in the ecliptic reference system, we note that Ω values show overpopulation centered at ≃100° value, and underpopulation shifted by about 180°, i.e., centered at ≃270° value. This result is not new (see, e.g., JeongAhn & Malhotra 2014, and references therein). The reason for this nonuniformity in Ω is due to planetary perturbations. The distribution of Ω transformed to the Laplace plane shows similar nonuniformity, only shifted by ~180°; this is due to a slight but significant ≃1.58° tilt between the ecliptic plane and invariant plane of planets. For small-inclination orbits (i.e., whose proper inclination value is small), this effect becomes larger, as also shown in Fig. 5. Having learned about the nonuniformity of osculating nodal longitudes of asteroids in the main belt we should now examine, whether the nonuniform distribution of their pole longitudes λ is not a simple implication of the primary effect in nodes.

Fig. 5 Distribution of the longitudes of ascending node Ω (from AstOrb catalog) for the asteroid main belt. The red curve represents the distribution in the ecliptic plane; the black curve indicates the distribution in the invariant (Laplace) plane. The green and blue lines represent the distribution in the invariant plane for asteroids with the inclination of the orbit I< 10° and I< 5°, respectively.

First, we ran the following experiment. We divided the asteroid population according to their value of Ω to 18 equal bins (each 20° wide). We found the bin that contains the smallest number N of asteroids, and from all other bins we randomly selected N objects. That way, we had a sample of asteroids whose distribution of nodes was uniform. We examined distribution of rotation poles of this subsample, in particular the distribution of their λ values, and we found it is still nonuniform, resembling that in Fig. 4. The KS test of compatibility of the λ distributions obtained from our subsample and the whole sample of asteroids gave us a likelihood Q_KS ≃ 0.90 that they have the same parent distribution. We repeated our experiment several times, creating new subsamples, and obtained the same results. We also ran the same experiment in the Laplace reference system, but the choice of reference plane does not influence the results. These experiments suggest that nonuniformity in the distribution of orbital nodes does not play fundamental role in the nonuniform distribution of pole longitudes of asteroid spins.

Fig. 6 Distributions of ecliptical longitudes λ of poles for different longitudes of ascending node Ω.

Yet, we would expect some relation between Ω and λ should exist. For instance, plotting λ distributions for asteroids in each of the Ω bins described above, we obtained data shown in Fig. 6. The results indicate that in each of the bins of restricted Ω values, distribution of pole longitude peaks at ≃(Ω − 90°). This is actually understandable in the simplest model, in which the spin axis of each asteroid just uniformly precesses about the normal to its osculating orbit due to solar gravitational torque. We have quantitatively tested how much a simple geometrical effect of such precession could contribute to the observed dependence of ecliptical longitude λ on node Ω. To that goal we assumed the pole position in the moving orbital plane is set with the obliquity ϵ and we chose the inclination of orbit I and the longitude of ascending node Ω. This initial set up was changed several times, specifically, we tested values of inclination sinI = 0.10, 0.15, 0.30, values of node Ω = 10°, 170°, 250° and values of obliquity ϵ<I, ϵ>I. Assuming a simple regular precession, we randomly picked many values of longitude ϕ, uniform in 0° to 360°. We then transformed poles to the ecliptic system, determined appropriate λ and construct a model distribution of the ecliptical longitudes. Results of these simple simulations satisfied our hypothesis of geometrical effect; for ϵ<I, the distribution of λ was only a tight interval of values near ≃(Ω − 90°), and for ϵ>I, the λ values ranged the whole interval from 0° to 360°, but with a peak at ≃(Ω − 90°). However, when we summed distributions of λ for values of ϵ from assumed distribution n(cosϵ) for a fixed I and Ω, we reached an almost uniform final distribution, which is far from the distributions shown in Fig. 6. We tested n(cosϵ) uniform and also some unrealistic distributions, for example, we assumed there were ten times more bodies with ϵ< 45° than with ϵ> 45°, but with insignificant effect on the final distribution.

Therefore, using two lines of evidence we show that the nonuniformity of the ecliptic Ω values together with only simple geometric (projection) effects cannot explain the nonuniformity in the distribution of pole ecliptic longitudes. However, the flow of pole orientation in the orbit frame may be much more complicated than just a simple steady precession about the orbital angular momentum vector. This is because of a possibility of resonant, spin-orbit effects described by Cassini dynamics (e.g., Colombo 1966; Henrard & Murigande 1987; Vokrouhlický et al. 2006). In fact, the large-asteroid subgroup in the Koronis family, the Slivan sample, has actually been identified as being captured in the most prominent s₆ Cassini resonance resulting in a common orientation of their pole longitudes near the stationary point at ecliptic longitude ≃35° (e.g., Vokrouhlický et al. 2003). Therefore, we examined whether such resonant effects could help us to explain the nonuniformity in the λ distribution.

However, we found the answer is negative. First, if the capture in the aforementioned Cassini resonance played a dominant role population wise, the pole longitude distribution would be peaked at the stationary point of the resonance (shifted by some 35°−40° from the maximum seen in Fig. 4). Next, Vraštil & Vokrouhlický (2015) have shown that the capture in this resonance is generally unstable (especially in the inner part of the main belt), and that its phase volume is small (few percent at maximum). The latter implies that expecting the spin pole located in this resonance by chance is very small. In order to verify these preliminary conclusions, we used the software described in Vraštil & Vokrouhlický (2015) to probe the expected effect. This is basically much more sophisticated variant of our previous Monte Carlo experiment in which we assumed a steady precession in the orbit frame. Here we propagated orbit and spin evolution of the first 10 000 main belt asteroids, giving them random initial rotation state parameters, such as rotation period, pole orientation, and dynamical ellipticity. We then numerically propagated orbit and spin evolution for tens of millions of years and monitored distribution of simulated ecliptic longitudes of the sample. We found the sample quickly forgets given initial conditions and fluctuates about a steady-state situation with basically uniform distribution of ecliptic longitudes of rotation poles. We repeated the numerical experiment several times with different initial conditions but always obtained very similar results.

3.2. Distributions of λ for groups of asteroids

Our next step was to study the distributions of λ for various groups of asteroids, specifically for asteroids with different sizes, different spectral types, dynamical families, and asteroids in different parts of the main belt. Distributions were again compared using the KS test.

3.2.1. Asteroids with different sizes

We divided asteroids into eight groups according their diameters: 0–3; 3–6; 6–9; 9–12; 12–15; 15–25; 25–50; and 50–1000 km; the number of asteroids decrease with higher diameters, therefore, we chose wider ranges of bins. We preferentially used diameters derived from the observations of the WISE satellite (Masiero et al. 2011)³. For asteroids not included there, we used diameters from AstOrb catalog. We compared distributions with each other and found that the differences are not significant, which means that the data do not reveal any dependence of λ on size. We also studied the dependence of ecliptical latitude β on size and tried to confirm the result from Hanuš et al. (2011). In Fig. 7, we can see that, even though the distributions show preference for sinβ ≳ 0.9, with decreasing diameter D, there is a visible depopulation of spin axes close to the ecliptic plane, which is in agreement with findings of Hanuš et al. (2011).

Fig. 7 Distribution of ecliptical latitudes β for groups of asteroids with different sizes. All distributions are divided by the distribution for synthetic data with noise 0.15 and model noise 0.06 (Fig. 2 on the right, black curve). This is the simplest debiasing procedure we can use.

This is yet another interesting hint about the origin of the nonuniformity of λ distribution. The affinity of latitudes toward extreme values for small asteroids clearly shows that the YORP effect has been affecting the population in a size selective way exactly predicted by the theory. However, the distribution of the longitudes does not indicate this size-selectivity, implying the YORP effect is not the primary mechanism in the longitude story. Indeed, the theory of the YORP effect so far has not predicted any significant effects for the pole longitude.

3.2.2. Taxonomic classes

We compared distributions of λ between asteroids belonging to the taxonomic class C and S (using the AstOrb catalog), which are the largest groups. The result of KS test, Q_KS = 0.45, indicates that there is no significant difference. The fact that the distribution of λ is independent of the sizes and taxonomic classes simplified tests with other subpopulations.

3.2.3. Different parts of the main belt

We also studied the distributions of λ for groups of asteroids located in different parts of the main belt. Specifically, asteroids with different inclinations sinI, eccentricities e, and semimajor axes a of their orbits. We found that the distribution of λ is not dependent on the eccentricity, however it is strongly dependent on the inclination (see Fig. 8). For sinI< 0.02 there is a huge excess of asteroids with λ from 60° to 100°, there are more than four times more bodies than for λ ~ 150°. With increasing I the distributions are closer to the uniform distribution. This result is surprising and it actually goes against the ideas about simple geometrical (projection) effects discussed in Sect. 3.1.2, suggesting that perhaps some unidentified yet dynamical effects are at play.

We also studied the dependence of the distribution of λ on the inclination of orbit in the invariant plane. Although the maximum of distribution for sinI< 0.02 is slightly lower, there is still strong dependence on the inclination.

Fig. 8 Distributions of ecliptical longitudes for asteroids with different inclinations sinI of their orbits.

We constructed distributions of λ for individual Ω bins (as in Fig. 6) for asteroids with sinI< 0.04. The peaks of all distribution were for λ from 60° to 100°, which corresponds with the distribution of λ for small inclinations. This means that for orbits with small inclination, the dependence of λ on Ω is suppressed.

We then constructed distributions for asteroids with different semimajor axes a. We separated the main belt into four parts (see Fig. 9), the inner, middle, pristine⁴, and outer belt, which are separated by mean-motion resonances with Jupiter. To eliminate the dependence on the inclination of orbit we divided asteroids of each part into bins with different inclination (we used the same bins as in Fig. 8) and we randomly chose such number of asteroids to have the same number of asteroids in corresponding bins of two populations. In other words, the distributions of inclination of orbit for the compared populations were the same. The results of KS tests show that only the pristine zone, bracketed by the powerful mean motion resonances 5/2 and 7/3 with Jupiter at ≃2.82 au and ≃2.96 au, has significantly different distribution from the middle and outer belt (Q_KS< 3 × 10^-6); specifically, the nonuniformity is more significant in pristine zone than in other parts. For the pair inner belt and pristine zone, the KS test gives Q_KS = 0.00013.

Fig. 9 Four parts of the main asteroid belt defined according to the proper semimajor axis a (we used proper values of a and I from Asteroids Dynamic Site; Knežević & Milani 2003). The number of objects plotted here is 100 000. The horizontal line separates the region sin I< 0.04, for which the maximum of distributions of λ for λ ~ 80° is highest. The locations of some more populous asteroid families are emphasized.

3.2.4. Dynamical families

Finally, we studied dynamical families. The family membership of asteroids was taken from Nesvorný et al. (2015). Distributions of individual families were compared with the distribution of corresponding background formed with asteroids from the same part (inner, middle, pristine, and outer) as the family and with inclinations of orbit from the interval defined by the members of the family. The KS test did not reveal any significant difference between any family and its background. We also compared families located approximately in the same interval of inclination (see Fig. 9) with each other, specifically: Themis with Massalia; Vesta with Eos, Hygiea and Flora; Hygiea with Flora; and Koronis with Nysa Polana. Again, the KS test showed no difference for these pairs of families. The distributions of λ for six selected families are shown in Fig. 10. The differences we can see between the distributions are caused only due to the dependence on the inclination of orbit.

Fig. 10 Distributions of λ for Flora, Nysa Polana, and Massalia (left). The distributions of λ for Eos, Themis, and Koronis (right).

The strong maximum of the λ distribution between ≃(80°−110°) in the Koronis family does not fit the interval of expected librators in Slivan states described by Vokrouhlický et al. (2003) which would be shifted by about 40° to 50° degrees toward lower values.

3.3. The bootstrap method

Formally, it is always possible to find the best ecliptical longitude λ and latitude β of the pole, i.e., the lowest χ². However, the minimum can be flat and in that case λ is not well determined. To estimate the errors of determined longitudes we applied the bootstrap method (Davison & Hinkley 1997) on the set of measurements for each asteroid; we used the first 10 000 numbered asteroids from the Lowell Observatory database, of which for 9774 there were enough data points. From the set we randomly selected data to get the same number of measurements, but some of them were chosen more than once and some of them were missing. We repeated this ten times, therefore, we obtained ten modified sets of measurements and thus ten possible longitudes for each asteroid. We considered that the longitude was well determined when the maximum difference among ten values of λ was ≤50°. This was satisfied for 3930 from 9774 asteroids; the mean value of the largest differences for these bodies is 30°. The dependences of λ on the longitude of ascending node Ω and on the inclination of orbit I for this new sample of 3930 asteroids did not significantly change, which means that the poorly constrained models did not cause any systematic effect to distribution of λ.

4. Distributions of the ratio of axes a/b

In this section, we study shapes of asteroids (specifically the ratios of axes a/b and b/c) derived from our model. We tested our model on synthetic data as described above in Sect. 2.1, the assumed noise was 0.15. The values of ratios a/b and b/c obtained with our model are compared with values from DAMIT models derived from the principal moments I₁, I₂, I₃ of the inertia tensor (assuming uniform density), (18)Since the values of ratios computed from our model were obtained from synthetic data based on DAMIT, they should be the same as values derived from the inertia tensor. The result is shown in Fig. 11. We calculated the linear (Pearson) correlation and Spearman correlation for both ratios, the coefficients ρ are summarized in Table 1. We obtained a good correlation for the ratio a/b, while the ratio b/c is not so well determined.

Fig. 11 Comparison of the values of a/b and b/c from DAMIT database with values calculated from synthetic data based on the models from DAMIT. One model was used more than once. The error bars are calculated as the dispersion from the mean value for one model from DAMIT.

For the real data, the setup was the same as described in Sect. 3: weight w = 5; model noise σ_model = 0.06; and data noise σ_L = 0.07, respectively, σ_L was calculated according Eq. (17) for asteroids less bright than 80. We compared resulting ratios a/b and b/c of 765 asteroids included in DAMIT with a/b_DAMIT and b/c_DAMIT and calculated correlation coefficients (see Table 1 and also Fig. 12). The correlation coefficients for the ratio b/c are lower than 0.1, which implies that b/c is not well determined and in following tests we will study only the ratio a/b. The problem to determine the ratio b/c is linked with our previous result that the distribution of ecliptical latitudes β, especially for small bodies (see Fig. 7), shows a preference for high values of β because for a spin axis with high latitude (small obliquity) we have observations only from limited range of polar aspect angles. The determination of b/c, however, requires observations from wide range of aspect angles.

Fig. 12 Comparison of the values of ratio a/b calculated from the real data from the Lowell database with values from DAMIT database.

As in Sect. 3.3 we used the bootstrap method to estimate errors of the ratio a/b. The allowed maximum difference among ten calculated values of a/b was 0.25, 3819 remain from 9774 asteroids, and the mean value of the largest differences is 0.18.

4.1. Distributions of a/b for asteroids with different sizes

As in the case of the ecliptic longitude λ, we studied distributions of a/b for several groups of asteroids. The test for asteroids with different diameters showed that larger asteroids (D> 25 km) are more spheroidal (values of a/b closer to 1) and smaller asteroids are more elongated (a/b ~ 1.6) as is shown in Fig. 13. The differences between distributions in Fig. 13 on the right are much bigger than the uncertainties estimated from bootstrap method 0.18.

This dependence of a/b on D also remained for a smaller sample of 3570 asteroids, which were obtained from the bootstrap method as having well-determined a/b. Since such dependence on diameter can influence the comparison of distributions of a/b of other populations of asteroids, we have to eliminate it in the following tests.

Fig. 13 Distributions of a/b for groups of asteroids with different sizes (from 0 to 12 km on the left; from 12 to 1000 km on the right).

McNeill et al. (2016) determined an average axial ratio for asteroids with diameter D< 8 km from PanSTARRS 1 survey as 1:0.85, i.e., a/b = 1.18, which is not in agreement with our findings. For the corresponding range of diameters, we calculated the average value of a/b for asteroids in DAMIT (derived from the principal moments of the inertia tensor). The average value is a/b_DAMIT = 1.47, however, the sample of asteroids from DAMIT with D< 8 km is biased. The DAMIT sample contains preferentially more elongated asteroids, for which is easier to find the solution of the light curve inversion method. We also checked our result, that larger asteroids are more often spheroidal, against the asteroids in DAMIT; for D> 50 km we obtained the average value of a/b_DAMIT = 1.23 (our model gives a/b = 1.29). We can conclude that the dependence of a/b on diameter is real, however our model gives higher values of a/b for asteroids with D< 25 km. Considering the dispersion of values σ_L in Fig. 3, this could be because of the underestimated data noise for smaller and less bright asteroids.

4.2. Different parts of the main belt

Next, we studied distribution of a/b for asteroids with different inclinations of their orbits. To remove the dependence of a/b on diameter, we used only asteroids with D< 20 km. The differences between resulting distributions of a/b are not so distinct as when we studied the dependence on diameter, and they are comparable with the uncertainties in a/b.

We also compared distributions of a/b for asteroids with different semimajor axes, specifically inner, middle, pristine, and outer belt, using only asteroids with diameters D< 20 km. The differences between distributions of a/b are not significant and are again comparable with the uncertainties in a/b.

4.3. Dynamical families and taxonomic classes

As in Sect. 3.2.4, we compared dynamical families with their backgrounds, using again only asteroids with D< 20 km. We did not reveal any significant differences between distributions of a/b of families and corresponding backgrounds. Also the comparison of families with each other did not show any differences larger than uncertainties in a/b.

Szabó & Kiss (2008) derived distributions of a/b for eight asteroids families using data from the Sloan Digital Sky Survey (SDSS). However, their distributions are different from ours; they are often bimodal (Figs. 4−6 therein) and the maximum is for a/b ~ 1.2 (our distributions have maximum around 1.6). They also suggest a possible dependence on the age of families (old families contain more spheroidal members), but we do not observe that in our distributions. We believe that they used an assumption that could influence the results. First they assumed that the rotational axes of all asteroids are perpendicular to the line of sight. Then they also tested fixed value β = 50° for all asteroids.

The last populations of asteroids we compared were different taxonomic classes, specifically C and S types. The result of KS test, Q_KS = 0.17, did not show any difference between these two groups.

5. Conclusions

We developed a new method that allows us to determine the orientation of rotational axes and equatorial axes ratio a/b, assuming a triaxial shape model, using sparse data obtained by all sky surveys. The goal of our approach is to provide a distribution function of the solved-for parameters for a large sample of main belt asteroids rather than detailed rotational state of individual objects. A limitation of our method is that it provides, first, longitude λ of the rotation pole in the interval (0°,180°) only, with values in (180°,360°) transformed to (0°,180°) by λ = λ − 180° rule; and second, absolute value of the ecliptic latitude β, instead of β itself. The latter means that we cannot determine whether the asteroid rotates in a prograde or retrograde sense. Additionally, our model also does not provide rotational period.

We first justified our method by applying it to a synthetic sample of asteroids and also to a known sample of objects with rotational states resolved by more detailed methods that can be found in the DAMIT database. We observed that our method reproduces well the distribution of the ecliptic longitudes and the equatorial axes ratio a/b in a statistical sense. The uncertainties, estimated using the bootstrap method, are ≃30° in λ and ≃0.2 in a/b without systematic effects on the mean value. The determination of ecliptical latitudes β shows bias toward finding preferentially spin axes near the ecliptic pole. Our tests may, however, provide a rough approximation of the bias function.

We then applied our method to 69 053 main belt asteroids for which a suitably rich and good quality set of observations were obtained from the Lowell Observatory database. The main results are as follows:

• 1.

  The distribution of λ is nonuniform, with an excess of asteroids with λ values between 60° and 100°. Similarly, there is a deficiency of asteroids with λ values between 130° and 160°. Curiously, our tests revealed a correlation of this nonuniformity with orbital inclination: asteroids with very low-inclination orbits (sinI ≤ 0.04) show the effect more significantly than asteroids with higher inclination orbits.

• 2.

  While not a primary result from our paper, we also determined distribution of the absolute value of sine of ecliptic latitude | sinβ |. We confirm previously reported results that asteroids with size D ≤ 25 km have their pole latitude tightly clustered about the poles of ecliptic. This is due to the YORP effect that makes the pole latitude to approach the extreme values asymptotically.

• 3.

  We also found that small main belt asteroids (D ≤ 25 km) are more elongated, with a median of ratio a/b ≃ 1.6, compared to the large asteroids (D ≥ 50 km), which have a median of ratio a/b ≃ 1.3.

• 4.

  We also analyzed our results for populations in different asteroid families. As to the λ distribution, they mainly derive from their inclination value of the aforementioned inclination dependence. For instance, the low-inclination families such as Massalia or Themis have the strongest nonuniformity of the λ distribution in our results.

Using a more detailed method, we confirmed the previously reported unexpected nonuniformity in distribution of ecliptic longitude of spin axes of the main belt asteroids. We tested various hypotheses of its origin, but we had to reject them, proving that the proposed processes would not lead to a significant enough nonuniformity. Therefore, this result remains enigmatic and requires further analysis. In particular, it would be very useful if more detailed methods of spin state and shape inversion from astronomical data confirmed this result and provided more details. We note, for instance, that methods both in Bowell et al. (2014) and here are not able to discriminate between the prograde- and retrograde-rotating asteroids. It would be important to see, if the excess in λ values at about 80° concerns equally well both classes, or whether it is preferentially associated with one of them. This could hint about the underlying processes that cause the effect. In the same way, all methods used so far fold the whole range of ecliptic λ values to a restricted interval (0°,180°).This is because of their intrinsic drawback of not distinguishing data for λ and λ + 180° cases. Yet, breaking this uncertainty may also help to disentangle the underlying physical causes of the nonuniformity.

Justifications of reliability of our method, by running blind tests against synthetic populations of asteroids and limited datasets for which complete models are already available, make our method a solid tool for further studies. It would be interesting to apply it to more accurate photometric data provided by Large Synoptic Survey Telescope (LSST).

Acknowledgments

H. Cibulková and J. Ďurech were supported by the grant 15-04816S of the Czech Science Foundation. We are grateful to an anonymous referee for useful and constructive comments.

References

All Tables

All Figures

Fig. 1 Photometric data of the asteroids (511) Davida corrected by the influence of the solar phase angle (black points). Red points with vertical lines denote the observed mean brightness and its dispersion in individual apparitions, green color denotes the same but calculated quantities for the best-fit model. To normalize L, we divided each value by mean value of L calculated over all apparitions.

In the text

	Fig. 2 Distributions of ecliptical latitudes β calculated from synthetic data with noise σL = 0 (left) and σL = 0.15 (right) for two best values of σmodel.
In the text

	Fig. 3 Dependence of the mean value of noise level σL on the mean value of Lmean; Lmean is mean brightness over all apparitions. Each mean value was calculated from a sample of 500 asteroids (as described in the main text). Gray lines denote dispersions of σL among corresponding 500 bodies. The red line denotes the fit that was applied in the model.
In the text

	Fig. 4 Distribution of λ derived for 69 053 asteroids from the Lowell Observatory photometric database with model noise σmodel = 0.06 and weight w = 5. The values of λ larger than 180° map to the values λ − 180°, because of the symmetry of the model.
In the text

Fig. 5 Distribution of the longitudes of ascending node Ω (from AstOrb catalog) for the asteroid main belt. The red curve represents the distribution in the ecliptic plane; the black curve indicates the distribution in the invariant (Laplace) plane. The green and blue lines represent the distribution in the invariant plane for asteroids with the inclination of the orbit I< 10° and I< 5°, respectively.

In the text

	Fig. 6 Distributions of ecliptical longitudes λ of poles for different longitudes of ascending node Ω.
In the text

	Fig. 7 Distribution of ecliptical latitudes β for groups of asteroids with different sizes. All distributions are divided by the distribution for synthetic data with noise 0.15 and model noise 0.06 (Fig. 2 on the right, black curve). This is the simplest debiasing procedure we can use.
In the text

	Fig. 8 Distributions of ecliptical longitudes for asteroids with different inclinations sinI of their orbits.
In the text

Fig. 9 Four parts of the main asteroid belt defined according to the proper semimajor axis a (we used proper values of a and I from Asteroids Dynamic Site; Knežević & Milani 2003). The number of objects plotted here is 100 000. The horizontal line separates the region sin I< 0.04, for which the maximum of distributions of λ for λ ~ 80° is highest. The locations of some more populous asteroid families are emphasized.

In the text

	Fig. 10 Distributions of λ for Flora, Nysa Polana, and Massalia (left). The distributions of λ for Eos, Themis, and Koronis (right).
In the text

	Fig. 11 Comparison of the values of a/b and b/c from DAMIT database with values calculated from synthetic data based on the models from DAMIT. One model was used more than once. The error bars are calculated as the dispersion from the mean value for one model from DAMIT.
In the text

	Fig. 12 Comparison of the values of ratio a/b calculated from the real data from the Lowell database with values from DAMIT database.
In the text

	Fig. 13 Distributions of a/b for groups of asteroids with different sizes (from 0 to 12 km on the left; from 12 to 1000 km on the right).
In the text