Reconstruction of asteroid spin states from Gaia DR3 photometry

Gaia Data Release 3 contains accurate photometric observations of more than 150,000 asteroids covering a time interval of 34 months. With a total of about 3,000,000 measurements, a typical number of observations per asteroid ranges from a few to several tens. We aimed to reconstruct the spin states and shapes of asteroids from this dataset. We computed the viewing and illumination geometry for each individual observation and used the light curve inversion method to find the best-fit asteroid model, which was parameterized by the sidereal rotation period, the spin axis direction, and a low-resolution convex shape. To find the best-fit model, we ran the inversion for tens of thousands of trial periods on interval 2-10,000 h, with tens of initial pole directions. To find the correct rotation period, we also used a triaxial ellipsoid model for the shape approximation. In most cases the number of data points was insufficient to uniquely determine the rotation period. However, for about 8600 asteroids we were able to determine the spin state uniquely together with a low-resolution convex shape model. This large sample of new asteroid models enables us to study the spin distribution in the asteroid population. The distribution of spins confirms previous findings that (i) small asteroids have poles clustered toward ecliptic poles, likely because of the YORP-induced spin evolution, (ii) asteroid migration due to the Yarkovsky effect depends on the spin orientation, and (iii) members of asteroid families have the sense of rotation correlated with their proper semimajor axis: over the age of the family, orbits of prograde rotators evolved, due to the Yarkovsky effect, to larger semimajor axes, while those of retrograde rotators drifted in the opposite direction.


Introduction
The first photometric measurements of asteroids from the ESA Gaia mission (Gaia Collaboration et al. 2016) were released in April 2018 as part of the Gaia Data Release 2 (DR2, Gaia Collaboration et al. 2018a). That dataset contained astrometry and photometry for ∼14,000 asteroids (Gaia Collaboration et al. 2018b) covering a time interval of about 22 months. The recent Gaia Data Release 3 from June 2022 (DR3, Tanga et al. 2022;Babusiaux et al. 2022) provided a significantly larger number of asteroid measurements: more than 3,000,000 photometric data points for about 150,000 asteroids covering a time interval of about 34 months. Typically there are fewer than 20 individual measurements per asteroid in DR3. The largest number of measurements is 96.
Although the Gaia DR2 photometric dataset was rather limited, we successfully derived physical models for almost 200 asteroids, consisting of mostly new solutions, and we published the results inĎurech & . We used the standard convex inversion method ) based on the inversion of photometric measurements. This physical model included the sidereal rotation period, the orientation of the spin axis, and a convex 3D shape model.
Thanks to a significantly larger amount of data and more extended temporal coverage of the DR3 data compared to DR2, we were expecting a dramatic increase in the number of successful shape model determinations leading to an unprecedented insight into the distribution of physical properties of asteroids. This paper describes our application of the same procedure as inĎurech  to process and analyze the DR3 data.
In Sect. 2 we describe the DR3 data downloading and processing and the inversion technique we applied, including the verification tests. In Sect. 3 we present our analysis of the spin properties of asteroids. We conclude our work in Sect. 4.

