© The Authors 2026

1 Introduction

The modern asteroid belt is primarily composed of collisional fragments of planetesimals. Groups of fragments that share a common parent body are known as asteroid families (Hirayama 1918). The long-term dynamical evolution of asteroid families in the Solar System is significantly influenced by non-gravitational forces, such as the Yarkovsky effect, which arises from the anisotropic emission of thermal photons from the surface of a rotating body and results in an acceleration that secularly modifies the body’s semi-major axis.

The Yarkovsky-induced orbital drift of asteroids in the main belt is ~ 10⁻⁵−10⁻⁴ au Myr⁻¹ for 1 km objects and scales as D⁻¹, with D being the asteroid diameter (Bottke Jr et al. 2006; Spoto et al. 2015), and it leads to the orbital dispersion of the asteroid family. The predictable nature of the Yarkovsky effect has made it an invaluable tool for estimating the ages of asteroid families. Over time, this size-dependent dispersion creates a characteristic V-shape in a plot of a family’s semi-major axis versus the inverse diameter of its members. The analysis of this V-shape, combined with an independent evaluation of the Yarkovsky drift rate, enables estimation of the time elapsed since the formation of the family (e.g., Vokrouhlickỳ et al. 2006; Walsh et al. 2013; Milani et al. 2014; Delbo’ et al. 2017).

However, the evolution is complicated by the dependence of the Yarkovsky effect on several parameters that are difficult to determine, including thermal inertia and bulk density. Assuming a similar composition and surface properties across a given family, the remaining dependencies are spin rate and obliquity. The rotation of asteroids is mainly affected by collisions and a thermal torque known as YORP (Yarkovsky–O’Keefe–Radzievskii– Paddack) that is induced by the scattering and reemission of solar radiation from an asteroid’s irregular surface (Paddack 1969; Rubincam 2000). The YORP effect applies a torque that, over long timescales, can significantly modify its spin vector, altering both the spin period and the spin-axis obliquity. Collisions dominate large and slowly rotating asteroids, whereas YORP is more effective for small and fast rotators (Pravec et al. 2008).

In recent years, data from the Gaia mission (Gaia Collaboration 2023) have led to the computation ofa vast number of new asteroid spin vectors (Durech & Hanus 2023; Cellino et al. 2024). This wealth of data has revealed new and unexpected trends. For example, Durech & Hanus (2023); Cellino et al. (2024); Vavilov & Carry (2025) analyzed the spin period distribution and confirmed an excess of slow rotators that are separated by a “valley” from the population of fast rotators in the period-diameter diagram. This excess of slow rotators is explained by Zhou et al. (2025) as being the result of the gathering of the slowly tumbling asteroids (i.e., asteroids in non-principal axis rotation) that suffer from a weaker YORP torque.

The aim of this work is to investigate the rotational evolution of asteroids among asteroid families from the perspective of observation. We introduce a new parameterization built on general features of the Yarkovsky drift that are common to different families. As described in Sect. 2, by exploiting the estimated family ages, we introduce a normalized timescale that allows us to build a common framework of interpretation. In Sect. 3, we discuss our results, which show that the current dataset of rotational properties for the members of 28 families provides useful new constraints on the evolution timescales of the asteroid rotations. In Sect. 4 we summarize our findings.

2 Methods

We selected 28 asteroid families with ages ranging from 14.34 Myr (Nele) to 3024.13 Myr (Themis) (Spoto et al. 2015). The families and their properties are listed in Table A.1. The membership of each asteroid to a family and the proper elements were retrieved from the Asteroids Dynamic Site (Milani et al. 2014). We are aware that the estimation of family ages could be affected by methodological limitations and uncertainties in the Yarkovsky calibration. These problems are addressed in Sect. 2.3, which also contains an assessment of the impact of errors, when estimating the diameter, on the final result.

The diameter, D, and albedo, p_V, were extracted from the WISE/NEOWISE survey (Masiero et al. 2011). For objects observed multiple times by the survey, we took the weighted average of their measurements using as weights the errors on D and p_V. Since the absolute magnitude, H, does not have an error associated with the measure, we just calculated the mean.

