Thermophysical Investigation of Asteroid Surfaces II: Factors Influencing Grain Size

Asteroid surfaces are subjected to mechanical weathering processes that result in the development and evolution of regolith. Two proposed mechanisms--impact bombardment and thermal fatigue--have been proposed as viable and dominant weathering processes. Previously, we compiled and estimated thermal inertias of several hundred asteroids (mostly in the main-belt) for which we determined dependencies on temperature, diameter, and rotation period. In this work, we estimate grain sizes of asteroid regoliths from this large thermal inertia dataset using thermal conductivity models. Following our previous work we perform multi-variate linear model fits to the grain size dataset and quantify its dependency on diameter and rotation period. We find that the best-fit model fit indicates that asteroid grain sizes are inversely dependent on object size for<10 km asteroids and exhibits no relationship above this size cutoff. Rotation period and grain size show a positive relationship when the rotation period is greater than ~5 hr, and an inverse relationship below this rotation period. We conclude that both impact weathering and thermal fatigue are significant regolith evolution mechanisms. Furthermore, we run post-hoc t-tests between spectral groups to identify compositional differences among our asteroid set. Notably, suspected metal-rich, M-type and E-type asteroids have larger than expected grain sizes, and P-types have distinctly smaller grains than other groups


Introduction
The study and characterization of asteroid regolith-the unconsolidated, heterogeneous, rocky material covering the surface of planetary bodies (Shoemaker et al., 1969)-is an important part of understanding the processes and evolution of airless bodies of the Solar System. Generally speaking, asteroid surfaces evolve from poorly-sorted, blocky mixtures to well-sorted fine-grained regolith (Hörz et al., 2020). Thermal inertia, Γ, is a thermophysical property that can be used to characterize asteroid regolith. Thermal inertias lower than that of bare rock (Γ ∼ 1500-2500 J m −2 K −1 s −1/2 Jakosky, 1986;Bandfield et al., 2011) indicate the presence of regolith, with lower values signaling a finergrained surface. Although some meteorites have measured thermal inertia lower than this (Γ ∼ 1000 J m −2 K −1 s −1/2 ; Opeil et al., 2010;Opeil et al., 2020), the range of many asteroid thermal inertias (< 150 J m −2 K −1 s −1/2 MacLennan and Emery, 2021) is significantly lower. The bulk thermal inertia can be expressed as Γ = k eff ρ grain c s (1 − φ) in which k eff is the effective thermal conductivity of the regolith, ρ grain is the density of a rock with no void spaces, c s is the specific heat capacity of the material, and the regolith porosity, φ, is explicitly accounted for. Some mechanisms have been suggested as the primary drivers of regolith mechanical weathering on airless bodies: meteoroid impacts and thermal cycling.
Small meteoroid (< 1 m) impacts can cause breakdown of surface material over time. Some of the energy from an impact is partitioned into fragmenting near-surface material of the target body, albeit the exact fraction of energy that goes into this mechanical work is uncertain (Horz and Cintala, 1997). Some of the energy is partitioned into launching particles off the surface on various trajectories (ejecta), which partly depend on the proximity to the impact site and on the mechanical properties of the impactor and target. Particles that do not reach escape velocity return to the surface as newly-formed regolith. Other factors, such as the target porosity and strength, also play relevant roles in the production of craters and of regolith. For example, Housen and Holsapple (2003) found a clear inverse relationship between the ejecta:impactor mass ratio and porosity of the target. Impacts that excavate into fractured bedrock will generate more ejecta blocks. Thus, regolith evolution driven by impacts is dependent on on the size of the body, strength of the material, among other factors.
Internal stresses caused by differential thermal expansion, as a result of cyclic heating and cooling, lead to regolith breakdown and erosion (Eppes et al., 2015). The efficiency of this thermal cycling process has been modeled under the thermal environments of airless solar system bodies (e.g., Molaro and Byrne, 2012;Molaro et al., 2017), and its feasibility has been experimentally demonstrated on meteorite samples (Delbo' et al., 2014;Libourel et al., 2020).  demonstrated that large cracks observed on rocks seen across Bennu's surface were consistent with thermal fatigue model predictions. In general, increasing the heterogeneity of a rock (by changing the mineral grain boundaries or pores of empty space) can raise the peak stresses by up to a factor of three, compared to a homogeneous rock . Thermal fracturing acts on spatial scales spanning many orders of magnitude, and modeling work has shown that the efficiency depends on heliocentric distance, rotation period, and material thermomechanical properties (Ravaji et al., 2019;El Mir et al., 2019). The propagation of the thermal wave within a boulder initiates one or more cracks at the micro-scale throughout the boulder interior and throughout the diurnal cycle (Hazeli et al., 2018). These micro-cracks most likely originate near mineral grain boundaries, preferentially grow perpendicular to the local surface, and ultimately intersect other fractures reaching sizes of several centimeters (Delbo' et al., 2014;Molaro et al., 2015). Growth of large-scale cracks occurs in the direction of the heat flow throughout the boulder, and are thus able to transverse rocks that are several times larger than the initial fracture(s). Although rapidly rotating asteroids exhibit small diurnal temperature differences, Delbo' et al. (2014) showed that even a short heating cycle of 2.2 hrs is sufficient to generate and grow cracks. Lastly, Molaro and Byrne (2012) claim that fast-rotators ought to facilitate larger maximum thermal gradients, thus predicting a rotation period dependence of the overall efficiency.
This work (part II) continues our thermophysical investigation of asteroid surface properties by estimating thermally-characteristic grain sizes from thermal inertias (part I; MacLennan and Emery, 2021). We use these grain sizes to study regolith evolution in the context of two proposed regolith evolution mechanisms: meteoroid impact degradation and thermal fatigue cycling. These processes have been proposed as relevant for the creation and subsequent evolution of regolith on asteroid surfaces. The relevance of each weathering mechanism can be examined by comparing the relative correlation between regolith grain size and asteroid diameter or rotation period. Specifically, we claim that if regolith development and evolution is highly dependent on meteoroid impacts then we expect a statistically-significant correlation between the grain sizes and asteroid diameters. Similarly, if thermal cycling is effective then we should observe that grain sizes are correlated with rotation periods. These two mechanisms are not mutually exclusive, thus we consider the possibility that both are relevant for asteroid surfaces by using a multi-variate linear model in our analysis (Sec. 3.2). We use this multi-variate model to test the following hypotheses: 1) larger asteroids exhibit evolved regoliths, quantified by smaller grain sizes, and 2) slowly-rotating asteroids have a poorly-developed regolith, characterized by larger grain sizes. We also investigate how the regolith properties of main-belt asteroids are related to spectral classification using post-hoc tests of the fitted model residuals (Sec. 3.3), and orbital properties of near-Earth asteroids (Sec. 3.4). Finally, we note that it is currently unclear in what ways these grain sizes characterize a heterogeneous mixture of e.g., regolith grains and boulders. Thus, although the grain sizes estimated herein may be inaccurate in the absolute sense (Sec. 2.4) we utilize the modeling results to interpret the trends (Sec. 4) within the asteroid population in the context of regolith evolution mechanisms.

Methods
A handful of models have been developed for estimating grain size from thermal inertia of planetary surfaces by experimentally observing and/or theoretically modeling heat flow in granular media. For example, laboratory experiments have been able to observe empirical effects of particle size on k eff (Watson, 1964;Presley and Christansen, 1997;Presley and Craddock, 2006;Sakatani et al., 2018). Since, for airless planetary regoliths, Γ is primarily, but not exclusively, influenced by changes in k eff , an estimate of grain size can be made from Γ by using an appropriate thermal conductivity model for a granular medium (e.g., Gundlach and Blum, 2013;Sakatani et al., 2017;Woods-Robinson et al., 2019;Wood, 2020). In general, for airless bodies k eff has a solid-state component that describes heat conduction through grains and across grain contacts, and a radiative component that describes heat radiated across pore spaces (e.g., Watson, 1964): The k 1 and k 2 coefficients in Eq. 1 are dependent on the material properties of the regolith, such as grain size, packing fraction, and amount of contact between the grains (Watson, 1964). Note that the solid component of thermal conductivity is temperature-dependent and varies with composition (Wood, 2020, and references therein). These coefficients are approximated in the works of Gundlach and Blum (2013) and Sakatani et al. (2017), from which measurements of Γ (combined with the compositional information about the asteroid) can be used to estimate a characteristic grain size for an asteroid. We thus employ both these models, with an anisothermality correction from Ryan et al. (2020) applied to the radiative conductivity term, to estimate characteristic regolith grain sizes for each asteroid. The Gundlach and Blum (2013) formulation of k 1 and k 2 are calibrated using laboratory heat-flow measurements of lunar regolith (Chan and Tien, 1973;Gundlach and Blum, 2012). On the other hand, Sakatani et al. (2017) used mixtures of powdered glass beads of varying sizes and porosities. The surface temperature is required input for these models, from which we use the color temperature reported in Part I (MacLennan and Emery, 2021). These models also need compositional information about the material, which we infer from the asteroid's spectral taxonomic type when available. For objects for which there is no spectral information available, we infer the spectral type from the geometric albedo. Many asteroids belong to a dynamical family, in which case we can infer the spectral type and perform a cross-check with the geometric albedo or reflectance spectrum, if available. We assign a meteorite analog to each spectral type as outlined in Sec. 2.2 and account for uncertainties in the thermal inertia and various material properties by using a Monte Carlo method (Sec. 2.3) when implementing the model.