Inversion of Gaia asteroid photometry
We downloaded the data from the Gaia archive. 1 We selected only the relevant parameters: the Gaia-centric JD in TCB (epoch), the calibrated G-band magnitude (g_mag), the G flux (g_flux), and the error in the G flux (g_flux_error). Because the JD epoch is given for each CCD position but the magnitude is the same over the transit, we averaged the epoch values over the transit. 2 We converted JD from the original Gaia-centric TCB to TDB according to the formula given in Tanga et al. (2022). We also computed the relative flux error as g_flux_error / g_flux. In this way we obtained 3,069,170 photometric data points for 156,789 asteroids (Tanga et al. 2022 report 156,801 in their Table 1 because some asteroids have only NULL magnitudes). Then we computed the heliocentric and Gaia-centric ecliptic coordinates for each observation using the API interface of the JPL Horizons ephemeris service. 3 We computed the lighttime correction, corrected the brightness to 1 au distance from the Sun and Gaia, and converted magnitudes to relative flux.
We limited our sample to 60,945 asteroids that had 21 or more observations. For asteroids with fewer data points, the number of model parameters would be higher than the number of observations, which would often lead to unrealistic fits with zero residuals. The distribution of the number of observations per asteroid is shown in Fig. 1; the mean number of observations is 20, and the median is 17. There are about 2000 asteroids with more than 50 data points, 473 with more than 60, and only 89 with more than 70.
We performed a period search with a convex shape model parameterized by spherical harmonics on the order and degree of three. The shape was discretized with 288 surface elements of areas σ i with normals n i isotropically distributed in spherical coordinates ϑ and ϕ. Instead of directly optimizing σ i , we used the expansion into spherical harmonics series with associated Legendre polynomials P m where a lm and b lm were subject to optimization. The parameter a 0 corresponds to the size of the shape model. Because we do not know the albedo, the size cannot be determined and we set a 0 to some arbitrary fixed value and only the remaining 15 a lm , b lm parameters are optimized. Periodograms were computed with the errors of individual photometric measurements taken into account. To avoid errors that are too small and that are below the resolution of the convex-shape model, we set up the minimum relative error to 0.01. An asteroid's rotation state was described by its spin axis direction in ecliptic coordinates (λ, β) and the sidereal rotation period P. The last parameter of our model was the slope k of the phase curve. For the light scattering model we used a combination of Lambert S L = µµ 0 and Lommel-Seeliger S LS = µµ 0 /(µ + µ 0 ) multiplied by a phase function f (α), µ 0 and µ being cosines of the angles of incidence and reflection, respectively. So the scattering model was where c was fixed at 0.1 and the phase function had the form (Kaasalainen et al. 2002) The distribution of solar phase angles for all DR3 asteroid observations is shown in Fig. 2. Most of the observations were carried out at phase angles between 10 and 30 deg where a linear function can approximate the phase curve. Therefore, we optimized only the linear parameter k in eq. (3), while the other two parameters describing the increase in brightness near opposition were fixed at typical values a 0 = 0.5, d = 0.1. The values of k for the sample of reconstructed models (Sect. 2.5) were distributed between −1.5 and 0 with the mean value −0.88 and standard deviation 0.20. The shape of the distribution was close to the Gaussian distribution. The periodograms were computed on an interval of 2-10,000 h. The step size in the period was set to 0.8 ∆P, where ∆P = 0.5 P 2 /T is a typical distance between local minima in the period (Kaasalainen 2004), which means that the time step increases quadratically with increasing period. In frequencies this corresponds to a uniform sampling with the step 0.4/T . For the longest Gaia DR3 dataset with T ≈ 34 months, ∆P ≈ 6 × 10 −5 h for P = 2 h and ∆P ≈ 1600 h for the longest period of 10,000 h. Histogram showing number distributions of observations per asteroid in the whole DR3 dataset (blue) and for asteroids for which we derived a spin and shape model (orange). The black curve shows the success rate of deriving a model for each data bin with at least ten counts, which is the ratio of the orange to blue bins.
The total number of period steps on the whole interval was tens of thousands. A globally unique period solution was defined as having the lowest χ 2 min with all other periods giving χ 2 higher than χ 2 tr = (1 + √ 2/ν) χ 2 min , where ν is the number of degrees of freedom, and the root mean square (RMS) residuals of all local minima had to be higher than 0.01. We obtained a sample of 14,192 asteroids with a unique period defined this way. Then we performed a pole search with the same χ 2 limit and selected only asteroids with one or two pole solutions, which reduced the sample to 11,854 asteroids. Starting from these poles, we created final models and reconstructed their 3D shape by the Minkowski procedure (Lamberg & Kaasalainen 2001). We selected only asteroids for which this conversion was successful and for which the ratio of their moment of inertia along the principal axis to that along the actual rotation axis was less than 1.1 (Ďurech et al. 2016), otherwise the shape would be unrealistically elongated along its rotation axis. This resulted in shape models for 8820 asteroids. We further rejected all asteroids with two pole solutions that had differences in pole latitudes larger than 50 deg and with differences in longitudes smaller than 120 deg (Ďurech et al. 2016). These limits were set arbitrarily to filter out suspicious solutions with pairs of poles that were too far from the expectation that the ambiguity in the spin axis direction leads to two poles with the same latitudes and longitudes difference of 180 deg (Kaasalainen & Lamberg 2006). The number of asteroids then reduced to 8230.
As inĎurech & Hanuš (2018), we also used a model of a geometrically scattering triaxial ellipsoid to construct the periodograms The step in the period was set to 0.5 ∆P in this case. The advantage of this simpler shape model approximation is that the number of model parameters is smaller than in the case of convex models (only two parameters describing the shape: semiaxes a and b with the shortest axis normalized to unity), the shape model always rotates in a physically correct way (the ellipsoid rotates along its shortest axis), and false half-period solutions are absent. When assuming that the scattering is geometric, the disk-integrated brightness of a triaxial ellipsoid under a general viewing and illumination geometry is proportional to the projected area of the visible and illuminated surface, which can be computed analytically (Ostro & Connelly 1984). It makes this approach about a hundred times faster than when using convex shapes. However, ellipsoidal shape models are insufficient when an asteroid's light curve significantly differs from a simple double sinusoidal curve. The analysis of ellipsoid-based periodograms resulted in 16,010 unique periods. After finding the unique period, we modeled the data with the standard convex shape approximation and applied the same selection criteria described above. We derived solutions with one or two poles for 9910 asteroids, 7873 of which had physically plausible inertia tensors, and 7122 passed the check for the longitude and latitude difference between two pole directions.