To increase the statistics for objects that do not have a measurement for the diameter, we calculated D using $$D=\frac{1329}{\sqrt{{\widehat{p}}_V}}\cdot {10}^{-0.2H}[\text{km}],$$(1) where $\widehat{p}v$ is the weighted average albedo of the family. Before the calculation of $\widehat{p}v$, possible outliers were discarded using the interquantile range. The absolute magnitude, H, was retrieved from the Asteroid Family Portal (Novaković et al. 2022). To calculate the error associated with $\widehat{p}v$, we performed a Monte Carlo bootstrap sampling of the errors associated with p_V and took the median of this set. We repeated the operation 5000 times. In the end, we calculated the standard deviation of the set of the medians calculated above.

Spin vectors were retrieved from Cellino et al. (2024) and Durech & Hanus (2023). In addition, we also included data from the Asteroid Lightcurve Database (LCDB). From this last source, we selected only observations with a quality flag of Q ≥ 3−. For objects with more than one suitable measurement, we selected the most recent one (Warner et al. 2009). In the end, we obtained 8739 objects with a spin period and 3794 objects with obliquity estimation. The complete dataset is reported in the Zenodo repository¹.

The lower statistics regarding the obliquity can be attributed to the difficulty in reconstructing light curves. The perioddiameter distribution of our sample is presented in Fig. 1, which further shows that this “valley” feature is present not only in the general asteroid population but also within individual families. It is important to note that the dataset is affected by observational biases. For example, it is easier to determine the spin state of fast rotators with respect to slow rotators because the latter require a longer observation time to cover a full rotation. Detecting super-fast rotators (P < 2.2 h) can also be challenging, particularly in sky surveys, due to an insufficient observation cadence (Novaković & Gutiérrez 2025). However, in our dataset only two asteroids, namely 1999 RF226 and 2008 VN1, can be classified as a super-fast rotator. In addition, the dataset is biased toward brighter asteroids, which are easier to observe.

We defined a dimensionless time, t, for each asteroid: $$t=\frac{{\tau }_{\text{age}}}{{\tau }_{\text{YORP}}},$$(2) where τ_age is the age of the family. Spoto et al. (2015) estimated the ages for both sides of the V-shape because families can be non-symmetrical in the (a, 1/D) space (e.g., in young families due to differences in the ejection velocities or in old families due to different impacts). For simplicity, we adopted the weighted average of the two age estimations, using as weights the errors associated with the measurements. The differences were taken into account in the error bars and are be addressed in Sect. 2.3.

The normal YORP timescale, τ_YORP, can be expressed as follows (Rubincam 2000; Vokrouhlickỳ & Čapek 2002; Golubov & Scheeres 2019): $${\tau }_{\text{YORP}}\sim 1{\left(\frac{D}{1\text{km}}\right)}^2{\left(\frac{a_p}{2.5\text{AU}}\right)}^2[\text{Myr}].$$(3)

Here, D is the diameter of the asteroid, and a_p is the proper semi-major axis. The YORP timescale is defined as the time it would take for an asteroid to double its rotation rate under the YORP torque (Rubincam 2000). The distribution of calculated r is shown in Fig. 2

Fig. 1 Observational data from Gaia (in gray) showing the perioddiameter distribution for asteroids. As purple squares we show the distribution of the Hygeia family, and as green triangles we show the distribution of the Massalia family. The dashed red line is the fit line that identifies the gap (Zhou et al. 2025). Spin period values were calculated by Durech & Hanus (2023) and Cellino et al. (2024). The ages in the legend are the weighted average ages computed from Spoto et al. (2015).

Fig. 2 Distribution of calculated t values for 8739 asteroids. The minimum t is tmin ≃ 0.003, while the maximum is tmax ≃ 1182.8.

2.1 Fraction of slow rotators

The excess of slow rotators was discovered in the past century (e.g., Harris & Burns 1979; Pravec & Harris 2000; Durech & Hanus 2023). Zhou et al. (2025) shows that they could be in a tumbling state caused by YORP de-spin. As a result, this excess of slow rotators is expected to emerge on ~ 2 YORP timescales. To examine this condition, we investigated the fraction of slow rotators and its evolution among asteroid families, as their ages are known, and we assumed that the initial spin state was set immediately after the family forming event.

We binned the data with respect to t. For each time bin, t_i, we calculated the fraction of slow rotators (i.e., asteroids with a rotation period P ≥ 30 hours): $$f_{\text{slow},i}=\frac{N_i(P\ge 30\text{h})}{N_i}.$$(4)

