Thermophysical Investigation of Asteroid Surfaces I: Characterization of Thermal Inertia

The thermal inertia of an asteroid is an indicator of the thermophysical properties of the regolith and is determined by the size of grains on the surface. Previous thermophysical modeling studies of asteroids have identified or suggested that object size, rotation period, and heliocentric distance (a proxy for temperature) as important factors that separately influence thermal inertia. In this work we present new thermal inertias for 239 asteroids and model all three factors in a multi-variate model of thermal inertia. Using multi-epoch infrared data of a large (239) set of objects observed by WISE, we derive the size, albedo, thermal inertia, surface roughness, and sense of spin using a thermophysical modelling approach that doesn't require a priori knowledge of an object's shape or spin axis direction. Our thermal inertia results are consistent with previous values from the literature for similarly sized asteroids, and we identify an excess of retrograde rotators among main-belt asteroids<8 km. We then combine our results with thermal inertias from the literature to construct a multi-variate model and quantify the dependency on asteroid diameter, rotation period, and surface temperature. This multi-variate model, which accounts for co-dependencies between the three independent variables, identified asteroid diameter and surface temperature as strong controls on thermal inertia.


Introduction
The thermophysical characterization of regolith-the unconsolidated, heterogeneous, rocky material covering the surface of other planetary bodies-is an important part of understanding the processes and evolution of airless bodies of the Solar System. By comparing thermal observations to thermophysical models, the regoliths of asteroids can be characterized by their thermal inertia (Γ). Thermal inertia is defined as Γ = √ kρc s , where k is the effective thermal conductivity of the regolith, ρ is the bulk density, and c s is the bulk specific heat capacity. Thermophysical models (TPMs) are often used to derive thermal inertia of a body by comparing the observed fluxes to those estimated from the model.
In addition to asteroid size, the rotation period has been suggested as a factor that controls asteroid thermal inertia. Harris and Drube (2016) used a thermal inertia estimator based on the near-Earth Asteroid Thermal Model (NEATM; Harris, 1998) and found a positive correlation between rotation period, P rot , and thermal inertia. This correlation was attributed to deeper penetration of the thermal wave into subsurface material that was higher in thermal conductivity and/or bulk density (caused by smaller porosities). The thermal wave can be expressed in terms of the thermal skin depth, l s , which is the length scale over which the diurnal temperature variation changes by a factor of e ≈ 2.71828: l s = kP rot /2πρc s . On the other hand, Marciniak et al. (2019) used a TPM to derive thermal inertias of slow-rotators and found no correlation between the two variables.
For airless bodies, the effective thermal conductivity is comprised of a solid and radiative component that correspond to different heat transport mechanisms in a regolith (e.g., Piqueux and Christensen, 2009). If radiation is the dominant form of heat propagation then thermal inertia is expected to vary as T 3/2 (Delbo' et al., 2015) because the radiative conductivity is proportional to T 3 (Vasavada et al., 1999). Rozitis et al. (2018) characterized the thermal inertia variation with heliocentric distance (as a proxy for temperature) for three individual asteroids and found a wide range of scaling dependencies for each. In particular, they found that the thermal inertia of two of the studied asteroids had stronger dependencies on heliocentric distance than the T 3/2 scaling law. On the other hand, the thermal inertia of Bennu as measured by the OSIRIS-REx spacecraft showed no evidence of temperature dependence (Rozitis et al., 2020).
Generally speaking, the temperature distribution of a surface is influenced by the thermal inertia. All other factors kept constant, higher thermal inertia surfaces have a smaller difference in sunlit and nighttime temperature, while low thermal inertia surfaces exhibit a greater diurnal temperature differences. As demonstrated in MacLennan and Emery (2019), thermal emission observations taken at pre-and post-opposition (multi-epoch) are sensitive to differences in the temperature distribution and thus can be used to estimate the thermal inertia. This technique is effective at estimating the thermal inertia when no shape or spin axis information about an asteroid is known a priori.
The sense of spin (i.e. retrograde or prograde) for an asteroid can also be estimated using multi-epoch observations (Mueller, 2007;MacLennan and Emery, 2019). Because there exists a time-lag between the period of maximum heating at local noon and when the maximum surface temperature is reached, a morning-afternoon dichotomy is present on surfaces with a non-zero thermal inertia. Thus, the morning and afternoon sides will correspond, respectively, or inversely, to pre-and post-opposition viewing aspects, depending on the object's sense of rotation.
In this work we use multi-epoch observations and the methods of MacLennan and Emery (2019) to derive thermal inertia and size estimated for 239 asteroids. In some cases, we constrain the roughness and object's sense of spin. Comparing our results to the benchmark study of MacLennan and Emery (2019) and other literature works, we assess the ability of this technique to estimate these TPM parameters. We then incorporate diameter, rotation period, and temperature into a unifying multi-factor thermal inertia model that simultaneously accounts for these variables. In a follow-up work, we produce grain size estimates from the thermal inertia values presented here and investigate compositional differences in regolith properties.

Observations & Thermophysical Modeling
Data from the Wide-field Infrared Survey Explorer (WISE) are used for model fitting. Absolute magnitude (H V ) and slope parameter (G V ) from Oszkiewicz et al. (2011), and P rot from the Asteroid Lightcurve Database (ALCDB; Warner et al., 2009) are used as TPM input values for each object (Table 1) along with mean and peak-to-trough fluxes calculated from sparse lightcurve data. The thermophysical modeling approach presented in MacLennan and Emery (2019) is used, as briefly summarized below, and we thus select objects that were observed by WISE at pre-and post-opposition. In MacLennan and Emery (2019) we extracted and used the mean and peak-to-trough flux quantities from thermal light curves via simple geometric averaging and subtracting the maximum and minimum values, respectively. Although those simplistic calculations are useful for dense lightcurve data, they can be problematic when used on sparsely sampled lightcurves for reasons discussed in Sec. 2.2.

Data Description
In 2010, WISE mapped the entire sky at four photometric filters: referred to as W1, W2, W3, and W4 with wavelength centers near 3.4, 4.6, 12, and 22 µm, respectively (Wright et al., 2010). WISE was designed as an astrophysics all-sky mapping mission, but its infrared sensors detected the thermal emission from warm asteroids in the inner solar system. A data-processing enhancement (NEOWISE; Mainzer et al., 2011a) to the nominal pipeline was thus designed and implemented to identify and measure the emission from these solar system objects. Since the initial mapping phase in which all four bands were operational (the cryogenic phase) the WISE telescope operated at shorter wavelengths and was later reactivated (NEOWISE-R). In this work, we only use the cryogenic phase of the mission.
Since WISE does not target moving objects, the asteroids were only observed for a relatively brief (typically less than a couple days) period of time, referred to as an epoch. Each epoch of observations nominally yielded between 10 to 20 individual measurements that were separated by ≈1.6 hr-the orbital period of the spacecraft. NEOWISE flux data are stored at the Infrared Processing and Analysis Center 2 (IPAC) and each detection of a moving solar system object was reported to the Minor Planet Center 3 (MPC), where the information regarding the sky position and time of observation can be retrieved. In downloading the data, we used the MPC observation file to parse the WISE All-Sky Single Exposure (L1b) catalog on IPAC's Infrared Science Archive (IRSA) and select detections acquired within 10 s of that reported to the MPC, with a search cone of 10 . We shift the isophotal wavelengths of the filters and perform a color-correction to the fluxes (Wright et al., 2010) using a spectrum calculated from NEATM temperatures, as per the recommendation of the WISE Explanatory supplement Cutri et al. (2012). Since the criteria used to parse IPAC can potentially return contaminated (i.e., by a background star or galaxy) or unwanted (non-asteroid) infrared sources from the catalog, we employ Peirce's Criterion (Peirce, 1852;Gould, 1855) on the infrared color, W4-W3, as detailed in MacLennan and Emery (2019), to better ensure the inclusion of only uncontaminated observations of asteroids.

Sparse Lightcurve Sampling
Due to the nature of WISE's orbit and survey cadence, a given asteroid will be observed an average of a dozen times during each epoch. This sparse sampling does not allow for the construction of a well-characterized rotational lightcurve. Since observations are taken at irregularly-spaced rotational phases-depending both on the number of observations and the object's rotation period-information may be missing for crucial points of an object's lightcurve, such as the minima and maxima. The WISE telescope orbital cadence may over-sample certain rotational phases, which poses a challenge for extracting scientificallyimportant characteristics such as the mean and peak-to-trough range of the lightcurve. We present here a technique for extracting these parameters from a statistically-scant photometric set, given a priori knowledge of the object's rotation period. We note that applying this approach to the objects in MacLennan and Emery (2019) does not significantly change the results of that work. The formulations below are similar to the analytical solution of a least-squares sinusoidal fit to lightcurve data, with some differences. We show the results of this technique on (91) Aegina in Fig. 1 and report the fluxes computed from this method for all asteroids studied in this work, along with observing circumstances, in Table 2.
First, we step through each possible pair of flux measurement points and compute their average and difference so that for the ith and jth point the mean and range (absolute differ-ence) are m ij and r ij , respectively. The flux uncertainties (δf ) are summed in quadrature, so that the errors in each mean (δm ij ) and range (δr ij ) are given by: Note the factor of 2 associated with the mean, as per the rules of error propagation. Proceeding, we calculate a weighting factor, s ij , based on the separation in rotational phase (normalized to 2π radians) of the two points, φ i − φ j . In this weighting scheme, pairs that sample around the same rotational phase or half a turn (φ i − φ j = 0, 1/2), are given a weight of s ij = 0, and pairs separated by a quarter turn (φ i − φ j = 1/4, 3/4) have s ij = 1, with linear scaling of the weights between these two extremes (top right panel of Fig. 1).
The weighted flux mean (F ) and error (δF ) are then given by: and Pair means are shown in the bottom left panel of Fig. 1, along with the result of applying Eq. 3 to the data. In order to formulate the lightcurve range (♦F ) and error (δ♦F ) we employ a slightly different approach than that used for the mean. For the ith point, we iterate across every combination of differences between points, to select the jth point that which maximizes the range between the two: r ij . Difference pairs that are separated by a quarter-turn of the asteroid are given more weight based off the pair weight, s ij , from above (i.e., the factor (1 − cos(4πs ij )) 2 ): and The factor, 1 − cos(4πs ij ), is used to scale the s ij factor in Eq. 4 in order to create a weight function based off a sinusoid, as opposed to a linear relationship, because we wish to add weights to the error estimation that are appropriate for a rotating non-spherical object. However, this factor is not used in Eq. 5 since doing so would create a penalty for data that do not resemble a sine function, such as for shapes that significantly deviate from an ellipsoid. Employing Eq. 4 is essentially the same as extracting the peak-to-trough range of the the best-fit sinusoid to the original lightcurve. Pair ranges are shown in the bottom right panel of Fig. 1 with the function given in Eq. 5 in red.

TPM Implementation
The TPM and data-fitting approach used here is identical to that presented in MacLennan and Emery (2019) and is summarized briefly here. First, the surface temperatures are modeled across the surface of a spherical object constructed of discrete facets. The one-dimensional heat transfer equation (Fourier's Law) is numerically solved using the estimated insolation (incoming solar radiation) as the energy input. The discrete facets are characterized as planar faces and divided into latitude bins. A diurnal cycle is simulated by rotating the facets about the object's spin axis. Two types of surfaces are modeled: a perfectly smooth surface in which only direct insolation is considered, and a rough surface that is comprised of spherical-section craters, for which direct and multiply-scattered insolation and thermally re-radiated energy from other facets are calculated. Surface roughness is characterized by the mean surface slope (θ; Hapke, 1984), which is varied by differing both the opening angle of the crater (γ) and the proportion of surface area that is covered by those craters (f R ); the latter is implemented when calculating the flux contribution of rough and smooth surfaces.
We use parameterized forms of the energy balance equation and heat diffusion equation (see MacLennan and Emery, 2019, for further details), which reduces the number of TPM variables that are necessary to calculate a unique surface temperature distribution, in order to construct temperature reference tables and reduce the computational time. In this scheme, the necessary information required for rough surface temperature calculation is the bond albedo (A b ), thermal parameter, and sub-solar latitude; whereas the smooth surface only requires the thermal parameter and sub-solar latitude. In Eq. 6, σ 0 is the Stefan-Boltzmann constant, ε B is the bolometric emissivity, and T eq is the sub-solar equilibrium temperature: In the case of smooth surfaces A b is implicitly accounted for in the T eq term and we thus do not need to specify it to run the TPM. In the case of rough surfaces, A b explicitly determines the amount of multiple-scattering within a crater. Thus we run the rough surface for various values of A b , as detailed in the next paragraph. The surface temperatures for both the smooth and rough surface TPMs are stored in reference tables, expressed as T = T /T eq . The smooth-surface TPM was run for 46 values of sub-solar latitude (0 • to 90 • in 2 • increments) and 116 values of the thermal parameter (spaced equally in log 10 space, from 0 to 450) whereas the rough-surface TPM was iterated across 3 values of γ = {45 • , 68 • , 90 • } and run for 46 values of sub-solar latitude (0 • to 90 • in 2 • increments), 116 values of the thermal parameter (spread out in log 10 space, from 0 to 450), and 7 values of A grid = {0, 0. 1, 0.2, 0.3, 0.4, 0.5, 1}. These parameters are chosen to ensure an accuracy within 1% between the surface temperature values interpolated from the grid and those calculated using the exact model parameters.
Surface temperatures calculated for spheres are mapped to prolate ellipsoids (b/c = 1, where a ≥ b ≥ c) using closed-form algebraic expressions (i.e., Appendix B in MacLennan and Emery, 2019) in order to model elongated bodies of differing a/b axis ratio. Fluxes are calculated for the given observing circumstances by interpolation of the flux calculated using the tabulated temperatures. The flux calculated from the interpolated grid are within 1% of the flux calculated by running the TPM with the exact thermophysical and observing parameters. Finally, thermal flux is calculated by a summation of the individual flux contributions from smooth surface (B smooth ) and crater elements (B rough ) and using a grey-body approximation, with spectral emissivity (ε λ = 0.9): e ∠ and e R∠ are the emission angle of the flat facet and crater element. The facet area is a f , and the crater element areas are a R = 2/(m(1 + cos γ)) where m = 40 is the number of crater elements . We note that the flux calculation formula stated in MacLennan and Emery (2019) neglects the latter two parameters for a rough surface and that the one presented here is correct. For rough (crater) elements, v is used to indicate if it is visible (v = 1) or not (v = 0) to the observer, and Λ is a correction factor used to adjust fluxes that deviate from the pre-computed A grid values (for more details, see MacLennan and Emery, 2019). In our data-fitting approach, the shape, spin vector (λ eclip , β eclip ), roughness, and thermal inertia are left as free parameters that we select from a pre-defined sample space. A sphere and prolate ellipsoids with a/b axis ratios of 1.25, 1.75, 2.5, and 3.5 are used. For each of these shapes, we sample 25 predefined thermal inertia values, 3 default roughness (mean surface slope;θ) values, and 235 spin vectors. We search for the best-fit D eff for each combination of these parameters. Each individual value of γ is paired with a value of f R = {1/2, 4/5, 1} that correspond to default mean surface slopes ofθ = {10 • , 29 • , and 58 • }. The thermal inertia points are uniformly distributed in log 10 space from 0 to 3000 J m −2 K −1 s −1/2 , and the spin vectors are spread evenly throughout the celestial sphere, which is achieved by constructing a Fibonacci lattice in spherical coordinates (e.g., Swinbank and Purser, 2006). For each shape/spin vector/θ/Γ combination we use a routine to find the D eff value which minimizes χ 2 . To place confidence limits on each of the fitted parameters, we use the reduced χ 2 statistic (χ 2 = χ 2 /ν) to express the solutions within a 1σ range asχ 2 <χ 2 min (1 + √ 2ν /ν) and consider solutions withχ 2 min < 8 to be acceptable.

Characteristic Temperature Calculation
Estimating the temperature of an asteroid at the time of observation is necessary to perform our multivariate analysis. Calculating a single value for the characteristic surface temperature of an asteroid, which exhibits wide temperature variations across the surface, can be approached in a few different ways. One approach is to rely on the estimation of the sub-solar temperature, based on the theoretical energy balance formulation (Eq. 7). This approach has two problems: the assumptions made in the energy balance equation will often lead to the overestimation of the true sub-solar temperature, and there is a low likelihood that the sub-observer point is close to the sub-solar point.
In order to overcome these possible problems we calculate the color temperature, T c , by independently fitting a blackbody curve (via least-squares minimization) to the asteroid's W3 and W4 thermal fluxes. This approach implicitly accounts for the spatial variation in surface temperatures and explicitly calculated from the data itself, as opposed to using the sub-solar temperature, T ss . We found that a blackbody assumption (ε B = 1) does not introduce uncertainty in the temperature, as any non-zero value would not shift the peak of a blackbody emission curve, which is related to the temperature through Wein's Law. Fitting a blackbody function directly to the WISE dataset is straightforward, but retroactively applying this approach to thermal inertias found in the literature is less so. To estimate T c for asteroids in the literature we calculate the relationship between the NEATM sub-solar temperature T ss 4 (Eq. 9) and T c for our set of asteroids: T ss = 0.777 × T . This best-fit equation, along with the data, is depicted in Fig. 2 by the dotted red line and blue dot-dash lines showing the 1σ uncertainty bounds in the exponent. We invert this equation and use it on the TPM results from previous works, since it is often not the case that a temperature is reported with the thermal inertia. For these objects we use, and assume a beaming parameter of η = 1.1 (the approximate mean for main-belt objects Mainzer et al., 2011b), to compute T ss and then T c .

Results and Analysis
The TPM was run for 239 objects: 3 near-Earth asteroids (NEAs), 2 Mars-crossers (MCs), and 234 main-belt asteorids (MBAs). Table 3 shows the best-fit and 1σ uncertainties for the effective diameter (D eff ), geometric albedo (p V ), thermal inertia (Γ), surface roughness (θ), elongation (a/b; prolate ellipsoid axis ratio), and sense of spin (⇑ for prograde and ⇓ for retrograde) for all 239 objects, with the results of the 21 object from MacLennan and Emery (2019) included at the top. Diameter errors for D eff > 10 km are below 15% of the diameter value, but can be as high as 40% for objects smaller than 10 km. Upper and lower thermal inertia uncertainties are, on average, 180% and 67% of the reported value, respectively. Surface roughness could only be estimated for 97 of the 239 (41%) objects. Informed by the model validation tests of MacLennan and Emery (2019), we use spherical shapes to make an estimate of the sense of spin, which could be unambiguously estimated for all but 17 of the 239 (93%) objects. In some cases, TPM fits only allowed for a lower or upper bound on the surface roughness. Note that objects with TPM fits havingχ 2 min > 8 are marked in Table 3 and should be used with caution.
We combine our sense of spin results with that in MacLennan and Emery (2019) and compare to the spin poles of object shapes that are in the DAMIT 5 database (Database of Asteroid Models from Inversion Techniques; Durech, 2010). In total, there are 101 objects both datasets and 77 of them have sense of spin estimate that agree. If we assume that the DAMIT spin axis estimates have 100% accurate sense of spin, then the TPM has a 76.2% ± 4.3% success rate, based on binomial probability distribution. MacLennan and Emery (2019) demonstrated that the sense of spin success rate is dependent on the thermal inertia, with a success rate of 65 − 80% in the range Γ = 40 − 150 J m −2 K −1 s −1/2 when using spheres. The fact that this agrees with our comparison in this work is encouraging, yet more investigation into model development should be performed in an effort to improve the success rate of constraining the sense of spin using TPMs.
Based on the TPM results, we observe a correlation between the retrograde/prograde ratio and asteroid size. We bin our set of objects by diameter in Fig. 3 and assign the uncertainty (shown as vertical lines) of the bins to be the number of objects with indeterminate spins in that diameter bin, or the ∼ 76% success rate based on our check with DAMIT spins-whichever is larger. The horizontal lines that transect some of the spin/diameter bins indicate the number of NEAs in that bin. The fraction of prograde to retrograde rotators in most of the bins in our sample are statistically-indistinguishable. Only asteroids with D eff < 8 km show a statistically-significant excess of retrograde rotators, which we discuss in Sec. 4.
The diameter estimates from NEATM fits presented by the WISE team (i.e., Mainzer et al., 2011b;Masiero et al., 2011) are reported for each sighting, or epoch, which can have pre-or post-opposition geometry. Because the NEATM assumes a spherical shape, it is most-useful to compare the volume-equivalent, effective diameters of the ellipsoid to their values. We present a comparison of these diameter pairs (D WISE eff ) to the diameter values (one per object; D TPM eff ) obtained here, and plot them (colored by observing geometry) in Fig. 4. There is a general agreement of within ± 15% between the two datasets, with a few important notes. Firstly, for objects ∼30 km and above, our TPM diameters are slightly higher than the NEATM model estimates of the WISE team. This discrepancy is likely due to the inherent model differences between our TPM approach and the NEATM used by the WISE team. Secondly, objects smaller than 20 km exhibit, on average, 5% lower diameters from our TPM analysis than from the WISE NEATM analysis. Lastly, we highlight an interesting trend seen in Fig. 4 for different observing geometries: pre-opposition (upright triangles) NEATM diameters are more similar to the TPM-derived diameters for objects smaller than ∼ 8 km while the post-opposition (downward triangles) diameters remain consistently offset from our TPM diameters at smaller sizes. From this result we can conclude that the majority (over 50%) of small diameter asteroids are retrograde rotators-which serves as an independent check on our sense of spin results (Sec. 3).
A handful of asteroids with thermal inertias presented here have previous estimates from the works of Hanus et al. (2018), Marciniak et al. (2019), and Pravec et al. (2019). We depict all these estimates in Fig. 5. In several instances, two thermal inertias were reported because two shape/spin solutions were used, for which we show both values. In nearly all cases there is good agreement between our estimate and the previously-reported value (i.e., the error bars overlap). Only (1741) Giclas shows a significant difference between our estimate and the previous estimates )-our estimate is smaller by around a factor of three and there is no overlap at the 1σ level. We note that our retrograde sense of spin estimate for Giclas is opposite to that of the prograde shape solution in Pravec et al. (2019), and our roughness estimate is much higher (≈ 58 • compared to 38.8 • ). If we were to only use the prograde solutions from out fitting, our estimate would not change, but if we consider higher roughness values that have a higherχ 2 min then the uncertainty in our estimate would overlap with Pravec et al. (2019). This difference in thermal inertia may most likely be caused by the ellipsoidal shape assumption used in TPM fitting.
In addition to comparing the diameters and themal inertias of individual objects from our dataset to the findings from the WISE team, we compare thermal inertia results of 220 asteroids from previous TPM works. Combined with the results from MacLennan and Emery (2019), we present thermal inertia estimates for 250 asteroids (19 of which have previous determinations in other works), an approximate doubling over the tally of literature values (Table 4)-mostly in the 5-50 km size range. We highlight previous authors and works that have presented thermal inertia estimates for 5 objects or more, notably Alí-Lagoa et al.  Hanus et al. (2015) and Hanus et al. (2018) have collectively modeled dozens of asteroids that were observed by WISE. Marciniak et al. (2018) and Marciniak et al. (2019) specifically targeted asteroids with longer rotation periods (P rot > 24 hrs); a group of objects that have lacked thermal inertia estimations. We refer to object thermal inertias presented in papers with less than 5 objects as "miscellaneous literature".

Multivariate Regression Model
We implement a forward stepwise multivariate regression model (Draper and Smith, 1998) on the thermal inertias presented here and in previous works (Table 3 and Table 4) in order to characterize the controlling factors. The independent factors in this model are color temperature (T c )-an approximation of the surface temperature; Sec. 2.4-object diameter (D eff ) and rotation period (P rot ). All variables, including Γ are transformed into log 10 space when included in the model. We use the inverse of the uncertainty in Γ as a weighting factor for each object in the model. The forward stepwise regression algorithm permits a factor to enter the model when the relationship with the dependent variable is statistically-significant (p < .05).
The regression model selected all input variables as statistically-significant explanatory variables. The equation is given by, log 10 Γ = W + X(log 10 T c ) + Y (log 10 D eff ) + Z(log 10 P rot ), with best-fit intercept and coefficient values: W = −4.65 ± 0.70, X = 2.74 ± 0.29, Y = −0.17 ± 0.03, and Z = 0.12 ± 0.05. This best-fit model and data are shown in Fig. 6 Previous studies quantified the thermal inertia dependence on diameter (e.g. Delbo' and Tanga, 2009), rotation period (Harris and Drube, 2016), and heliocentric distance (as a proxy for temperature; Rozitis et al., 2018) seperately, but no study to-date has attempted to simultaneously account for the effect of all three of these variables on thermal inertia. Performing a multivariate regression, as we have done here, accounts for confounding effects between variables, such as the codependency between diameter and surface temperature. This is particularly important because the smallest objects tend to be observed at smaller heliocentric distances and thus have larger surface temperatures.

Noteworthy Objects
Notably, two objects in this study have higher estimated thermal inertias than any other asteroid to-date: (3554) Amun. Discovered in 1986, this Venus-crossing, Aten NEA has an estimated size of D eff = 2.71 ± 0.02 km. Its rotation period of 2.53 hr places it close to the theoretical spin barrier limit, and near-infrared reflectance observations show a red and featureless spectrum yielding an ambiguous classification as a X-or D-type (Thomas et al., 2014). Our moderate albedo (p V ≈ 0.241) estimate and its very high thermal inertia are highly suggestive of a metal-rich surface, which may help explain the thermal inertia of ∼ 1400 J m −2 K −1 s −1/2 estimated here. Additionally, our low roughness estimate is interesting to note, as it suggests a surface that is relatively smooth at the cm scale (i.e., on the order of l s ).
(5604) 1992 FE. This V-type NEA is also an Aten and has been flagged as a Potentially Hazardous Asteroid (PHA) by the MPC. This sub-kilometer object has a very high thermal inertia, but with large error bars: 1100 +2200 −600 J m −2 K −1 s −1/2 . Its high optical albedo and radar circular polarization ratio 6 are consistent with its V-type taxonomic classification and having surface properties similar to Vesta (Benner et al., 2008).

Discussion
Our results show a slight excess of prograde spins at larger sizes ( Fig. 3), which is generally consistent with previous findings of spin vector distributions estimated from lightcurve inversion methods (Kryszczynska et al., 2007;Hanus et al., 2011;Durech et al., 2016). The estimated success rate of our sense of spin determinations places some uncertainty on this claim, however. We can be most confident about an prograde/retrograde difference in the > 60 km size bin, which is consistent with the findings of Hanus et al. (2011);Durech et al. (2016). But the large, overlapping uncertainties in the 16-30 km and 30-60 km places doubt on any claim of excess prograde rotators. The prograde excess for these large objects is likely a remnant of the primordial spins of large protoplanets due to the accretion direction of pebbles into planetesimals (Johansen and Lacerda, 2010).
Our results show a curious overabundance of small (< 8 km) retrograde rotators. An overabundance of retrograde rotators among NEAs was presented by Spina et al. (2004), with the cause attributed to a dynamical selection effect: retrograde MBAs are more likely to feed into resonances, via the Yarkovsky effect, that alter their orbits into near-Earth space (Bottke et al., 2002). Properly investigating and explaining this result is beyond the scope of this work, but we suspect that modeling of YORP spin obliquity evolution (e.g., Vokrouhlický et al., 2003) and/or the spin alteration due to collisions (Ševeček et al., 2019) should be used to investigate this topic. Yet, the MBAs studied here have not yet been subjected to this dynamical selection effect. We also note that, because we only consider asteroids with rotation periods in the ALCDB, our object set is subject to the observational biases inherent in the determination of rotation periods. This includes, but is not limited to, the skew of known rotation periods to less than Earth's rotation period and object shapes that depart from spherical shapes.
The multivariate model of asteroid thermal inertia indicates that temperature is a strong controlling factor (p < .001). The best-fit coefficient in Eq. 10 can be written as the proportionality: Γ ∝ T 2.74±0.29 c . Because surface temperatures generally scale with the inverse square of heliocentric distance (Eq. 7), it follows that Γ ∝ r α . Our results can be expressed in terms of this proportionality by using T X c = R −2α au , which gives: α = −1.37 ± 0.14. If only the radiative component of thermal conductivity on thermal inertia is considered, the expected coefficient would be α = −0.75. Rozitis et al. (2018) calculated α for three objects ranging from −2.5 to −1, with each object having a different best-fit α. Our result is remarkably consistent with all three objects (see Fig. 8 in Rozitis et al., 2018), although the asteroids studied in that work exhibited vastly different thermal inertia dependence on heliocentric distance. Similar to Harris and Drube (2016), Rozitis et al. (2018) suggested that increased solar heating would allow the thermal wave to sample higher thermal inertia material in the sub-surface due to an increase in l s . This would result in an increase in thermal inertia for warmer objects if the thermal conductivity and/or bulk density increases with depth-as is the case for the Moon (Keihm and Langseth, 1973).
In the 200-350 K range-the approximate range for most asteroids-the heat capacity is also temperature-dependent (c s ∝ T ; Opeil et al., 2012;Macke et al., 2019) and should also contribute to the thermal inertia temperature-dependence. The overall dependence of thermal inertia on temperature should be stronger than that predicted using only the radiative component of thermal conductivity, namely Γ ∝ T 2 . When combining the temperaturedependence of the radiative component of thermal conductivity and heat capacity together into the thermal inertia dependency of temperature, the expected relationship is still weaker than the observed dependence presented here. Our multivariate regression model for thermal inertia also selected the diameter (p < .001) and rotation period (p = .011) as statistically significant factors. The trend of increasing thermal inertia with smaller asteroid size was established by Delbo' et al. (2007) and has been supported by subsequent works that increase the overall number of thermal inertia estimates (Delbo' et al., 2015). Delbo' and Tanga (2009) found a power-law exponent of −0.21 ± 0.4 between diameter and thermal inertia for NEAs and MBAs with sizes < 100 km. The Delbo' and Tanga (2009) value is consistent with, but somewhat smaller than, our value of −0.17 ± 0.03, implying a stronger relationship. It is not unexpected that our estimate is larger because we account for the temperature-dependency. For example, Rozitis et al. (2018) found that the diameter and temperature power-law exponents are inversely correlated, and that α = −1.32 corresponds to a power-law diameter exponent of −0.18, which is consistent with our estimate of −0.17 ± 0.03.
Whereas the dependency on diameter is statistically robust, our findings show that the relationship between thermal inertia and rotation period is barely distinguishable from a slope of zero (Z = 0.12±0.05). In the work of Harris and Drube (2016), who used a NEATMbased thermal inertia estimator, found a significant correlation between thermal inertia and rotation period for asteroids with rotation periods spanning 2-200 hrs. However, the works of Marciniak et al. (2018Marciniak et al. ( , 2019 found an abundance of low-Γ slow-rotators (P rot > 12 hrs) using a TPM that explicitly accounts for thermal inertia. Considering the results in this paper and from these previous works we claim that the relationship between thermal inertia and rotation period, if present, is very weak.
Future thermophysical modeling efforts should target more slow rotators to better characterize their thermal inertia and understand the its relationship (or lack thereof) with asteroid rotation period. Higher thermal inertias could be indicative of the increase in thermal conductivity (or bulk density due to compression) for objects with large l s values. These objects with large l s , which include asteroids with high surface temperatures , can be used to investigate possible changes in regolith properties as a function of depth (i.e., Harris and Drube, 2016).

Conclusions and Follow-Up Work
In this work, we applied the method of MacLennan and Emery (2019) to WISE multiepoch observations in order to estimate the effective diameter, geometric albedo, thermal inertia, and surface roughness for 239 asteroids (Table 3). Additionally, we report the shape and sense of spin for a large fraction of these objects Sec. 3. Our thermal inertia estimates are consistent with previous literature values for individual objects (Fig. 5) and for objects with similar size and rotation period. From our results, we conclude that surface temperature, asteroid size (inverse relationship), and rotation period are controls of thermal inertia of asteroids. We find that the relationship between thermal inertia and size is present, but less pronounced than suggested in previous works that do not also consider the influence of temperature (Sec. 3.1). The temperature dependence (Γ ∝ T 2.74±0.29 c ) is larger than the theoretical prediction of Γ ∝ T 1.5 if only the temperature-dependence of the radiative component of thermal conductivity is considered, and of Γ ∝ T 2 if the temperature-dependence of heat capacity is additionally considered. Instead, this relationship between thermal inertia and temperature is consistent with temperature-dependency of both the heat capacity and thermal conductivity (Sec. 4). The thermal inertia dependence on object rotation period is weak and increased statistics of slow-rotator thermal inertias in the future could either support or negate this finding.
In a follow-up work, we will utilize a thermal conductivity model to estimate characteristic grain sizes for each object in this thermal inertia dataset. These grain sizes will then be used to investigate plausible regolith development mechanisms such as impact erosion and thermal fatigue cycling. We will then run the grain sizes through a multi-variate regression model similar to that performed here in order to explore the controlling factors of regolith evolution on asteroids.     Masiero et al. (2011) andMainzer et al. (2011b) to our reported TPM values. We plot the difference between the individual pre-and postopposition diameters of the WISE team and our TPM diameter as a function of the diameter from our TPM. Purple, upward-facing and green, downward-facing triangles are data collected at pre-and post-opposition, respectively. the green line shows a running mean of the post-opposition data, the purple line shows a running mean of the pre-opposition data, and the black line shows a running mean of the relative diameter difference for all objects.  Table 4) for individual objects to estimates in this work (filled circles; Table 3). All objects with previous estimates are from Hanus et al. (2015), except (538) Fredicke  and (1741)   : Thermal inertia dependence on surface (color) temperature, object diameter, and rotation period. The black lines show the best-fit multivariate regression model to the data (dashed) and 1σ uncertainty (dotted). Red points are objects from this work (Maclennan'21; Table 3) and MacLennan and Emery (2019). Other colors indicate the source from other works (Table 4), as indicated in the upper right panel, where misc. literature refers to papers that present less than 5 thermal inertias.   Note. † Indicates a Prot value that has been truncated from the reported value at four decimal places.                              Object