Selecting stable period solutions
To detect and reject solutions that are not stable with respect to small perturbations of input photometric data, we performed a similar test to that performed inĎurech & Hanuš (2018). We randomly divided each dataset into ten parts, each part containing about one-tenth of the data points. We then removed one part, so about 10% of data, and repeated the period search and determined the period with the lowest χ 2 . We repeated this ten times, removing 10% of the points in each run and obtaining ten bestfit periods. We then compared these "jackknife" periods with the original best-fit period and selected only asteroids for which all ten jackknife periods were the same (within 0.1%) as the original period or for which there was at most one disagreement between the periods (i.e. there were nine jackknife periods that were the same as the original). After this selection, we obtained a set of 6269 models based on a convex-shape period search and 5784 models based on a period search with ellipsoids. There was some overlap between these two sets. For 3431 asteroids, we had a period solution from convex shape models and ellipsoids. Among them we found 20 asteroids for which the periods were different; we excluded them from further analysis, so the final set contained 8602 reliable models.

Comparison with Lightcurve Database
For 3690 asteroids in our sample, there was an independent period estimate in the Lightcurve Database (LCDB, version from 14 Dec 2021, Warner et al. 2009) that could be compared with our results. To consider only reliable LCDB periods, we selected asteroids with U = 3 uncertainty codes (848 cases). We identified 20 cases where the periods differed by more than 10%. They are listed in Table 1, where we also comment on each case.
Based on our inspection of original publications from which LCDB periods were compiled, we concluded that only five periods were incorrectly determined from Gaia DR3 data. With 848 periods in total, this represents a false solution rate of just ∼0.6%. We removed the five incorrect solutions from our dataset.

Comparison with Gaia DR2 results
The DR2 data are now part of the DR3; however, the data processing between DR2 and DR3 differs. Therefore, the reported magnitudes can be slightly different. Moreover, DR3 data usually contain additional epochs for each asteroid compared to DR2. Out of 129 models reconstructed from DR2 and published byĎurech & Hanuš (2018), 108 are among our final solutions from DR3. In four cases the rotation periods do not agree within their errors: asteroids (1540) Kevola, (2760) Kacha, (14410) 1991 RR1, and (21904) 1999 VV12. Apart from (2760) Kacha, the three asteroids were marked as less reliable inĎurech & Hanuš (2018) because when 10% points were removed, the period was not unique. For Kacha, the period of 53 h derived by an ellipsoidal model from DR2 is two times longer than our new value 26.5 h derived with a convex model.