The error associated with the fraction is the 95% confidence interval of a binomial distribution. The error model used to compute the uncertainties on the bins is described in Appendix B. Based on a visual inspection of the data, it was clear that there is an increasing trend up to a certain value of t – a “break point.” After this point, the fraction of slow rotators remains roughly constant. For the analysis, we fixed the number of bins at 20: ten bins before and ten bins after the break point. This is a good compromise between the number of bins and the number of objects in each bin. The detailed fitting procedure is explained in Appendix B. The result for the best break point is $t_{\text{bp}}^{*}={20.225}_{-1.040}^{+1.063}$. The final fit of the data was done using an exponential-linear function described as $$f(t)=-\left(a{e}^{-jt}+mt+k\right).$$(5)

The parameters for the best-fit result are presented in Table B.1. Figure 3 shows the fit function with the original asteroid data and synthetic data from Zhou et al. (2025) for comparison. The general trend of f_slow increases and asymptotically converges to ~0.25, a result that we discuss in more detail in Sect. 3.

We note that the distribution of data points in t is not uniform. In particular, Fig. 2 shows an exponential decay of objects in the interval t ~ 200–450, which results in a visible gap in the binned trends shown in Fig. 3. This feature arises in the first place from the fact that we used a quantile function to bin the data (see Appendix B) and in the second place from the combined effect of the discrete age distribution of asteroid families, the strong dependence of t on the asteroid diameter (t ∝ D⁻²), and the observational incompleteness for small objects (≲ 6 km), especially in older families.

Fig. 3 Fraction of slow rotators as a function of the dimensionless time t. The green points are the binned-collected observational data in this work. Black triangles are simulation data from the Zhou et al. (2025) model, where the YORP timescale is manually increased by ten times compared to the normal YORP timescale. The red line is the fit function in Eq. (5). The shaded area is the 95% confidence interval of the fit. We do not show individual families because data are not sufficient for a statistical analysis.

2.2 Polarization fraction

Another interesting feature to study is the distribution of spin poles inside a V-shape. It is known that the Yarkovsky effect also depends on the obliquity of the spin vector (Bottke Jr et al. 2006). The Yarkovsky effect is responsible for the V-shape distribution of an asteroid family in the (a − 1/D) space (Vokrouhlický et al. 2015). Prograde rotators drift outward in semi-major axis, while retrograde rotators drift inward. We therefore classified an asteroid as aligned when its spin sign, relative to the orbital plane, matches the position according to the expected drift direction. We also included spins with obliquities ε ~ 90°. We are aware that this may introduce some ambiguity due to the errors on the spin axis direction (several degrees) that may result in similar probabilities for both spin directions. However, this concerns a small number of cases, and they do not change the global statistics – the asteroids with 80^° < ε < 100^° correspond to 1.82% of the dataset.

The polarization fraction in each bin is then $$f_{\text{pol},i}=\frac{N_i(\text{aligned\hspace{0.17em}objects})}{N_i}.$$(6)

We again binned the data into ten bins before and ten bins after the break point.

The result for the best break point is $t_{\text{bp}}^{*}={18.958}_{-1.095}^{+1.228}$. The final values for the parameters were found fitting the following function to the data: $$f(t)=-a{e}^{-jt}+b{e}^{-gt}+\mathrm{0.5.}$$(7)

The second half of the equation is an exponential decay function that approaches the value of 0.5, which is the expected value of the polarization fraction for a random distribution of spin vectors. The final result for the parameters is presented in Table B.1. The asteroid data and our fit are shown in Fig. 4, which shows that the function reaches the maximum value of ≃80% at t ≃ 15.07.

A similar gap in the t ~ 150–400 range is visible in the polarization fraction (Fig. 4) and has the same origin discussed in Sect. 2.1.

Fig. 4 Polarization fraction as a function of the dimensionless time t. The green points are the binned data; the red line is the fit function in Eq. (7). The black triangles are the distribution of the polarization fraction for the Eos family. The general trend is visible in single asteroid families. The shaded area is the 95% confidence interval of the fit.

2.3 Assessing uncertainties

Each asteroid in the dataset is associated with uncertainties in its diameter, D_i,err, and in the estimated age of its parent family, τ_age,err. Accounting for these sources of error may lead to variations in the inferred value of the optimal break point, $t_{\text{bp}}^{*}$. In addition, the uncertainty on the family age reflects possible methodological limitations in the age determination and in the calibration of the Yarkovsky drift.