Thermal Conductivity Modeling
The solid thermal conductivity component, k 1 , is modeled by Gundlach and Blum (2013) via computing the efficiency of heat transfer through the surface contacts within a network of grains: k grain is the thermal conductivity of a grains with zero porosity, ψ = 1−φ is the volume filling factor (packing fraction). The empirically-derived constants f 1 = (5.18 ± 3.45) × 10 2 and f 2 = 5.26 ± 0.94 encapsulate information about the path of contact chains in the direction of heat flow (Gundlach and Blum, 2012). Thermal conductivity measurements of Apollo 11 & 12 samples are used to estimate Ξ = 0.41 ± 0.02, which incorporates and accounts for the irregular shapes of the particles and heterogeneity of regolith on the whole. The factors contained within brackets in Eq. 2 approximate the adhesive forces between regolith grains, which dominate over gravity on small bodies, to estimate the contact area between them.

4
The contact radius between grains, r c , is calculated by JKR (Johnson-Kendall-Roberts) theory (Johnson et al., 1971), assuming that adhesive forces are dominant: where ν is Poisson's ratio, E is Young's Modulus, and γ(T ) = T · 6.67 × 10 −5 Jm −2 is the specific surface energy of each grain-a measure for the adhesive bonding strength between grains.
In Sakatani et al. (2017), the solid conductivity component is modeled in a similar way to the formula above: k Sak where r c is the contact radius, also represented by Eq. 3. The coordination number, is the average number of particles that are in contact with each other and is a function of the regolith porosity. Here, f = 0.07318 + 2.193φ − 3.357φ 2 + 3.914φ 3 . The factor ξ is dependent on the shape and smoothness of the particles and is equal to unity for perfectly smooth spheres. We use ξ = 0.4, which approximates rough, non-spherical particles and is well within the range of experimentally-derived values from Sakatani et al. (2017). The radiative thermal conductivity (k rad ) coefficient is calculated by Gundlach and Blum (2013) to be: with the bolometric emissivity, , Stefan-Boltzmann constant, σ 0 ≈ 5.6710 −8 W m −2 K −4 , and the empirically-derived coefficient e 1 = 1.34 ± 0.01 (Dullien, 1979;Gundlach and Blum, 2012). The mean free path of a photon between regolith grains is directly proportional to their size, hence the r g factor in Eq. 6. In Sakatani et al. (2017) the radiative heat transfer coefficient is calculated via: where ζ is an enhancement factor for which we use ζ = 1.7-an average of the experimentallyderived values in Sakatani et al. (2017). However, we note that k rad values presented in Sakatani et al. (2017) show a clear inverse dependency on r g , which may be related to the effect of short-range (immediately-adjacent) versus long-range (non-adjacent) radiation exchange that is dependent on the particle size (van Antwerpen et al., 2012). Thermal gradients that exist within individual regolith particles cause a non-linear relationship between r g and k rad . Such thermal gradients will arise in large particles, when the solid component of thermal conductivity is low, or for particles that have significant intra-granular porosity. In these cases, the radiative component of thermal conductivity will be less than the T 3 theoretical relationship. Building off the work of van Antwerpen et al. (2012), Ryan et al. (2020) use the dimensionless parameter Λ s in an updated version of the non-isothermal correction factor, f k : and f k = a 1 tan −1 (a 2 Λ −a 3 s ) + a 4 .
When Λ s >25 then f k is set to unity, otherwise f k can be calculated via Eq. 9 where a 1 =-0.568, a 2 =0.912, a 3 =0.765, and a 4 =1.035. To account for this non-isothermality, we multiply the coefficients calculated in Eq. 6 and Eq. 7 by f k .

Input Parameters
Some of the required input parameters to the equations listed above are dependent on the composition of the material. Since meteorites are samples of asteroids, we attempt to establish a meteorite analog that is most appropriate for each asteroid in this study. The most direct approach can be made when high quality spectra have been acquired for an object. Various taxonomic systems have been defined based photometric colors and albedos (Tholen, 1984), or absorption features with spectral slopes of reflectance spectra (Bus and Binzel, 2002;DeMeo et al., 2009). The taxa defined in these works relate to different compositions, although connections can be ambiguous in many cases; particularly, for featureless spectra. Visible spectra and color information can be used to broadly distinguish between the S-complex (any spectral taxon beginning with an "S" and Q-types that are the unweathered endmember), K-type, V-type, B-type, C-complex (any taxon beginning with a "C"), and Bus-DeMeo X-complex. For the X-complex we adopt the Tholen (1984) E, M, and P-type system, which are distinguished by their geometric albedos. We use V-band geometric albedo cutoff of p V 0.42 for E-type objects, 0.12 < p V < 0.42 for M-types, and p V 0.12 for P-types. Many asteroids in our sample have no color or spectral information available, thus we are left to infer their compositional information by other means. For this exercise we use the geometric albedo reported and compiled by MacLennan and Emery (2021), and V-band slope parameters (G V ) from Oszkiewicz et al. (2011), which have both been shown to correlate with spectral class (e.g. Oszkiewicz et al., 2012;Veres et al., 2015). We found that using cutoff criteria only for p V was sufficient, but still utilize G V values as a consistency check that each object is within the range of expected values for its spectral class. All objects' p V and G V values are included in Fig. 2. It is interesting to note the positive correlation between p V and G V for the entire dataset (across taxonomic groups), yet, within each taxonomic group this correlation is absent or exhibits an inverse relationship. Such a scenario is an example of Simpson's statistical paradox 1 (Simpson, 1951;Yule, 1903). We briefly note this paradox group and analyze data from many sources, when appropriate, to account for temperature and porosity effects. We note that, in our albedo-based classification, M-types and lowalbedo V-types could potentially be misclassified as S-types and, similarly, some P-types may be misclassified as C-types. However, we see no obvious reason to re-assign any of these S-types to the M-type group, as their thermal inertias are inconsistent with a high thermal conductivity which would otherwise suggest a metal-rich surface. On the other hand, some P-types might be misclassified as C-types. We change the classification of (4003) Schumann from C-type to P-type, based on its location in the outer part of the Main Belt (≈ 3.4 au) where P-types are more abundant over C-types Carry, 2013, 2014). The difference in assumed material properties between C and P-types that we use is not large and would not alter the reported grain size by more than a few percent.