Comparison with DAMIT
There are 1324 asteroids for which we have a model from DR3, and an independent model exists in the Database of Asteroid Models from Inversion Techniques (DAMIT, version from 7 Jun 2022,Ďurech et al. 2010). For 33 asteroids, their DAMIT periods are different from those we derived from DR3. We list these cases in Table 2. After checking the original data from which DAMIT models were reconstructed, we realized that in most cases, our DR3 solutions are more reliable, and DAMIT solutions are likely incorrect. These incorrect DAMIT solutions are usually based on sparse data only, often on the Lowell Observatory dataset with poor photometric quality. Only two of our 33 DR3 solutions are clearly wrong; we do not include them in our final dataset.

Final models
Out of the 8602 asteroids that passed all the checks (Sect. 2.1), 6 were excluded based on the period comparison with LCDB and DAMIT (Tables 1 and 2), and so there are 8596 final models. Their spin solutions are available at the CDS; the first and the last part of the table are shown in Table 3 as an example. New models that are not in DAMIT will be uploaded there. The false DAMIT models identified in Table 2 will be updated. As expected, the success rate of deriving a unique spin solution and a corresponding shape model depends on the number of observations. In Fig. 1 the orange histogram shows the distribution of the number of observations N in our final sample of models. The plot also shows the success rate as the ratio of the number of models we derived to the total number of asteroids in the DR3 sample for each bin. It starts at a few percent for N = 21 and increases to about 40% for N around 60. For even larger N, the number of asteroids per bin is small, and therefore the ratio fluctuates significantly.
Article number, page 3 of 14 A&A proofs: manuscript no. 45889 Table 1. Asteroids whose rotation period derived from DR3 was different from that in the LCDB.

Uncertainty of shape models
Compared to ground-based surveys, the accuracy of Gaia asteroid photometry is much higher. This high photometric accuracy enables us to reconstruct asteroid rotation states reliably from fewer measurements than we would need with ground-based data. With Gaia, we need tens of data points, while with less accurate photometry we would need hundreds of them (Ďurech et al. 2005, 2020). However, with only tens of observations the reconstructed shape model is sensitive to the number of observations, their distribution in time, and the model parameterization. We demonstrate this in Fig. 3, where we show shape models of asteroid (43) Ariadne reconstructed from a different number of observations going from 20 to 53 (the full DR3 set). The shape is also sensitive to the model resolution described by the order and degree of spherical harmonics series that describe the surface curvature . While spins derived from DR3 are reliable and stable with respect to the perturbations of input data (e.g., the dispersion in spin directions for different shape models in Fig. 3 is ±2 • in ecliptic latitude β and ±3 • in ecliptic longitude λ), the shape models are not, and any results based on shape parameters should take this into account. We do not report uncertainties of spin parameters in Table 3. We expect them to be smaller than 20 • , which was a mean difference between pole directions of models derived from Gaia DR2 and independent DAMIT models (Ďurech & Hanuš 2018). Uncertainty in the sidereal period P is typically a fraction (∼ 1/10) of the distance between local minima ∆P.

Results
In this section we analyze the spin and shape properties of asteroids we reconstructed from Gaia photometry (Sect. 2). We also search for correlations with other physical properties adopted from the literature. This includes diameters and family classification from the MP3C 4 database, and family ages from Nesvorný et al. (2015).
Due to the symmetry of the inverse problem (Kaasalainen & Lamberg 2006), we usually have two pole solutions with similar pole latitudes and differences in longitudes of ∼180 • . Their order in Table 3 is given by the quality of the fit, so always selecting the first one should not bias the analysis of physical properties. In some cases we plot both solutions in the figures. The poles are expressed in the ecliptic coordinate frame (λ, β). Computing the pole obliquity ε for each model (the angle between the spin vector and the normal to the orbital plane) requires a standard transformation into the orbital reference frame.
Due to the low number of optical measurements, the shape is only loosely constrained (Sect. 2.6). So instead of using the whole shape information, we utilize only the elongation a/b, which is determined as the axis ratio of the dynamically equivalent ellipsoid. It is the primary and most stable parameter of the shape.