To assess the sensitivity of our results to these uncertainties, we performed a Monte Carlo resampling of the dataset. For each asteroid, we independently resampled the diameter and family age according to D_i,ew ~ 𝒩(D_i, D_i,err) and τ_age,new ~ 𝒩(τ_age, τ_age,err), and we recomputed the corresponding dimensionless time, t. The resampled dataset was then analyzed using the same fitting procedure described in the previous sections and detailed in Appendix B, yielding anew estimate of the best break point. This procedure was repeated 500 times to build a distribution of break point values, from which we computed the median and the 95% confidence interval.

The resulting best estimations of the uncertainty-propagated break points are shown in Fig. 5, where they are compared with the original estimates. For the slow rotators fraction, we obtained $t_{\text{bp,err}}^{*}={19.883}_{-1.260}^{+1.116}$, while for the polarization fraction, we found $t_{\text{bp,err}}^{*}={18.926}_{-1.356}^{+1.129}$.

To quantify the agreement between the original and uncertainty-propagated break point estimates, we employed a two-tailed probability-to-exceed (P.T.E.) test (Raveri & Hu 2019). We first computed the absolute difference between the median values, $\left|t_{\text{bp}}^{\ast }-t_{\text{bp,err}}^{\ast }\right|$, referred to as the “observed difference.” We then resampled with replacement both the original and the uncertainty-propagated break point distributions, and we computed the difference between their medians for each realization. The P.T.E. was finally obtained as the fraction of resampled differences exceeding the observed difference in absolute value.

The inclusion of uncertainties in diameter and family age produced only minor variations in the inferred break point locations. The uncertainty-propagated estimates show substantial overlap with the original determinations, with largely overlapping 95% confidence intervals and median values consistent within the corresponding uncertainty bands (Fig. 5). The resulting P.T.E. values are ≃48.4% for the slow-rotator fraction and ≃71.6% for the polarization fraction, indicating statistically compatible results and supporting the robustness of the original break point estimates.

Fig. 5 Comparison between the best break point estimated from the original dataset, $t_{\text{bp}}^{*}$ (circle marker), and the best break point derived from the resampled dataset, $t_{\text{bp,err}}^{*}$ (square marker). The break points related to the fraction of slow rotators (polarization fraction) are shown in blue (red). The break points calculated from the resampled datasets are consistent with the original break points.

3 Implications

3.1 Excess of slow rotators

Our studies confirm the excess of slow rotators for the overall asteroid population. Significantly, for the first time, we have captured the process of its occurrence. We observed a rapid evolution in the fraction of slow rotators before the transition (t < 20.5), after which it converged to ~25%. This fraction is consistent with a previous study by Pravec et al. (2008). In their work, f_slow ~ 20% was derived from 268 asteroids with diameters D=3–15 km. Of them, 79 are members of either the Hungaria or Phocaea family, and most of them are older than 4 Gyr.

In early studies, this excess was empirically modeled by artificially reducing spin acceleration by a factor of two Pravec et al. (2008) or by stopping rotational evolution Bottke et al. (2015) beyond a prescribed period. However, these approaches are not grounded in a self-consistent physical mechanism.

In the model of Zhou et al. (2025), tumbling was introduced for the first time, and the YORP torque acting on tumblers was reduced by a factor of η_tumbling = 10 due to isotropic radiation under a chaotic tumbling rotation. The free parameter η_tumbling was calibrated to reproduce the slow-rotator fraction measured by Pravec et al. (2008). Consequently, the tumblers spin up or down more slowly than spinners (i.e., those in principle-rotation), leading to an accumulation of asteroids in slow tumbling states and naturally producing an excess of slow rotators. Therefore, f_slow is directly tied to the fraction of tumblers and increases with both η and the probability of a negative YORP torque. Figure 3 illustrates the evolution of f_slow in the model of Zhou et al. (2025) for η_tumbling = 50.

3.2 The stochastic YORP timescale