Thermal Conductivity
Laboratory measurements of meteorite thermal conductivities reveal dependence on the material porosity and temperature. In their review paper of thermal conductivites, Flynn et al. (2018) show that laboratory measurements of thermal conductivity do not change significantly in the temperature range from 100 and 300 K, but approach zero at 0 K. Additionally, void spaces can significantly impede the net solid-state heat flow within a meteorite sample (Soini et al., 2020). To account for porosity effects, we estimate the thermal conductivity at zero porosity using data collected at the same temperature.
In order to estimate thermal conductivity at φ = 0, lines are fit to meteorite samples with porosities less than 12% from the literature (Table 1), as shown in Fig. 3. Since the carbonaceous chondrites with measured k eff had porosities exceeding this cutoff, we did not use any for our analysis. We found that enstatite chondrites have a larger thermal conductivity than the ordinary chondrite samples (and one Shergottite), likely due to mineralogical differences. The y-intercept of the best-fit lines give the thermal conductivity at zero porosity of 4.05 and 6.23 W m −1 K −1 for ordinary chondrites and enstatite chondrites, respectively. Because the thermal conductivity for terrestrial enstatite is 90% that of enstatite chondrites at 275 K we scale the function downward by 10% to approximate the aubrite thermal conductivity, which are mostly comprised of enstatite. We used the same line slope to extrapolate a single measurement of thermal conductivity (Opeil et al., 2010) of the FeNi meteorite Campo del Cielo to estimate the thermal conductivity at φ = 0 are used in the thermal conductivity model to compute grain size.
Data presented by Soini et al. (2020) shows that the thermal conductivities of ordinary and carbonaceous chondrites are in agreement when the sample porosities are taken into account. Because we account for sample porosities and are using grain densities, we use the same zero-porosity thermal conductivity value for carbonaceous chondrites and ordinary chondrites, with the exception of CM chondrites. Newer thermal conductivity measurements of CM chondrites from Opeil et al. (2020) are also consistent with other carbonaceous chondrites. However, the spacecraft observations of B-type Bennu and Cg-type Ryugu revealed surfaces with weak, low-conductivity boulders that have thermal inertias which are ∼2-3 times smaller than stronger, high-conductivity boulders (Grott et al., 2019;Rozitis et al., 2020). This stronger material may be identical to CM chondrites found on Earth, yet, the weak material is likely less abundant among meteorite collections (Popova et al., 2011). In order to account for the presence of this weak material that is likely present on Ch and B-type asteroids, we adopt a lower thermal conductivity for both. Because this weaker material is most likely primitive in origin, we assign the low thermal conductivity to CI chondrites as well.
Our computed thermal conductivity values agree with the Soini et al. (2020) model fits to essentially the same dataset as ours. We accounted for temperature-dependent thermal conductivity among these samples by bootstrapping our k grain calculations with the data compiled in Flynn et al. (2018). Ordinary chondrites, HEDs, and carbonaceous chondrites do not demonstrate a significant dependence on temperature (Opeil et al., 2012;Flynn et al., 2018) in the range relevant to this study, so we use constant values. For FeNi and enstatite chondrite curves, we scale the curves presented in Flynn et al. (Fig. 15; to our values of k grain at 200 K. FeNi meteorites show a roughly linear dependence, and enstatite chondrites (and by deduction, the enstatite-dominated aubrites) show an inverse dependence on temperature. These temperature-dependent k grain functions, which are valid for 200-350 K, are listed in Table 3 and are assigned a 10% uncertainty in the thermal conductivity model.

Grain Density and Heat Capacity
Meteorite specimens, while not regolith, still contain small pores and cracks that lower the bulk density of the sample (e.g., Cadenhead and Stetter, 1975). These micro-porosities, φ, range from 4% to 10%, but can be as high as 20% for weathered fines and as low as 2% for iron meteorites. We wish to use the grain density (ρ grain ), the density of a material with zero porosity, as input to the thermal conductivity models. We primarily use the findings Macke (2010), who measured the porosities of individual meteorites with various compositions and computed ρ grain . Table 3 lists grain densities for many meteorite analogs related to their respective spectral types. For the S-complex and V-types we use a uniform distribution for ρ grain that spans the values reported in Macke (2010), instead of a normal distribution.
The ability of material to store thermal energy per unit mass is quantified by its specific heat capacity, c s . Laboratory heat capacity measurements of several meteorites reveal a clear temperature-dependence (Beech et al., 2009;Opeil et al., 2012;Szurgot et al., 2012;Consolmagno et al., 2013;Wach et al., 2013;Macke et al., 2019;Opeil et al., 2020). We account for this dependence by performing an independent meta-analysis-similar to that on the thermal conductivity dataset above-in which we use literature data on meteorite heat capacities measured at various ambient temperatures. While Flynn et al. (2018) perform a similar analysis, they present fits for a temperature range of 75-200 K, which does not cover the full range of asteroid surface temperatures relevant to this work. Our empirical fits to the data are in agreement to the semi-empirical functions of Macke et al. (2019), who used a larger dataset. Our results are very similar to theirs over the same temperature range. Opeil et al. (2020) present temperature-dependent heat capacities for five CM chondrites, which we compare our results to below. We compute 2nd order polynomial fits to the data shown in Table 2 over a temperature range of 175-300 K and with a y-intercept fixed at zero (because at low temperatures the heat capacity approaches zero). Among our dataset, we find that two distinct trends emerge, forming two groupings: one with all the iron-nickel meteorites (metallic) and one comprised of all chondritic and achondritic meteorites (non-metallic). The best-fit heat capacity equations (and coefficient uncertainties) for non-metallic and metallic meteorites at different temperatures are given by c s = -0.0033 (± 0.0004)×T 2 + 3.39 (± 0.10)×T and c s = -0.0044 (± 0.0004)×T 2 + 2.84 (± 0.11)×T , respectively. Results are shown in Fig. 3 and the fits from Macke et al. (2019) and Opeil et al. (2020) are included for comparison. Because the Opeil et al. (2020) CM chondrite dataset clearly exhibit higher heat capacities than our fits to other non-metallic meteorites, we use the function c s = -0.0036×T 2 + 3.84×T for Ch and B-type asteroids, which is calculated as 115% of the fit to the non-metallic meteorites. As shown in Fig. 3, this scaled function falls in the range of formulas presented by Opeil et al. (2020).

Mechanical Properties and Emissivity
Young's modulus (E) and Poisson's ratio (ν) are used in the thermal conductivity model to calculate the contact area between grains within a regolith. In particular, these two quantities are used to approximate the deformation both along the axis of an applied force and in the orthogonal dimensions. For M-type asteroids, we use the values of these properties measured for Fe-Ni alloys, which are similar to iron meteorites, by Ledbetter and Reed (1973). Ibrahim (2012) reports Young's and Bulk Modulus (G) for many ordinary and carbonaceous chondrite meteorites. We calculate Poisson's Ratio using the relationship between the three variables: ν = E 2G −1. This equation assumes isotropic material properties and is also used to propagate the reported E and G uncertainties presented in Ibrahim (2012). The values for E-type asteroids are assumed to be the same as used for S-types, as indicated by italics in Table 3.
Emissivity values of meteorites and terrestrial rocks cluster around 0.9-the value that we adopt in this work for all spectral groups with the exception of suspected metal-rich asteroids. Laboratory and modeling efforts by Sih and W.Barlow (2004) place the emissivity of powdered and solid iron roughly in the range of 0.5 to 0.8 for temperatures of 200 to 300 K. Gundlach and Blum (2013) assume an emissivity for metal-rich asteroids of 0.66, which we adopt here.

Monte Carlo Implementation
We compute a thermal conductivity, k obs eff , that is empirically-derived from Γ and assumed values of c s , ρ grain , and a range of porosities. This value is equated to the theoreticallymodeled thermal conductivity calculated using the above procedure, in order to obtain an estimated regolith grain size. Most input variables are taken from a distribution that is generated based on the uncertainty of that parameter given in Table 3 and mentioned above. Most parameters have associated 1σ uncertainties from which we generate Gaussian probability distributions. The porosity is one exception, for which we use a uniform random distribution from 0.276 to 0.876; the lower value is the porosity of close-hexagonal packing scheme for identical spheres and the upper value represents the hypothesized porosity of cometary regolith (i.e., Sunshine et al., 2016, and references therein). For each iteration, a single grain size is produced, ultimately constructing a distributed set of 10,000 grain sizes.
For some trials, the combination of input parameters produced model thermal conductivity values that are incompatible with the observed value. These cases are therefore not included in the final distribution and are instead replaced by additional, successful trials. We note that failed trails are indicative of input parameters that fall at the tail ends of the adopted distributions and do not necessarily imply incorrect assumptions with the models.
Output grain sizes are transformed into the logarithmic scale developed by Krumbein and Aberdeen (1937) which allows for the comparison of sediment sizes across many orders of magnitude. These Krumbien phi-scale grain sizes, d gφ , are referenced from 1 mm. The output grain diameters, d g = 2r g , are transformed into this scale with: Note that smaller (and negative) values on this scale represent larger grain sizes. Output grain size distribution for an object does not necessarily, or even typically, represent a Gaussian distribution. Thus, we report the median of this distribution as the best-fit grain size, and report upper and lower uncertainties by computing the grain sizes that encompass one standard deviation from the median value. The diurnal e-folding thermal skin depth 3 , is also computed and the median and standard deviations from the output distribution are reported.

Caveats and Model Limitations
For some asteroids in our sample, the material properties inferred from spectroscopy, color, or albedo data may not be correct. For example, incorrect assignment of an asteroid as a metallic M-type when it is actually comprised of silicate S-type material is a possible case that could cause non-negligible change in the grain size estimate. To understand the effects of incorrect classification, we modeled the thermal conductivity model for (22) Kalliope (an Mtype) with both metal-rich and S-type properties. Kalliope's low radar albedo may suggest a silicate-rich instead of a metal-dominated surface (Lupishko and Belskaya, 1989;Ockert-Bell et al., 2010;Hardersen et al., 2011). Assuming that Kalliope has S-type properties instead of being metal-rich increased the grain size by ∼8%, which is far below the typical d gφ uncertainty reported in Sec. 3. From this case, we can reasonably assume that uncertainty in the assumed material properties do not significantly contribute any significant bias in grain sizes.
Both of the thermal conductivity models assume homogeneous, monodispersed grain sizes throughout the surface and constant thermophysical properties with depth. This is also true of reported thermal inertia values. Surface processes that sort grain sizes both vertically and spatially are likely present on asteroids (Richardson et al., 2020, and references therein). Thus, this assumption of regolith homogeneity is certainly not the case for any asteroid, and many spacecraft missions have revealed surfaces that are heterogeneous (Fig. 1). Specifically, the Hayabusa mission showed that the surface of Itokawa was consistent with the remotelydetermined thermal inertia (600-800 J m −2 K −1 s −1/2 ), yet still contained regions dominated by fine-grained regolith-counter to interpretation of a high thermal inertia (Müller et al., 2014;Yano et al., 2006). Overall, most asteroid regoliths are expected to be heterogeneous, and the assumption of a single grain size and homogeneous regolith is a practical simplification. However, estimates on the spatial heterogeneity of an asteroid surface are limited due to the fact that we are using disk-integrated (spatially-unresolved) observations.
The recent spacecraft visits to two primitive NEAs, Ryugu and Bennu, unexpectedly revealed surfaces that were different from the initial interpretation of their respective thermal remote observations (see Fig. 1). The surfaces of these objects exhibited surfaces comprised of a significant fraction of large porous boulders that were significantly larger than the grain size estimated from thermal inertia (Grott et al., 2019;Rozitis et al., 2020). Because these boulders have a low thermal conductivity that is similar to a coarse-grained regolith, they exhibit thermal properties indistinguishable from a coarse-grained regolith. Ryan et al. (2020) claimed that the underestimated grain size of Bennu from its thermal inertia can be explained by the inherently low thermal conductivity of Bennu's surface, which resembles CM chondrites and causes non-isothermal effects within the regolith. Grain size estimates are larger when accounting for this non-isothermality, thus can partly explain the large boulders on Bennu and Ryugu.
Because asteroid regoliths are assuredly a mixture of particle sizes the reported grain size estimates herein should be thought of as a thermally-characteristic grain size. Furthermore, it is unknown how this grain size relates to the size distribution of regoloith particles. Presley and Craddock (2006) modeled the thermal conductivity of known granular mixtures and deduced that the modeled grain size is most representative of the larger regolith grainsspecifically the 85th to 95th percentiles-rather than the mean or modal grain size. On the other hand, Ryan et al. (2020) showed that the particle size from thermal inertia corresponds to within 15% of the mean volumetric particle size for various size-frequency distributions. Additionally, because emitted thermal flux is strongly temperature dependent (∝ T 4 ) it is possible that these thermally-characteristic grain sizes are more representative of warmer regolith patches that consist of smaller grains. Another complication could exist due to the fact that warmer areas emit flux at shorter wavelengths as opposed to cooler areas that emit flux at longer wavelengths. Therefore, its possible that the thermal inertia and thermallycharacteristic grain sizes are correlated with the wavelength(s) used in the observation(s). Investigation into the relationship between the grain size and regolith particle distribution is beyond this work and we suggest it as a topic for future studies.
As remarked by Ryan et al. (2020), the apparent thermal inertia will be that of solid rock when the effective particle size is approximately larger than the skin depth, which could be the case for fast rotators with low k grain . If the size range of particles exceeds the skin depth, then the surface is best modeled as a lateral, checkerboard-like mixture of regolith and bare-rock (boulders) (e.g., Bandfield et al., 2011;Rozitis et al., 2020). Despite the caveats, limitations, and uncertainties of the thermal conductivity model, we claim that the grain sizes derived here are appropriate to identify potential trends among asteroid surfaces.

Results and Analysis
We report modeled estimates of grain size and thermal skin depth in Table 4, along with the model input parameters (with associated uncertainties) for each object-Γ, T surf and spectral type. In some cases, the estimated grain sizes exceed the calculated thermal skin depth (see Sec. 3.1). In order to caution readers we mark grain sizes that exceed the lower 1σ skin depth lower limit with a " " in Table 4. We find that the Sakatani et al. (2017) thermal conductivity model is incompatible with some of the thermal inertias in which the model over estimates the thermal conductivity. In such cases, only the grain size from the Gundlach and Blum (2013) model are reported. Similarly, we find that several C-type asteroids have very low reported thermal inertias that were incompatible with the assumed thermal conductivity. In these cases, we re-ran the model with B-type input parameters, which assume a lower value for the thermal conductivity. We mark these objects as "B ", and include them in the assumed B-type group for the post-hoc analyses. Lastly, we have found that the albedos of a few objects in our set that are inconsistent with the spectral type of the reported family. Similar to Masiero et al. (2013), we thus reject the family membership of these objects and indicate this by crossing through the family name in the rightmost column of Table 4.
We now aim to investigate possible explanatory factors and quantify their influence on the regolith grain size for our sub-sample of the asteroid population. We test the hypothesis that grain sizes are negatively correlated with asteroid diameter (indicating impact driven processes) and whether grain size is positively correlated with rotation period (which is caused by thermal fracturing processes). Various multiple linear (hereafter, multi-linear) regression models are fit to the MBAs to identify and characterize the dependencies of these factors (Sec. 3.2). Although the NEAs in this work offer insight into the regolith of very small asteroids, the chaotic nature of changes in their orbital parameters makes it difficult to hold other factors constant-such as thermal environment and impact flux-potentially complicating the multi-linear analysis. We perform a separate analyses on NEAs in Sec. 3.4. In our trend analyses, we use the grain sizes produced by the Gundlach and Blum (2013) thermal conductivity model because the Sakatani et al. (2017) model was unable to provide an estimate for some objects. In general the Sakatani et al. (2017) grain size estimates are consistently larger than those from the Gundlach and Blum (2013) model, but we claim that the trends we are investigating should still be found among the Sakatani et al. (2017) grain size dataset.

Thermal Skin Depth
First we review the results for the thermal skin depth, l s of each asteroid in order to place our grain size estimates in better context. As we have mentioned, asteroid surfaces are comprised of a unknown mixture of regolith and boulders of various sizes, and it is unknown how the thermal inertia (and characteristic grain size) are representative of this mixture.
It is often thought that l s ought to be considered, in effect, as an upper limit to the size of regolith grains that can be distinguished from bedrock and that particles that are larger than this have effective thermal inertias equal to bedrock (or a boulder that is much larger than l s ). Taking this statement to be truth, we can expect that asteroids with larger l s should have lower thermal inertias, on average, than those with smaller l s . However, this expectation is confounded by the fact that rotation period also influences l s (Eq. 11). To aid our understanding of the relationship between l s , P rot , and Γ we place the objects in our dataset into 5 rotation period bins (< 5 hrs, 5−10 hrs, 10−20 hrs, 20−40 hrs, and > 40 hrs). These rotation period bins are assigned different colors in all three panels of Fig. 4. The top panel in this figure shows the distributions of l s for each rotation period bin. The middle and bottom panels show the thermal inertias and calculated grain sizes as a function of skin depth.
When comparing across the rotation period groups in the top panel of Fig. 4 we can see that the two lowest P rot groups have indistinguishable l s distributions. For the other groups, there is a clear increase in the average l s among larger rotation periods. Considering this trend, we claim asteroids with P rot < 5 hrs have larger than expected skin depths. The cause for this is lower thermal inertias among the fast rotators, which is shown in Fig. 5 and quantified in our multi-variate model fits in Sec. 3.2. The question now becomes: is the higher thermal inertia among fast rotators a result of shrinking l s , or caused by an increase in the regolith grain size?
Comparing the thermal inertias and skin depths (middle panel of Fig. 4) within the rotation period bins, some interesting findings emerge. Among asteroids with P rot > 40 hrs, the relationship between thermal inertia and skin depth is roughly what is expected by the proportionality relationship implied by Eq. 11 which is shown as a dashed line with an arbitrarily-set intercept. Most of the objects with l s > 1 m have very long rotation periods and this subset of objects, on average, have thermal inertias lower than 300. Asteroids in other rotation bins exhibit a greater dependence of thermal inertia on skin depth than expected from Eq. 11, which is possible when the porosity decreases for larger thermal inertia values. Therefore, this trend can be explained if these surfaces are covered with more boulders that have lower porosity than the surrounding regolith. It is also possible that this trend is a result of larger grain sizes, which we consider below.
As the relative contact radii between grains decreases for larger grains, the radiative heat transfer becomes dominant over solid heat transfer through grain contacts. In MacLennan and Emery (part I; 2021), we showed that the variation in thermal inertias with temperature are consistent with the theoretical value when radiation is the dominant heat transfer mechanism. Additionally, both the Gundlach and Blum (2013) and Sakatani et al. (2018) models predict that the radiative thermal conductivity component is directly proportional to grain size. Therefore, we can expect that the thermal skin depth increases in proportion to the square-root of the grain size or, equivalently, that the grain size is proportional to the square of the skin depth. The lines showing this proportionality in the bottom panel of Eq. 11 shows that this is indeed the case across our entire dataset. Interestingly, it appears that objects with smaller rotation periods exhibit a stronger relationship between the predicted grain size and skin depth. Furthermore, objects with d g > l s are more likely to be 14 fast rotators as 10 out of 15 of these objects have P rot < 5 hr. In response to the question posed earlier in this section, we posit that the smaller skin depths caused by short rotation periods do indeed have a non-negligible effect on the thermal inertia and resulting regolith grain size. Yet, because most of the grain size estimates for these fast rotators are lower than the skin depth, we claim that this effect is minimal.

Grain Size Factors Model
Multi-linear regression is a method that attempts to model a dependent variable (the grain size, d gφ ) as a linear combination of several independent variables-in our case, diameter and rotation period. The fitted slope, or coefficient, of each independent variable and a y-intercept, along with 1σ uncertainties for each of these parameters is reported. We consider several segmented, or piece-wise, multi-linear regression models with different numbers and combinations of break points for both independent variables (D b eff , and P b rot ), as listed below: These segmented models partition the indicated independent variable into intervals for which a different function (slope) is fit to the data. Adding a break-point creates two lines that form a continuous, piece-wise function in lieu of a single linear fit. In addition to slope estimates of the lines, the location of the break-points are estimated. The diameter and rotation period are transformed by taking the log 10 when used in the multi-linear regression models. This variable transformation is done to best capture the wide variance in these variables which each span more than 2 orders of magnitude. Unlike the thermal inertia multi-linear analysis presented in MacLennan and Emery (2021), we do not consider temperature as an explanatory variable because it is already accounted for in the thermal conductivity model via temperature-dependent heat capacity values (Sec. 2.2) and radiative heat transport.
We compare the multi-linear model fits to one another by taking the adjusted r-squared (r 2 adj ) statistic and Bayesian Information Criterion (BIC). The r 2 adj is a determination of the degree to which the model explains the variance in the dependent variable (i.e., higher values indicate a better fit), while also accounting for the number of predictor parameters, w, in the model. The number of free parameters is calculated from the total number of fitted variables, which increases by 2 when a break-point is added 4 The adjusted r-squared is related to r 2 via: . The BIC is used to indicate which model maximizes the likelihood of matching the data (lower values indicate a higher likelihood), and accounts for the number of model parameters (more parameters increases the score). It can be calculated via BIC = N ln (RSS) − N ln (N ) + w ln (N ) in which RSS is the residual sum of squares (see, e.g. Feigelson and Babu, 2012, for more details). Both of these statistics are shown in Table 5, along with the number of model parameters.
When comparing the results amongst all models, we find that M-4 has the largest r 2 adj and lowest BIC. Statistically speaking, the difference between two model BIC values, ∆BIC, indicate a preference for one over another. Kass and Raftery (1995) state that ∆BIC > 6 indicates a strong preference for the lower BIC. Since M-4 has a ∆BIC = 34 between it and the second lowest BIC value (M-2), we use this as an indication that it is the preferred model. We plot the best-fit and preferred M-4 model in Fig. 5. Black symbols are estimated grain sizes for each object and colored symbols are values from the multi-linear models. Solid bars indicate the 1σ range of the diameter and rotation period breakpoints, respectively, at the top of each panel. The estimated intercepts, linear coefficients (slopes), and break-point(s) (D b eff , P b rot ), along with the associated uncertainties, are listed in Table 6. We remark that the model intercept value represents the predicted d gφ for a hypothetical object with D eff = 1 km and P rot = 1 hr, which is just under 1 m for M-4.

Compositional Effects
Here, we investigate if, and how much, the grain sizes vary by or depend on the surface composition of an asteroid. For this analysis we use spectral group as a proxy for composition. Simply comparing the means of the grain size distributions between the groups is not appropriate, as some of the independent explanatory variables are correlated with grain size: for example, primitive C-complex and P-type bodies are more represented at larger sizes, and the E/Xe-types are largely represented in the lower size range 5 . Instead, we perform post-hoc t-tests between the spectral classes by using the model residuals from M-4. Since the multi-linear analysis doesn't account for the spectral type, any differences in the group residuals can be used to indicate disparities in regolith grain sizes.
We apply Welch's t-test 6 (Welch, 1947) in a series of trials between each possible combination of spectral class, as well as between each class and the entire sample with that particular class removed. The null hypothesis that is tested is that the means between the groups do not differ. We report the p-values of these trials in Table 7: lower p-values indicate a higher probability that the null hypothesis is not supported. Table 7 shows the mean and standard deviation of the model residuals and the number of objects for each group. We analyze the S, C, and E-types both with and without the objects in which their spectral type was inferred by their albedo. Similarly, for B-types we analyze both the confirmed B-types and those marked as B in Table 3 for which we assigned CM meteorite properties based on their low-thermal inertias. Consistent with their low thermal inertias, the B objects have consistently lower grain size residuals compared to confirmed B-types. This may indicate that these objects are not true B-types and may actually be closer in composition to P-types, which have the lowest grain sizes of all the spectral groups included here.
The model residual distributions, grouped by composition, along with the mean model residuals are depicted in Fig. 6. Compared to the mean of the remainder of the sample, there is strong statistical evidence (p < .01) that suspected metal-rich asteroids on average exhibit grain sizes that are nearly twice as large as asteroids of the same diameter and rotation period. We interpret this difference as an indication that one of the thermophysical or material properties of FeNi metal affects the efficiency of regolith breakdown process. Additionally, the average grain sizes of M-types and E-types are statistically indistinguishable from one another (Table 7), and from the suspected metal-rich objects. On the other hand, P-types show a much smaller mean grain size than the remainder of the sample (p < .01). Finally, we note here that S-types and C-types have statistically smaller and larger grain sizes, on average, compared to the rest of the sample. Consistent with this finding, the carbonaceous K-types have coarser-grained regoliths compared to S-types. Potential explanations for these group differences in grain size are discussed in Sec. 4.

Near-Earth Asteroids
Thus far in our analysis we have excluded NEAs and Mars-Crossers (MCs) in order to mitigate the potential influence on regolith evolution caused by widely-varying thermal and impact environments. Because these asteroids exhibit a wide range of orbital parameters (varying in both semimajor axis, a, and eccentricity, e) they are subject to drastically varying external influence, whereas we wish to examine the factors inherent to asteroids themselves. We analyze the grain sizes (Sec. 3) calculated for the NEAs with thermal inertia estimates from other works (N = 21) and those estimated in MacLennan and Emery (2021) (N = 7). With grain sizes for these 28 asteroids, we seek to identify regolith dependencies on orbital factors.
We perform a multi-linear regression model, similar to the one presented in Sec. 3.2, but with different input factors. We do not consider segmented linear fits as the low number of objects in this NEA subset may result in model over-fitting. The independent explanatory variables we include here are the diameter and rotation period transformed into log 10 space. But now we also include the orbital semimajor axis (a), perihelion (q), and aphelion distance (Q) as possible explanatory variables. However, these orbital parameters are not transformed into logarithmic space in the multi-linear model and are left as-is.
The best-fit regression model for NEAs indicated that Q is the only statistically significant variable. Thus, the NEA grain sizes dataset does not appear to show the same trends with the diameter and rotation period that exist among main-belt objects. On further inspection, potential co-linearity between D eff and Q could raise some doubt of the significance of this result. However, employing a multi-linear model with only Q and D eff did not change the significance of either variable. We thus conclude that Q is a better predictor of grain size than diameter and rotation period for NEAs. The grain sizes of NEAs as a function of aphelion, perihelion, diameter, and rotation period are shown in Fig. 7. Symbol colors indicate the spectral type of the objects, with unfilled symbols indicating that the classification was used using the albedo as a proxy. It is interesting to note here that the grain sizes of objects with 0.9 < q < 1.1 au are somewhat higher (on average) than the asteroids outside this range and that asteroids with largest grain sizes (> 10 cm) are found having perihelia interior to Earth's orbit.
An important caveat to note here is the fact that NEA and MC asteroid diameters extend to a lower size range compared to MBAs. Because the overlap in NEA and MBA diameters occurs for 2 km < D eff < 40 km, it is difficult to draw conclusions about the difference between the two populations from these sets of multi-linear models. Grain sizes for NEAs with D eff > 2 km (Fig. 7) are similar to those of similarly sized MBAs and exhibit the same inverse dependency on asteroid diameter. On the other hand, NEAs smaller than 2 km don't show any discernible trend between grain size and asteroid diameter. It is possible that MBAs exhibit a similar lack of trend, but the MBAs in our dataset do not extend to this size range.
Additionally we note that, in general, the comparison between rotation period and grain size for NEAs appears to be consistent with MBAs, with a few exceptions. Asteroids with spin rates near the spin barrier (P rot ≈ 2.12 hr) exhibit a range of grain size spanning 3 orders of magnitude, which mimics the grain size span for the entire Main Belt sample. At the low end of the range is 1950 DA, with two other fast-rotating NEAs exhibiting very large grain sizes. The latter two appear to be very consistent with the inverse trend between rotation period and grain size for main-belt objects. Thus, 1950 DA appears to be an outlier in this respect and we further discuss this point below.

Discussion
The best-fit multi-variate model, M-4, indicates that a break-point value near ∼ 10 km occurs in the relationship between asteroid diameter and regolith grain size. The slope fit to objects smaller than this size shows a strong inverse relationship between asteroid diameter and regolith grain size. On average, the regolith grain size of 10 km bodies is 0.6 mm, with M-4 predicting a grain size of 25 mm for 2 km objects-an increase in d gφ of nearly a factor of fifty, which translates to a power slope of ∼ 6. The M-4 slope fit for for asteroids larger than 10 km is statistically indistinguishable from zero, which suggests no dependence of regolith grain size on asteroid diameter. Interestingly, there appears to be an inverse relationship between the upper limit grain size and asteroid diameter among objects larger than ∼ 80 km. This may indicate a different process or an additional factor that is unaccounted for in our model that influences regolith development on 100-km-scale asteroids, which are through to be primordial bodies.
The preferred multi-linear model predicts that objects with P rot ≈ 5 hrs, on average, have the smallest grain sizes when comparing across all rotation periods. As the rotation period increases, the model predicts that the grain size increases. Objects rotating faster than ∼ 5 hrs show a steep increase in grain size with decreasing rotation period. It is interesting to note that the model predicts that objects with 2.2 hrs and 300 hrs rotation periods should both have d g ∼ 8 mm. These grain size trends with rotation period could indicate at least one regolith evolution processes that depend on the spin rates of asteroids. We explore such processes in Sec. 4.2 and Sec. 4.3. Alternatively, the inverse trend at lower rotation periods could be a result of thermal skin depths that are comparable to the particle sizes, which we have shown this not to be the case in Sec. 3.1.

Meteoroid Impact Breakdown
Repeated impacts of smaller asteroids or meteoroids into the surface of an asteroid will create ejecta-some of which is retained at the surface as regolith. The general relationship between the mass of the ejecta (M ) above a certain velocity (u), as a function of impactor mass and velocity (m & U ) is: where τ , an empirically-derived exponent, is taken as 0.41-the value for sand (Housen and Holsapple, 2011). For reference, the probability distribution of relative velocities between asteroids in the Main Belt has a mode of 4.3 km s −1 and a mean value of 5.3 km s −1 (Farinella and Davis, 1992;Bottke et al., 1994). The velocities of incoming meteoroids are mostly independent of the size of target body (due to the extremely small gravitational attraction of an asteroid), which leaves the relative size of the impactor to the target asteroid as the dominant factor in Eq. 12. Additionally, the mechanical and structural properties of the target asteroid play a role in determining the outcome (Housen and Holsapple, 2003). For example, in targets with low porosity, the ejecta velocities are larger than that of an otherwise identical impact into a porous target (Nakamura and Fujiwara, 1991;Nakamura et al., 1994). Energy from the impact is concentrated near the impact site and partitioned into crushing the surface material (Flynn et al., 2015), instead of ejecting it at high velocities. Holsapple et al. (2002) points out that if the porosity of the target region is greater than 50%, only 10% of the crater mass is ejected, as most of the energy is partitioned into compacting the material. This experimental finding was verified by Cambioni et al. (2021) who showed that the presence of fine-grained regolith was less likely to be found surrounding porous boulders on Bennu's surface. The authors claim that the production of fine-grained regolith in the general asteroid population is frustrated by the presence of porous boulders. If the surfaces of B-types contain porous boulders similar to Bennu's then we might expect a larger grain size compared to S-types, assuming that the latter are comprised of low porosity boulders. Inconsistent with this claim, the grain sizes of S-types and B-types are statistically similar for the asteroids in our sample. This most likely indicates that other factors may be more significant than porosity in determining the production of regolith. Further investigation into this topic should be pursued in future studies.
Using the measured meteorite flux at 1 au, Basilevsky et al. (2013) calculated characteristic lunar boulder survival lifetimes from boulder size-frequency statistics on rims of lunar craters with known ages. They estimated that it takes 25-50 Myr to destroy 50% of a typical 1 m-sized boulder, and after 190-300 Myr, more than 99 % of boulders should be completely destroyed. Basilevsky et al. (2015) extrapolated the lunar timescale to Ceres and Vesta assuming that boulders are broken down exclusively via impact weathering. On these large MBAs the boulder weathering timescale is 3% of the lunar value (i.e., 0.75-1.5 Myr), which is mainly a consequence of the meteorite flux being ∼2 orders of magnitude greater in the Main Belt than at 1 au. In all cases, the remaining fraction of boulders exponentially decreases with time.
This theoretical timescale is significantly shorter than the ages of most asteroid families, although a few young (< 1 Myr) families have been identified (e.g., Table 1 in Nesvorný et al., 2015) in the Main Belt. In principle, asteroids in these young families could theoretically be used to investigate regolith properties of asteroid surfaces that have been reset from an impact. One consideration that must be made is the fact that the small asteroids that exist in these recently-formed families are less likely to retain fine-grain regolith particles because of their low surface gravity. Statistically-speaking, a 1 km asteroid is expected to experience an average of 5 catastrophic breakups every 1 Myr (Holsapple et al., 2002), or once every ∼ 200 kyr. This timescale is shorter than the estimated boulder survival timescale (0.75-1.5 Myr) of Basilevsky et al. (2015), so we expect that an MBA of this size would not survive long enough to have much of a developed regolith. On the other hand, a 10 km MBA will survive, on average, 16 Myr (∼ 80 times longer than a 1 km body). The implication that 10 km asteroids exist long enough to develop a regolith through impact degradation is consistent with our grain size results.
Our results show that the differences in grain size for NEAs and MCs are most likely dependent on their orbital parameters. This relationship could be explained by the drastic change in collisional probability and velocities of NEAs with MBAs in the Main-Belt (Bottke et al., 1993). Because an NEA that enter the Main-Belt is traveling relatively slower, the velocity difference between it and a typical MBA will be larger compared to two MBAs with similar orbits. This would theoretically result in more efficient impact weathering rate for NEAs with Q 2.2 au, compared to MBAs of a similar size. However, there is no obvious difference in grain sizes for NEAs and MBAs of similar object sizes. Thus, we suppose that another process is the reason for the grain size dependency on orbital parameters.
Evidence from meteorite impact breccias-fragmented samples of an asteroid collision (Burbine and Binzel, 2002)-offers insight into to the state of asteroid regolith. These meteorites are formed from the lithification of near-surface material from the heat and pressure of an impact and classified by their texture and presence/absence of clasts (Bischoff et al., 2006). A direct connection to asteroid regolith can be made if they are rich in solarwind gases that are identified isotopically (McKay et al., 1989). Clasts ranging in size from several hundred microns up to many centimeters have been found in HEDs, carbonaceous, ordinary, and enstatite chondrites, and other stony meteorite regolith breccias (Bischoff et al., 2006, and references within). Interestingly, Bischoff et al. (2006) point out that some aubrites are known to contain large enstatite clasts up to 10 cm in size, as well as metal grains up to 1 cm in size. If the aubrite parent body is matched by an E-type spectrum, then these systematically-larger clasts and grains are well-matched to the finding that E-types have larger regolith grains than other spectral types.

Thermal Fatigue Breakdown
The cyclic heating and cooling experienced by asteroid surfaces, as a consequence of time-varying insolation, can result in large spatial temperature gradients within the material. These gradients create heterogeneous thermal stress fields (expansion and contraction) across mineral grain boundaries of a rock . This thermal fatigue process involves the structural weakening of the material from thermal gradients and ultimately results in the breakdown of asteroid regolith. It has been debated whether or not thermal fracturing is a relatively significant weathering mechanism for terrestrial rocks (e.g., Molaro and McKay, 2010, and references within) yet Delbo' et al. (2014) experimentally demonstrated the effectiveness of this process, in vacuo, on a chip of the Murchison (CM2 carbonaceous chondrite) and Sahara 97210 (L/LL 3.2 ordinary chondrite) meteorites. Libourel et al. (2021) performed a follow-up study with these two meteorites and Allende (CV3) where they demonstrated that larger temperature gradients, such as those experienced for low perihelion asteroids, result in greater crack growth rates. Images from the encounter of Bennu by OSIRIS-REx provided widespread evidence of thermal cracking across the surface . Convincing evidence from the Moon and Mars (Ruesch et al., 2020;Eppes et al., 2015) indicates that thermal breakdown is efficient on other solar system bodies. Both Molaro et al. (2017) and Ravaji et al. (2019) posit that the efficiency of thermal fracturing is controlled by the maximum thermal stress experienced by the rock, as opposed to the average, so we focus on this quantity in our discussion. Molaro and Byrne (2012) numerically simulated the amount of internal stress experienced by a rock on Vesta, Mercury, and the Moon's surface. In particular, they estimated the amount of temperature change, per unit time ( ∆T /∆t) and compared to the temperature gradient ∇T within rocks having orientations at different latitudes on the surface. They found that while a rock on Vesta experienced the largest temporal temperature change, it also had the lowest ∇T-leading to a smaller stress field. This case is opposite to that of the Moon and Mercury, in which ∇T and ∆T /∆t were highly correlated. The modeling efforts of Molaro and Byrne (2012) demonstrated that greater thermal stresses were experienced for rocks experiencing faster sunrises (i.e., shorter rotation periods). The stresses were more pronounced for rock surfaces that faced the rising sun. Daytime shadowing, especially when occurring just after local sunrise or before local sunset, was also a major contributor to increasing the temperature gradients. The authors also predict that pre-existing cracks would contribute even more to the crack growth rate and shorten breakdown timescales.
Our asteroid grain size modeling results indicate that objects with P rot = 5 hrs are most likely to exhibit the fine-grained regoliths and coarse-grained regoliths (d g 1 cm) are found among the asteroids with P rot > 100 hrs. This result coincides with the findings of Molaro et al. (2017), El Mir et al. (2019, and Ravaji et al. (2019) as they all predict less-efficient thermal fatigue for objects with longer rotation periods. The rotation period breakpoint in the multilinear fit at ∼ 5 hrs could be indicative of an optimal rotation rate that maximizes thermal fracturing for asteroid surfaces. The Molaro et al. (2017) model do not predict any rotation period that is optimized thermal fracturing efficiency, yet El Mir et al. (2019) and Ravaji et al. (2019) found that breakdown timescales for a 10 cm rock are minimized for rotation periods in the range 10 − 15 hrs. It is not immediately obvious why these works predict a slightly larger optimal rotation rate for thermal breakdown. Future works should seek to reconcile thermal fatigue models with our empirical findings. Molaro et al. (2015) studied the sensitivity of their micro-stress model to changes in the material properties of the rock and found the thermal expansion coefficient and Young's modulus (which describes tensile elasticity) to be the most influential properties, and that thermal conductivity, perhaps counter-intuitively, does not significantly affect the rate of crack growth. The thermal parameter, may intuitively at first seem to be useful in predicting the effectiveness of thermal fatigue on asteroid surfaces, as it can be used as a proxy for the amplitude of the diurnal temperature range. However, its formulation doesn't capture the essential information regarding the heating and cooling rates of a boulder/rock at the surface. For example, regolith on very slow rotators will experience slow sunrises that result in relatively slow rate of temperature change. This behavior may seem paradoxical if one were only to examine the thermal parameter for a hypothetical slow-rotator: as P rot increases, Θ becomes smaller and the diurnal temperature range is maximized. However, the rate of temperature change is relatively small for a hypothetical point on the surface of a slow-rotator, and the insolation changes experienced during long sunrises and sunsets do not result in spatial temperature gradients large enough to cause significant thermal fracturing.
Our results indicate that the grain size of E-type (including assumed E-types), M-type, and suspected metal-rich (Met) asteroids are, on average, ∼ 2 times higher than asteroids of the same diameter and rotation period (Fig. 5). This suggests that either regolith generation is inefficient for these bodies or that smaller regolith grains are preferentially lost (see Sec. 4.3). Considering the former, we consider the idea that compositional differences may lead to less efficient breakdown. For example, the relatively high values of Young's moduli for FeNi metal should result in more efficient breakdown via thermal fatigue  because stiffer materials are less able to accommodate the differential thermal expansion of grains. This idea is supported if we consider the heterogeneous mesosiderite meteorites, which consist of a silicate-metal mixture, as the appropriate analog for metalrich bodies. We suspect that greater thermal stresses surrounding silicate grains will lead to efficient thermal fatigue breakdown. Therefore we suspect that E-types, M-types, and metal-rich bodies are subject to one or more processes that removes smaller particles from their surfaces.
Finally, we remark that thermal fatigue may be a self-limiting process, as long as other processes are insignificant. In this case we claim that an insulating layer of regolith may form over a boulder as a result of thermal fatigue, effectively shielding it from larger temperature gradients. In actuality, impact processes disrupt asteroid surfaces over time and overturn the topmost layers (i.e., impact gardening; Hörz and Schaal, 1977). This may lead to regolith stratification over long periods of time, depending on the impact gardening timescale and efficiency relative to thermal breakdown. Therefore, asteroid regolith may become stratified at mm or sub-mm length scales with progressively smaller grains near the surface.

Regolith Retention & Loss
The discovery and characterization of active asteroids-those which are observed to exhibit mass-loss-is direct evidence that asteroids are subject to processes that remove regolith from their surface. Several potential mechanisms have been proposed and some remain as only potential explanations the activity. Jewitt (2012) reviews such mechanisms, which include impact, centrifugal (inertial), and electrostatic ejection. The OSIRIS-REx cameras directly observed particle ejection from Bennu's surface, which Lauretta et al. (2019) posited were caused by thermal fracturing, volatile dehydration, or meteoroid impacts. Other plausible particle ejection mechanisms such as electrostatic charging were ruled out for Bennu but may act on other objects' surfaces. Depending on the physical properties of an asteroid and its orbital location, all of these processes may or may not be relevant or efficient. We also note that each mechanism may alter asteroid surfaces in different ways because their effectiveness is dependent on the grain size.
In addition to generating regolith, meteoroid impacts may also remove regolith grains from an asteroid surface. As noted above, the energy from a meteoroid impact upon an asteroid surface is partly transferred to individual regolith grains, resulting in a velocity distribution that varies with grain size (Takasawa et al., 2011). As remarked by Lauretta et al. (2019), high energy impacts into asteroid surfaces (large impactor mass and/or velocity) do not produce ejecta, but instead partition their energy into modifying the target material (e.g. Fiege et al., 2019). The smallest regolith particles have higher velocities and are thus the most likely to be ejected. The upper size limit for particle ejection depends on factors such as the impactor mass and target mass (Eq. 12). These variables can take on a wide range of values varied in order to predict the mass fraction of ejecta that is retained as a function of velocity. The active asteroids (596) Scheila (D eff ∼ 113 km) and P/2010 R2 (La Sagra) (D eff ∼ 1.4 km) are suspected to have been subject to impact-induced mass loss due to their single and sudden brightening events (Jewitt, 2012).
The material strength of the target plays a vital role in the ejecta outcome. In general, stronger surface materials produce more ejecta, which can escape the gravitational well of the body more easily than an asteroid comprised of weaker material. In addition, Matsui and Schultz (1984) show that impacts into brittle (low temperature) metal will result in fracture and spallation compared to ductile metal at higher temperature. For high energy impacts, the strength and/or porosity of the object can play a similar role as the surface material does for low energy impacts. We expect that ∼ 100 km objects that exhibit larger regolith grains could be stronger and/or less porous than those with smaller regolith grains. This could explain why the suspected metal-rich asteroids in our sample, all of which are > 10 km, have larger than average grain sizes. Future efforts in mass determination for many asteroids can be used to investigate a potential relationship between the macro-porosity of asteroids and regolith development.
Our results show that, on average, grain sizes for asteroids larger than ∼ 10 km don't show any dependence on object size. On the other hand, the average grain size increases with decreasing asteroid size for objects smaller than 10 km. This size cutoff could be indicative of an abrupt change in the process of regolith breakdown, or alternatively, may indicate a change in object physical properties (i.e., interior strength) that controls the mass of impact ejecta. Using crater scaling laws, Housen et al. (1979) predict that asteroids in the 1-10 km range should not harbor a significant regolith ( 1 mm thick) due to a decrease in both the strength and gravitational field of bodies of this size. Asteroids with higher strengths have an increased ejection velocities, for similar impact scenarios, compared to low strength objects (Nakamura et al., 1994). Strength could be coincident with the material composition, which could explain why metal-rich bodies, which are presumably stronger, have less fine-grained regolith compared to S-complex and C-complex asteroids.
For fast-rotators, the outward centripetal acceleration near their equators is comparable to the downward gravitational acceleration. Regolith grains in these locations therefore exist in a precarious state, and a small perturbation can transfer enough energy to cause the ejection of equatorial material where the effective gravity is near-zero (Guibout and Scheeres, 2003). It is difficult to theoretically predict whether this kind of mass-loss occurs incrementally, on a grain-by-grain basis, or catastrophically, with large portions ejected at a time (Scheeres, 2015). Lauretta et al. (2019) ruled out centrifugal ejection as an ejection mechanism for Bennu and, to date, centrifugal ejection has not been definitively determined for any active asteroid (Jewitt, 2012).
Among small asteroids, the cohesive (e.g., van der Waals) forces between regolith grains smaller than ∼ 1 cm have theoretically (Scheeres et al., 2010) and experimentally (Murdoch et al., 2015) been shown to dominate over gravitational, inertial, and, in some cases, electrostatic forces. Evidence for cohesive forces present on an asteroid surface was found for the 1.3 km NEA (29075) 1950 DA (Rozitis et al., 2014;Gundlach and Blum, 2015). Cohesive forces are inversely dependent on the surface area of the regolith grains (Scheeres et al., 2010). This suggests that a regolith comprised of smaller grains may be stronger than one comprised of larger grains, which has implications for the effectiveness and behavior of potential mass-loss mechanisms. Particles or rocks with low cohesion could be individually ejected, yet inter-particle forces can cause the regolith to be more susceptible to large-scale structural failure. In this scenario planes of weakness preferentially form around massive clumps of grains (Sanchez and Scheeres, 2020). As these grain structures are held in a higher energy state, even a small perturbing force that exceeds the yield stress could trigger a landslide (Scheeres et al., 2010), potentially resulting in a catastrophic loss of fine-grained regolith or, more likely, a mixture of regolith and boulders.
The best-fit multilinear model predicts that faster rotators (with P rot < 5 hrs) have larger regolith grains. Generally speaking, this finding is consistent with thermal breakdown of regolith and/or with centrifugal ejection of smaller grains. Because smaller grains have greater cohesive properties, we slightly favor the former, although it is possible that centrifugal mass loss can affect a wide range of particle sizes. Interestingly, two fast rotators in the sample, (29075) 1950 DA (D eff ≈ 1.3 km, P rot ≈ 2.12) and (3554) Amun (D eff ≈ 2.7 km, P rot ≈ 2.53 hrs), posses regolith characteristics at the extremes of the sample: Amun has a high thermal inertia, suggesting it may be completely devoid of fine-grained regolith, in contrast to 1950 DA (Rozitis et al., 2014). Amun may have always been lacking any significant regolith or it could had been similar to 1950 DA in the past. In the latter case, it was subsequently and periodically spun-up such that centrifugal forces exceeded any cohesion holding the regolith together resulting in the loss of fine-grained regolith held together in massive clumps, as described above. In this scenario, we might expect 1950 DA to lose clumps of fine grained material as its rotation rate incrementally increases and structural failure occurs within the regolith. More thermal inertia estimates of rapid rotators (P rot < 3 hrs) may provide greater detail into the efficiency and probability of centrifugal ejection and the cohesive strength of regolith (Rozitis et al., 2014).
Regolith particles can become electrostatically charged when exposed to ionized solar wind particles, which are primarily comprised of protons (Lee, 1996). Solar UV radiation releases photoelectrons from the surface and a negatively-charged "sheath" builds up on the sunlit side of an asteroid. The electrostatic forces that regolith grains experience in this field can potentially overcome the cohesive and gravitational forces (Hartzell, 2019). For smaller grains, cohesive forces will always dominate and electrostatic forces are not likely to remove them from the surface. The largest grains, on the other hand, are too heavy to be affected by electrostatic forces. However, there is an optimal size-proportional to g 1/4 a -for which electrostatic forces overcome cohesive and gravitational forces, resulting in grain levitation (Hartzell and Scheeres, 2011). The optimal grain size is estimated by Hartzell and Scheeres (2011) to be around 2-4 cm for an Eros-sized object (g a ∼ 5 × 10 −3 m s −2 ) and 10-30 cm for an Itokawa-sized object (g a ∼ 8 × 10 −5 m s −2 ).
If electrostatic forces dominate over cohesive forces, we would then expect a positive correlation between asteroid size and regolith grain size as opposed to an inverse relationship. Although the entire NEA grainsize dataset doesn't show any correlation with asteroid size, the sub-km NEAs in our dataset do exhibit a positive correlation between effective asteroid diameter and grain size that is consistent with this expectation. We suspect that regolith processes in this size range become confounded by multiple competing forces and future modeling could show which process(es) dominate on different bodies.

Conclusions and Future Work
In this work we have presented regolith grain size estimates for 452 asteroids (Sec. 3) and performed a few multi-linear models fits to the dataset and post-hoc tests between compositional groups and asteroid families. From our results we conclude the following: 1) Regolith grain size across our sample of main-belt asteroids is inversely dependent on object size when the effective diameter is below ∼ 10 km (Sec. 3.2). We suspect that this relationship is most likely due to the removal of smaller regolith particles via impact processes due to the lower surface gravity of smaller objects. For MBAs larger than ∼ 10 km, regolith grain size does not exhibit any dependence on asteroid size. Additionally, grain size shows a positive dependence on rotation period for P rot > 5 hrs and an inverse relationship below this value. Nearly all the asteroids with rotation periods above 100 hrs are coarse-grained with grain sizes > 1 cm. 2) We show that potential metal-rich asteroids have consistently larger regolith grains than objects belonging to other spectral types of similar size and rotation period Sec. 3.3. We find that M-types and E-types have coarsergrained regoliths, which may be related to their enstatite content. Finally, the compositionally-primitive P-types exhibit lower-than-average regolith grain sizes (more fine-grained) than any other spectral type.
3) The regolith grain sizes of NEAs, which represent asteroid diameters smaller than the MBAs in our sample, are inversely dependent on aphelion (Sec. 3.4) but not diameter or rotation period. Additionally, the grain sizes of NEAs with aphelia in the Main-Belt are consistent with MBAs with similar aphelia, and the grain sizes of > 2 km NEAs are consistent with that of MBAs at the same sizes. Thus, it is not possible to say if NEAs have different regoliths than MBAs, or if aphelion, rather than size, is a controlling factor of NEA grain sizes. 4) Finally, we infer from our analyses that evidence for both impact weathering and thermal fatigue/fracture exists in the asteroid population. In future work, processes that remove regolith (Sec. 4.3), which may preferentially act on smallest regolith particles or largest boulders, should be accounted for when modeling asteroid regolith evolution.
Although the thermal conductivity equations in this work are useful to estimate grain size from thermal inertia, more effort is needed to improve modeling efforts to better resemble asteroid surfaces and regoliths. In particular, considerations for mixed surfaces comprised of heterogeneous grains, and both solid and porous boulders should be sought (e.g. Cambioni et al., 2019).
Future thermophysical investigations of asteroids should target objects belonging to categories for which there are a relatively low number of thermal inertia estimates. For example, MBAs smaller than ≈ 3 km, super slow rotators (P rot > 100 hrs), and members of very young asteroid families (< 10 Myr). In addition, thermophysical investigation of the regoliths of object groups not studied in this work (e.g., contact binary asteroids, asteroids with satellites, objects with close planetary encounters, low perihelion asteroids, etc.) would further elucidate regolith evolution processes and their respective timescales.                 28.1 ± 7.9 Notes & References. 1 SMASS I; 2 SMASSS II; 3 S3OS2 (Lazzaro et al., 2004); 4 MITHNEOS; 4b Tholen (1984); 5 Gietzen and Lacy (2007) t The reported grain size is larger than the estimated thermal skin depth. a Thermal inertia values are in SI units (J m −2 K −1 s −1/2 ). C Colors from (Tholen, 1984) are used. e The mutual eclipse duration of 4 hr is used to compute thermal skin depth. F Family association (Table 8) is used. O Family outlier, with suspected family crossed out in last column. See text for details. A Object albedo (Fig. 2) is used to infer meteorite connection. S Reflectance spectrum is used. See above for sources. T Large thermal inertia is used to infer a metal-rich surface. U A uniform, rather than a Gaussian, distribution was used for the thermal inertia.       Table 1) and specific heat capacity as a function of temperature (bottom; Table 2). Dashed lines from this work show linear fits in the top panel, and parabola fits through the origin in the bottom panel, with grey regions representing 95% confidence to the fits. The dashed-dotted line through the Campo del Cielo datum in the top panel assumes a slope of 0.3. Note the break in the y-axis in the top panel.