Spin distribution in the main belt and beyond
The distribution of pole obliquities in the main belt is shown in Fig. 4. For visualization purposes, the scale of the vertical axis is linear in cos ε, so an isotropic distribution of spins would have a uniform distribution in cos ε. We also use this scale in other figures when plotting the distribution of obliquity ε. As expected (see, e.g., Hanuš et al. 2011), the distribution is far from uniform. Poles are clustered toward extreme obliquities with fewer asteroids with poles close to their orbital plane (ε ∼ 90 • ). Although there is a selection bias against asteroids with poles close to the Notes. The table lists for each asteroid the period P Gaia we derived from Gaia DR3 data, the period P DAMIT from the DAMIT database, the number of points N in DR3, the method used for computing periodograms: C -convex shape models, E -ellipsoids, CE -both methods provided the same period, and the reference to the DAMIT model. The last column gives our conclusion about the discrepancy between the models. ecliptic plane, it cannot be responsible for the observed distribution. Moreover, the distribution also depends on the asteroids' size. This has been interpreted as the evolution caused by the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect. The plot also nicely demonstrates the evolution due to the Yarkovsky effect, which causes the retrograde asteroids (ε > 90 • ) to migrate inward to the smaller semimajor axis, while prograde rotators (ε < 90 • ) migrate outward. This causes regions close to mean-motion resonances with Jupiter to be depleted; prograde rotators lean to the resonance from the left, get scattered when crossing it, and the space on the right does not contain prograde rotators. The situation is symmetric for retrograde rotators. We illustrate this coupling in more detail near the ν 6 , 3:1 and 7:3, and the outer edge of the main belt in Fig. 5. The concentration of asteroids near a p ∼ 3.174 au is caused by the 5-2-2 three-body resonance causing the chaotic diffusion in the Veritas collisional family (Tsiganis et al. 2007). This group has a smaller fraction of extreme obliquities compared to the main belt. As the Veritas family is quite young (∼ 8 Myr, Carruba et al. 2017), the spin vectors of its members are likely less evolved by the YORP.
Another fact that is demonstrated on the plots is that the prograde-retrograde distribution is not symmetric; spins of retrograde models are more tightly clustered toward the perpendicular orientation than those of prograde models. This might be caused by secular spin-orbital resonances, which primarily af-fect only the prograde rotators. However, we do not provide any further analysis here.
The asymmetry in positions of prograde and retrograde rotators is present, perhaps to a lesser extent, also in the Cybeles, Hildas, and Jupiter Trojans groups. However, we note two possible caveats here. These groups, on average, contain larger bodies than main-belt asteroids, due to their larger heliocentric distance and dark albedos (i.e., comparable S/N of the DR3 fluxes are for different sizes). In addition, due to the short time span of only about three years of DR3 data, there might be some observational bias correlated with the semimajor axis because the viewing geometry has changed less than for main-belt asteroids. For this reason, spin and shape determination for distant asteroids could be less reliable.