We find that at t ~ 20.23, f_slow starts to converge. In the models of Pravec et al. (2008) and Zhou et al. (2025), this corresponds to approximately two effective YORP timescales. The “static” YORP timescale, τ_YORP (Eq. (3)), refers to the time required for a body to evolve from its initial spin state to an end state under a fixed, unchanging YORP torque. In reality, however, YORP torques, and thus the evolution of the spin vector of an asteroid, often undergo a random walk (Bottke et al. 2015) because YORP is highly sensitive to small-scale surface features, such as craters and boulders, that may evolve (Statler 2009; Golubov & Krugly 2012; Cotto-Figueroa et al. 2015; Golubov & Lipatova 2022; Zhou et al. 2022; Zhou & Michel 2024).

To account for these stochastic processes, the effective YORP timescale is better represented by the “stochastic YORP timescale” rather than the “static” YORP timescale (Bottke et al. 2015), which can be estimated by $${\tau }_{\text{YORP,stoc}}={\eta }_{\text{stoc}}\cdot {\tau }_{\text{YORP}},$$(8) where η_stoc describes the random-walk evolution of the YORP torque over the body’s collisional history. This modifies not only the rotation evolution rate but also the Yarkovsky-driven spreading and V-shape evolution of asteroid families. However, η_stoc remains largely uncertain due to the complexity of the models.

This study provides the first direct observational validation of these theoretical predictions through our examination of the rotational evolution of asteroid families spanning a broad range of ages. The rotational model predicts that f_slow converges at ~2τ_YORP,stoc (Pravec et al. 2008; Zhou et al. 2025). Figure 3 shows both the observational asteroid data collected in this study and the simulation data from the model of Zhou et al. (2025), where the YORP timescale was manually increased by a factor of ten. The agreement between the simulation data and observations supports the stochastic YORP timescale being roughly an order of magnitude longer than the classical YORP timescale (i.e., η_stoc  10).

Given such a long stochastic YORP timescale, in combination with the interplay of collisions, YORP self-limitation could be triggered. Significantly, this phenomenon could explain why small asteroids evade transformation into top shapes and/or binaries. This is achieved by their rotational rates being constrained, thereby regulating the amount of angular momentum that can be imparted to a deformable body (Cotto-Figueroa et al. 2015). With YORP evolution proceeding about 20 times slower than classically predicted, many asteroids simply do not experience sufficient YORP cycles within their dynamical lifetimes to undergo catastrophic reshaping.

3.3 Polarization fraction

The evolution of the polarization fraction provides insight into the interplay between YORP-driven spin evolution and collisional reorientation. At early evolutionary stages (t ≲ 18.96), the polarization fraction increases, reaching a maximum of ≃80% at t ≃ 15.07. This rise indicates that Yarkovsky initially dominates the orbital evolution in the semi-major axis space, systematically aligning the asteroid spin axes with their orbital drift direction. In this initial stage, the YORP effect could increase the efficiency of the Yarkovsky drift by driving asteroids toward extreme obliquities (0° or 180°), thereby enhancing the correlation between spin orientation and semi-major axis position within the V-shape. However, the fact that polarization reaches only 80% rather than complete alignment suggests a combination of stochastic reorientation due to collisions and errors in the currently determined spin directions, the precise proportions of which remain unknown.

The subsequent decay of polarization to 0.5 at a larger t indicates complete randomization of spin axes relative to the orbital position. This transition suggests that collisional reorientation becomes increasingly important relative to Yarkovsky- and YORP-driven evolution at extended timescales. The cumulative effect of non-catastrophic collisions progressively breaks the connection between the current spin state and the Yarkovsky drift history that determined the asteroid’s position in the V-shape. This interpretation is consistent with theoretical models suggesting that for smaller objects (i.e., more evolved family members), collisions occur more frequently and with greater strength, causing bodies to transition between different spin states more rapidly (Marzari et al. 2011).

The randomization of spin axes through collisions has implications for the interpretation of asteroid family ages derived from the V-shape. The fundamental hypothesis underlying V-shape age dating is that the direction and drift rate of the semi-major axis for family members located at the edges of the V-shape remain constant over time (Milani et al. 2014). However, collisional spin reorientation introduces a random-walk component to the Yarkovsky drift. Under these conditions, the inverse slopes of the V-shape borders are no longer proportional to the family age. Instead, they scale with the square root of the family age (Marzari et al. 2011). The decay of polarization toward a random state at large dimensionless times indicates that collisional spin reorientation becomes increasingly important for family members older than t ≈ 20, marking the limit beyond which V-shape age estimates become unreliable. For older families, partial erasure of the Yarkovsky signature may bias V-shape results. In such cases, the degree of polarization could serve as an independent diagnostic or correction factor. Further work is needed to quantify these effects.

4 Conclusions

In this work, we have utilized the rapidly expanding dataset of asteroid spin states to empirically investigate the long-term rotational evolution of main-belt asteroid families. By introducing a dimensionless evolutionary timescale (t = τ_age/τ_YORP), we directly compared 28 asteroid families with a wide range of ages and provided new observational constraints on the coupled YORP-Yarkovsky collision processes.

We showed that the fraction of slow rotators increases rapidly with t and saturates at 0.25 near t ≈ 20. This behavior matches theoretical predictions for the onset of stochastic YORP selflimitation and yields an effective stochastic YORP timescale ten times longer than the classical formulation (η_stoc  10). Our results therefore provide the first direct population-level confirmation that YORP evolution is significantly slowed by surface-morphology-driven torque variability.

We further showed that the polarization fraction – quantifying the spin–orbit correlation within V-shapes – initially rises to ~80% as YORP drives bodies toward extreme obliquities and efficient Yarkovsky drift. However, beyond t ≈ 20, the polarization decays toward the random limit of ~50%, indicating that collisional reorientation increasingly dominates the spin evolution of older families. This transition marks the approximate age boundary beyond which V-shape estimates become systematically biased and suggests that the degree of polarization may be used as a diagnostic or correction factor for advanced families.

The rotational evolution trends of f_slow and f_pol can be further tested with the forthcoming Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST), which is expected to increase the available spin-state sample by at least an order of magnitude (Ivezić et al. 2019). With such datasets, these trends could provide an additional dimension for constraining asteroid family ages, for example, by fitting τ_age to the evolutionary tracks in Figs. 3 and 4 for individual families.

Acknowledgements

W.-H. Zhou acknowledges funding support from the Japan Society for the Promotion of Science (No. P25021). G. Bertinelli thanks G. Viterbo for helpful discussions on topics related to this work. The authors thank the anonymous referee for productive comments and suggestions.

References

Appendix A Asteroid families

Appendix B Fitting procedure

Since the distribution of t (Eq. 2) is not uniform (see Fig. 2), we used the quantile-based discretization function² to roughly have the same number of asteroids in each bin. The value of t at which each bin is centered is associated with an error, given by the standard deviation of the dimensionless times inside the bin. The fraction of slow rotators in each bin is calculated as the number of slow rotators divided by the total number of asteroids in the bin. The error associated with the fraction is the 95% confidence interval of a binomial distribution. In particular, we used the Wilson score interval with continuity correction (Newcombe 1998).