Distribution of rotation periods
In Fig. 6 we show the distribution of rotation periods with respect to the size and also include the reliable solutions from the LCDB (flag U = 3 or 2). Clearly, the two samples (i.e., Gaia and LCDB) are different due to the observational bias: Gaia lacks larger bodies as those were often saturated, but also small kilometer and subkilometer-sized objects that were either faint for Gaia or did not have frequent close approaches. The small asteroids in LCDB are almost exclusively NEAs that have rare A&A proofs: manuscript no. 45889 Notes. The table lists the first ten and the last ten asteroid models derived from Gaia DR3 data. The full table is available in electronic format at CDS. For each asteroid we report its spin axis direction in ecliptic latitude λ 1 and longitude β 1 (the second pole solution has ecliptic coordinates λ 2 and longitude β 2 ), the sidereal rotation period P, the number N of photometric measurements in DR3, and the method used for computing periodograms: C -convex shape models, E -ellipsoids, CE -both methods provided the same period. The uncertainty of the rotation period P is on the order of the last decimal place.
close approaches, but if they do, they can get bright enough even for small aperture telescopes. Moreover, it is relatively simple to obtain a reliable period during the close approach, but to derive a spin state and shape we need more approaches with different observing geometries. We do not have such data from DR3. DR3 data also allowed us to derive many solutions with rotation periods > 50 h. Not many such solutions are available in the LCDB, due to the observational bias, because long periods are challenging for ground-based observatories that produce dense light curves (Marciniak et al. 2015).
The most intriguing feature in the distribution of rotation periods is the group of slower rotators (P 50 h) separated by a gap from the faster rotators. A similar behavior was observed for the Jovian Trojans (Kalup et al. 2021). In the absence of the YORP effect, the period bimodal distribution was thought to be caused by the presence of slowly rotating bodies that are believed to be originally synchronized equal-sized binaries. This scenario is unlikely to explain the period bimodality in the population of the small main-belt asteroids. Instead, rotation periods of smaller asteroids (D < 30 km) are evolved by YORP; they gradually increase or decrease, sometimes even changing the sign of the evolution due to possible changes in surface topography (i.e., after a collision, due to close approaches with terrestrial planets, or after landslides due to fast rotation). If the period decreases, the asteroid can even start to rotate in the non-principle axis state (known as tumbling) as its rotation energy gets low enough to become excited even by a minor impact. However, this scenario predicts that the YORP-induced change in the period is smooth; there should not be any gap between the faster and slower rotators. In the group of slower rotators, their rotation period correlates with their size: larger bodies have longer rotation periods.
In Fig. 7 we illustrate the dependence of the pole obliquity on the rotation period. Bodies with periods P < 50 h are qualitatively similar to the general trends in the population; prograde and retrograde rotators both tend to have obliquities close to the YORP end states of ∼ 0 or ∼ 180 degrees, while prograde rotators exhibit a larger scatter. However, the slower rotators (P < 50 h) violate this behavior. The range of obliquities is larger for retrograde rotators and smaller for prograde ones. The behavior for prograde rotators with P > 100 h is particularly puzzling as basically all solutions have 0 < ε < 30 • .
We note a possible bias in our dataset concerning the slower rotators. Many tumblers were found in this population, and we should have many in our dataset as well. However, the convex inversion model we used assumes a relaxed rotation state and is inadequate for tumblers. This could be indicated by worse fits to the data; the average RMS of our solutions should be larger for the tumblers compared to the faster-rotating bodies. We observe larger RMS values for bodies with P > 50 h (see Fig. A.1), so it is likely that some asteroids for which we derived a model are not exactly in the state of principal-axis rotation. However, they are only slightly excited so that their sparse photometry can still be modeled with a single-period model.
The DR3 shape models exhibit a dependency between the size D and the shape elongation a/b (see Fig. A.2), which is in good agreement with previous reports (e.g., Cibulková et al. 2016;Szabó et al. 2022). Asteroids larger than ∼ 30 km are less elongated (i.e., more spherical) than the smaller ones.
Finally, Fig. A.2 illustrates the shape elongation dependence on the asteroid's rotation period. Highly elongated asteroids with a/b 1.5 are rare for rotation periods shorter than about 3 h. In general, the elongation increases with increasing rotation period with a peak near 6 h. The elongation then slowly decreases for periods in an interval of ∼ 6-50 h. After that, the elongation increases and has the largest mean values for this population of rotators. Our results agree with Szabó et al. (2022), who based their study on almost 10,000 rotation periods derived from TESS data. Contrary to our results, Szabó et al. (2022) did not report the increase in the mean elongation for asteroids with P > 50 h that we see in the DR3 models. However, it is unclear whether this is a real effect or a systematic bias related to the higher RMS values discussed in the previous paragraph.