We fit the data with respect to the slow rotators fraction (Sect. 2.1), with a piecewise function defined as $$f(t)=\{\begin{array}{ll}-\left(a{e}^{-jt}+k\right)\hfill & \text{if}t\le t_{\text{bp}}\hfill \\ -(mt+b)\hfill & \text{if}t>t_{\text{bp}},\hfill \end{array}$$(B.1) where t_bp is the break point and a, j, k, m, and b are the parameters to be fit.

Parameter estimation We employed a weighted Huber loss function to ensure robustness against outliers: $$L(r)=\{\begin{array}{ll}\frac{1}{2}{r}^2\hfill & \text{if}\left|r\right|\le \delta \hfill \\ \delta \left(\left|r\right|-\frac{1}{2}\delta \right)\hfill & \text{otherwise}\hfill \end{array}\text{with}r=\frac{y-f(t)}{\sqrt{{\sigma }^2}},$$(B.2) where δ = 1.345 (Huber 1964). The weighted residual r accounts for error propagation through $\sqrt{{\sigma }^2}=\sqrt{y_{\text{err}}^2+{\left({\partial }_{t}f(t)\cdot t_{\text{err}}\right)}^2}$, where y_err and t_err are the errors associated with the fraction of slow rotators (polarization) and with t, respectively. Since the error on the fraction of slow rotators (polarization) is asymmetric (from the Wilson confidence interval), we adopted the following scheme: $$y_{\text{err}}\text{=}\{\begin{array}{ll}{y}_{\text{upper}}-y\hfill & \text{if}r\ge 0\hfill \\ y-y_{\text{lower}}\hfill & \text{otherwise}\hfill \end{array}$$ otherwise.

The total loss was calculated as ℒ_Huber = Σ_i L(r_i). The Huber loss provides robust parameter estimates that are less sensitive to outliers compared to ordinary least squares, and does not require errors to be normally distributed (Huber 1964).

Break point selection We implemented an iterative scan over candidate t_bp values to identify the optimal model complexity· At each candidate break point, we fit the model parameters using the Huber loss, then evaluated model quality using the Akaike information criterion (AIC): $$\text{AIC}=2k+n\ln (\text{RSS}/n),$$(B.3) where k is the effective number of free parameters, n is the number of data points, and RSS is the weighted residual sum of squares: $$\text{RSS}={\displaystyle \sum }_i{\left(\frac{y_i-f\left(t_i\right)}{\sqrt{{\sigma }_i^2}}\right)}^2.$$(B.4)

Note that while parameter estimation used the Huber loss, model selection employed RSS for AIC calculation to maintain consistency with standard information-theoretic model comparison (Burnham & Anderson 2004)· The AIC penalizes model complexity and prevents overfitting, with lower values indicating better models· The break point that minimized AIC was selected as optimal·

The parameters were fit using the Sequential Least Squares Programming (SLSQP) algorithm from the SciPy package in Python, with an equality constraint enforcing the continuity of f(t) at t_bp.

The scan was initialized with a wide range of t_bp values and then refined by manually reducing, in an iterative fashion, the interval around the values that provided the lowest AIC value· To have a statistical estimate of the parameters, we performed a Monte Carlo bootstrap sampling of the data, i.e., we randomly sampled the data with replacement and fit the piecewise function to each set· We repeated the process 500 times· The values of the final parameters are the median of the distribution obtained from the procedure, while the errors are given by the 95% confidence interval of the distributions.

The fitting procedure for the polarization fraction (Sect. 2.2) is the same as explained above, and the best break point $t_{\text{bp}}^{*}$ was found by fitting the following function: $$f(t)=\{\begin{array}{cc}-\left(a{e}^{-jt}+k\right)& \text{if}t\le t_{bp}\\ b{e}^{-gt}+0.5& \text{if}t>t_{bp}.\end{array}$$(B.5)

The second part of the piecewise function is an exponential decay function that approaches the value of 0.5, which is the expected value of the polarization fraction for a random distribution of spin vectors.

In the end, for both the slow rotators fraction and the polarization fraction, we fit the data with a unique function (Eq. 5 and Eq. 7, respectively), fixing the break point at $t_{\text{bp}}^{*}$. The final parameters were again obtained from a Monte Carlo bootstrap sampling after 1000 simulations.

All Tables

All Figures

Fig. 1 Observational data from Gaia (in gray) showing the perioddiameter distribution for asteroids. As purple squares we show the distribution of the Hygeia family, and as green triangles we show the distribution of the Massalia family. The dashed red line is the fit line that identifies the gap (Zhou et al. 2025). Spin period values were calculated by Durech & Hanus (2023) and Cellino et al. (2024). The ages in the legend are the weighted average ages computed from Spoto et al. (2015).

In the text

	Fig. 2 Distribution of calculated t values for 8739 asteroids. The minimum t is tmin ≃ 0.003, while the maximum is tmax ≃ 1182.8.
In the text

Fig. 3 Fraction of slow rotators as a function of the dimensionless time t. The green points are the binned-collected observational data in this work. Black triangles are simulation data from the Zhou et al. (2025) model, where the YORP timescale is manually increased by ten times compared to the normal YORP timescale. The red line is the fit function in Eq. (5). The shaded area is the 95% confidence interval of the fit. We do not show individual families because data are not sufficient for a statistical analysis.

In the text

	Fig. 4 Polarization fraction as a function of the dimensionless time t. The green points are the binned data; the red line is the fit function in Eq. (7). The black triangles are the distribution of the polarization fraction for the Eos family. The general trend is visible in single asteroid families. The shaded area is the 95% confidence interval of the fit.
In the text

Fig. 5 Comparison between the best break point estimated from the original dataset, $t_{\text{bp}}^{*}$ (circle marker), and the best break point derived from the resampled dataset, $t_{\text{bp,err}}^{*}$ (square marker). The break points related to the fraction of slow rotators (polarization fraction) are shown in blue (red). The break points calculated from the resampled datasets are consistent with the original break points.

In the text