Spin vector distribution in asteroid families
The Yarkovsky-driven evolution is also seen in asteroid families, where the 1/D dependence of the Yarkovsky drift on the diameter D causes the typical V-shaped spreading in the proper semimajor axis. In Figs. A.3, A.4, and A.5, we show 14 asteroid families with the highest number of asteroid models. In general, the distribution of prograde-retrograde models agrees with the theoretical expectations that prograde rotators are to the right (larger a) from the center of the family, while retrograde ones are on the left (smaller a) (see the plots for the Eos, Eunomia, Themis, and Koronis families). However, this trend is not as evident in some other families (e.g., Dora and Euphrosyne). Fami-lies that are truncated by resonances can consist of just one wing to which correspond either prograde or retrograde rotators. Typical examples of truncated families are Phocaea and Maria, both having "well-behaved" spin directions.
There is no family where all bodies follow the ideal Yarkovsky-driven evolution. Specifically, bodies close to the family center or at its outskirts often have spin that is inconsistent with the expected Yarkovsky drift. The former is likely due to the slow evolution of some bodies, especially if their initial locations in semimajor axis a were the most extreme, on the opposite side with respect to their sense of rotation. Bodies at the borders of the family are then likely interlopers mistakenly associated with the family by the HCM method (Zappalà et al. 1990). In general, interlopers can be present at any place in the family, but their number should not be larger than ∼10%. Moreover, non-catastrophic collisions can randomize the spin states of the family members at any location. The importance of these collisions is given by the collisional timescale and the family age. These timescales are dependent on the asteroid size and are usually on the same order as the age of the family for bodies with sizes of ∼10-30 km.
So far, we do not see any clear sign of the stochastic YORP evolution (Statler 2009;Bottke et al. 2015) that should cause the spin orientation of small family members to be oriented randomly, nor following the prograde-retrograde dichotomy of larger bodies.

Rotation periods in families
The distribution of rotation periods in several asteroid families was recently studied by Szabó et al. (2022). The authors report that rotation period distributions differ between the cores and outskirts of some collisional families (e.g., Flora and Maria) while are consistent among some other families. Moreover, the authors also show that the lightcurve amplitude distributions in families could be correlated with the family age, probably due to the temporal evolution of asteroid shapes. Here we focus on the same correlations based on the properties derived from DR3 photometry. Figure A.6 shows the distribution of rotation periods and elongations in various families. We ordered the families according to their approximate age (parameter c 0 from Nesvorný et al. 2015). For the distribution of the rotation periods, we selected only asteroids with P < 24 h to remove the population of slower rotators. We see a clear trend in Fig. A.6: older families have a larger fraction of asteroids with P = 10-24 h and with rounder shapes (smaller a/b), which likely indicates a more evolved population in the period due to YORP spin down.

Conclusions
Although the DR2 already showed the scientific potential of Gaia asteroid photometry (Ďurech & Hanuš 2018;Mommert et al. 2018;Colazo et al. 2021;Wilawer et al. 2022), it is the DR3 that enables us to make a significant leap forward in asteroid modeling. We derived ∼8600 models from Gaia DR3, which is more than a factor of two higher than what is currently available in the DAMIT database (∼3500). Considering the overlap between our models from DR3 and those already in DAMIT (about 1300 asteroids), we obtained ∼7300 new shape and spin state models. So we now have information about spin axis direction, rotation period, and shape for more than ten thousand asteroids.  The analysis of the distribution of asteroid spins we present in this paper confirms previous findings and the expected trends. Specifically, the spins are affected by the YORP effect: they evolve toward extreme obliquity values. Some asymmetry between prograde and retrograde rotators might be related to spinorbital resonances. Prograde and retrograde rotators have opposite Yarkovsky drifts on the semimajor axis, which leads to the separation of these two groups near the prominent mean-motion resonances and in asteroid families. We also see correlations between obliquity, rotation period, and shape elongation.
The results presented in this paper are based solely on Gaia DR3 data. The number of new asteroid models is large compared to the number of models currently known, but still small compared to the number of asteroids for which DR3 photometry is available. The next step is to combine DR3 with photometry from other surveys. This will increase the number of data points for individual asteroids, enlarge the time span of observations, and eventually lead to thousands of additional asteroid models.
In future Gaia data releases, the number of asteroids will not increase dramatically. However, the number of observations per object increases as Gaia continues to collect data, so with a three times longer observing window, for example, most asteroids will have more than 30 detections (see Fig. 1), and we expect that it will be possible to reconstruct spin states for tens of thousands of asteroids. With more data points, the uncertainty of spin parameters and shape will decrease.