Thermophysical Investigation of Asteroid Surfaces. I. Characterization of Thermal Inertia

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 ${\rm{\Gamma }}=\sqrt{k\rho {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 the thermal inertia of a body by comparing the observed fluxes to those estimated from the model.

The largest asteroids (ranging from ∼500 to 1000 km) in the solar system—i.e., (1) Ceres, (2) Pallas, and (4) Vesta—exhibit global ⁵ thermal inertias (Mueller & Lagerros 1998; Capria et al. 2014; Alí-Lagoa et al. 2020; Rognini et al. 2020) that are comparable to that of the Moon (≈50 J m⁻² K⁻¹ s^−1/2; Wesselink 1948; Cremers 1975; Bandfield et al. 2011; Hayne et al. 2017). On the other hand, smaller asteroids such as (433) Eros, (162173) Ryugu, (101955) Bennu, and (25143) Itokawa have approximate, estimated thermal inertias of, respectively 150, 225, 350, and 700 J m⁻² K⁻¹ s^−1/2 (Mueller 2007; Dellagiustina et al. 2019; Shimaki et al. 2020). The observed correlation among asteroids of thermal inertia with size (Delbo' et al. 2015; Hanuš et al. 2018; Rozitis et al. 2018) suggests that only larger asteroids harbor a fine-grained regolith, whereas smaller asteroid regoliths comprise coarse sand-sized particles and/or a higher fraction of blocky material. This general trend suggests that asteroid size is a large factor in determining the state of asteroid regolith.

In addition to asteroid size, the rotation period has been suggested as a factor that controls asteroid thermal inertia. Harris & 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 the 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}=\sqrt{{{kP}}_{\mathrm{rot}}/2\pi \rho {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 comprises a solid and radiative component that corresponds to different heat transport mechanisms in a regolith (e.g., Piqueux & 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³ (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 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 greater diurnal temperature differences. As demonstrated in MacLennan & Emery (2019), thermal emission observations taken at pre- and post-opposition (multiepoch) 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 multiepoch observations (Mueller 2007; MacLennan & 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 nonzero 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 multiepoch observations and the methods of MacLennan & Emery (2019) to derive thermal inertia and size estimated for 239 asteroids. In some cases, we constrain the roughness and the object's sense of spin. Comparing our results to the benchmark study of MacLennan & Emery (2019) and other works in the literature, we assess the ability of this technique to estimate these TPM parameters. We then incorporate diameter, rotation period, and temperature into a unifying multifactor 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.

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 light-curve data. The thermophysical modeling approach presented in MacLennan & 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 & 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 light-curve data, they can be problematic when used on sparsely sampled light curves for reasons discussed in Section 2.2.

2.1. Data Description

2.2. Sparse Light-curve 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 light curve. Because observations are taken at irregularly spaced rotational phases—depending on both the number of observations and the object's rotation period—information may be missing for crucial points of an object's light curve, such as the minima and maxima. The WISE telescope orbital cadence may oversample certain rotational phases, which poses a challenge for extracting scientifically important characteristics such as the mean and peak-to-trough range of the light curve. 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 & 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 light-curve data, with some differences. We show the results of this technique on (91) Aegina in Figure 1 and report the fluxes computed from this method for all asteroids studied in this work, along with observing circumstances, in Table 2.

Zoom In Zoom Out Reset image size

Table 2. WISE Observation Circumstances and Fluxes

Object	UT Date a	Δtobs b	Nc	RAU d	ΔAU e	α⊙ (°) f	$\overline{{F}_{{\rm{W}}3}}$ g	⧫FW3 h	$\overline{{F}_{{\rm{W}}4}}$ g	⧫FW4 h
	17 Jan 2010	1.125	12	2.604	2.376	22.19	8604 ± 1238	2721 ± 2350	12920 ± 260	3602 ± 442
(91) Aegina
	5 Jul 2010	1.258	12	2.772	2.491	−21.43	5539 ± 797	1289 ± 1590	10360 ± 220	1922 ± 280
	30 Jan 2010	1.125	12	2.860	2.694	20.14	599.0 ± 6.9	155.2 ± 14.7	1388 ± 27	363.3 ± 55.6
(155) Scylla
	16 Jul 2010	3.902	20	3.252	2.999	−18.12	299.1 ± 3.7	92.30 ± 7.06	858.8 ± 19.0	241.3 +36.9
	20 Jan 2010	0.992	10	3.278	3.104	17.46	849.2 ± 9.2	152.4 ± 18.5	2411 ± 45	341.6 ± 73.0
(271) Penthesilea
	6 Jul 2010	1.258	15	3.306	3.056	−17.83	789.9 ± 9.3	169.0 ± 18.8	2296 ± 41	493.3 ± 76.6
	31 Jan 2010	1.125	11	3.122	2.977	18.39	198.6 ± 2.6	22.02 ± 5.41	547.7 ± 12.6	46.63 ± 24.87
(295) Theresia
	22 Jul 2010	0.992	12	3.257	3.023	−18.12	155.2 ± 2.2	16.47 ± 4.28	457.0 ± 10.0	46.24 ± 19.72

Notes. All mean flux and range values are in units of mJy = 10⁻²⁹ Wm⁻² Hz⁻¹.

^a UT date of the first observation. ^b Time spanned by observations (days). ^c Number of observations used. ^d Mean heliocentric distance. ^e Mean WISE-centric distance. ^f Mean solar phase angle. Positive and negative values indicate post- and pre-opposition observing geometry, respectively. ^g Light-curve-averaged mean flux. ^h Photometric range of the light curve.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

First, we step through each possible pair of flux measurement points and compute their average and difference so that for the ith and jth points, the mean and range (absolute difference) 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

$$\begin{eqnarray}&&2\delta {m}_{{ij}}=\delta {r}_{{ij}}=\sqrt{\delta {f}_{i}^{2}+\delta {f}_{j}^{2}}.\end{eqnarray} \tag{ 1 }$$

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 Figure 1).

The weighted flux mean ($\overline{F}$) and error ($\delta \overline{F}$) are then given by

$$\begin{eqnarray}&&\overline{F}=\displaystyle \frac{\sum _{i}\sum _{j\gt i}{m}_{{ij}}\delta {{m}_{{ij}}}^{-2}{s}_{{ij}}^{2}}{\sum _{i}\sum _{j\gt i}\delta {{m}_{{ij}}}^{-2}{s}_{{ij}}^{2}}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\delta \overline{F}=\sqrt{{\sum }_{i}{\sum }_{j\gt i}\delta {{m}_{{ij}}}^{2}{s}_{{ij}}^{2}/{\sum }_{i}{\sum }_{j\gt i}{s}_{{ij}}^{2}}.\end{eqnarray} \tag{ 3 }$$

Pair means are shown in the bottom left panel of Figure 1, along with the result of applying Equation (3) to the data.

In order to formulate the light-curve 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 maximizes the range between the two: ${r}_{\widehat{{ij}}}$. Difference pairs that are separated by a quarter-turn of the asteroid are given more weight based on the pair weight, s_ij , from above (i.e., the factor ${\left(1-\cos (4\pi {s}_{\widehat{{ij}}})\right)}^{2}$):

$$\begin{eqnarray}&&\unicode{x029EB}F=2\ \displaystyle \frac{\sum {r}_{\widehat{{ij}}}\delta {r}_{\widehat{{ij}}}^{-2}(1-\cos (4\pi {s}_{\widehat{{ij}}}))}{\sum \delta {r}_{\widehat{{ij}}}^{-2}{\left(1-\cos (4\pi {s}_{\widehat{{ij}}})\right)}^{2}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\delta \unicode{x029EB}F=\sqrt{\sum \delta {r}_{\widehat{{ij}}}^{2}{s}_{\widehat{{ij}}}^{2}/\sum {s}_{\widehat{{ij}}}^{2}}.\end{eqnarray} \tag{ 5 }$$

The factor, $1-\cos (4\pi {s}_{\widehat{{ij}}})$, is used to scale the ${s}_{\widehat{{ij}}}$ factor in Equation (4) in order to create a weight function based on a sinusoid, as opposed to a linear relationship, because we wish to add weights to the error estimation that are appropriate for a rotating nonspherical object. However, this factor is not used in Equation (5) as 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 Equation (4) is essentially the same as extracting the peak-to-trough range of the best-fit sinusoid to the original light curve. Pair ranges are shown in the bottom right panel of Figure 1 with the function given in Equation (5) in red.

2.3. TPM Implementation

The TPM and data-fitting approach used here is identical to that presented in MacLennan & 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 comprising spherical-section craters, for which direct and multiply scattered insolation and thermally reradiated energy from other facets are calculated. Surface roughness is characterized by the mean surface slope ($\bar{\theta };$ 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 & 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,

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{\rm{\Gamma }}}{{\varepsilon }_{B}{\sigma }_{0}{T}_{\mathrm{eq}}^{3}}\sqrt{\displaystyle \frac{2\pi }{{P}_{\mathrm{rot}}}},\end{eqnarray} \tag{ 6 }$$

and subsolar latitude, whereas the smooth surface only requires the thermal parameter and subsolar latitude. In Equation (6), σ₀ is the Stefan–Boltzmann constant, ε_B is the bolometric emissivity, and T_eq is the subsolar equilibrium temperature:

$$\begin{eqnarray}&&{\varepsilon }_{B}{\sigma }_{0}{T}_{\mathrm{eq}}^{4}=\displaystyle \frac{{S}_{\odot }(1-A)}{{R}_{\,\mathrm{au}}^{2}}.\end{eqnarray} \tag{ 7 }$$

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^{\prime} =T/{T}_{\mathrm{eq}}$.

The smooth-surface TPM was run for 46 values of subsolar latitude (0◦–90◦ in 2◦ increments) and 116 values of the thermal parameter (spaced equally in log₁₀ 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 subsolar latitude (0◦–90◦ in 2◦ increments), 116 values of the thermal parameter (spread out in log₁₀ 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 & 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 the smooth-surface (B_smooth) and crater elements (B_rough) and using a graybody approximation, with spectral emissivity (ε_λ = 0.9):

$$\begin{eqnarray}&&\begin{array}{l}F(\lambda )=\displaystyle \frac{{\varepsilon }_{\lambda }}{{{\rm{\Delta }}}_{\,\mathrm{au}}^{2}}\sum _{S}\{(1-{f}_{R}){B}_{\mathrm{smooth}}(\lambda ,T)\cos (e)\\ \quad +{f}_{R}{a}_{R}v{\rm{\Lambda }}{B}_{\mathrm{rough}}(\lambda ,T)\cos ({e}_{R})\}{a}_{f},\end{array}\end{eqnarray} \tag{ 8 }$$

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 \gamma ))$, where m = 40 is the number of crater elements (Emery et al. 2014). We note that the flux calculation formula stated in MacLennan & Emery (2019) ignores 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 & 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 predefined 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; $\bar{\theta }$) 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 corresponds to default mean surface slopes of $\bar{\theta }=\{10^\circ ,29^\circ ,$ and 58°}. The thermal inertia points are uniformly distributed in log₁₀ space from 0 to 3000 J m⁻² K⁻¹ 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 & Purser 2006). For each shape/spin vector/$\bar{\theta }$/Γ combination we use a routine to find the D_eff value that minimizes χ². To place confidence limits on each of the fitted parameters, we use the reduced χ² statistic $\left({\tilde{\chi }}^{2}=\tfrac{{\chi }^{2}}{\nu }\right)$ to express the solutions within a 1σ range as ${\tilde{\chi }}^{2}\lt {\tilde{\chi }}_{\min }^{2}\left(1+\tfrac{\sqrt{2\nu }}{\nu }\right)$ and consider solutions with ${\tilde{\chi }}_{\min }^{2}\lt 8$ to be acceptable.

2.4. 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 subsolar temperature, based on the theoretical energy balance formulation (Equation (7)). This approach has two problems: the assumptions made in the energy balance equation will often lead to the overestimation of the true subsolar temperature, and there is a low likelihood that the subobserver point is close to the subsolar 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 calculates from the data itself, as opposed to using the subsolar temperature, T_ss. We found that a blackbody assumption (ε_B = 1) does not introduce uncertainty in the temperature, as any nonzero 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 data set 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 subsolar temperature T_ss ⁸ (Equation (9)) and T_c for our set of asteroids: ${T}_{\mathrm{ss}}\,=0.777\times {T}_{c}^{1.063(\pm 0.005)}$. This best-fit equation, along with the data, is depicted in Figure 2 by the dotted red line and blue dotted–dashed lines showing the 1σ uncertainty bounds in the exponent. We invert this equation and use it on the TPM results from previous works because it is often not the case that a temperature is reported with the thermal inertia. For these objects we use,

$$\begin{eqnarray}&&{T}_{\mathrm{ss}}={\left[\displaystyle \frac{{S}_{\odot }(1-A)}{\eta {\varepsilon }_{B}{\sigma }_{0}{R}_{\,\mathrm{au}}^{2}}\right]}^{1/4},\end{eqnarray} \tag{ 9 }$$

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 .

Zoom In Zoom Out Reset image size

The TPM was run for 239 objects: 3 near-Earth asteroids (NEAs), 2 Mars-crossers (MCs), and 234 main-belt asteroids (MBAs). Table 3 shows the best-fit and 1σ uncertainties for the effective diameter (D_eff), geometric albedo (p_V ), thermal inertia (Γ), surface roughness ($\bar{\theta }$), 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 & 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 & 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 ${\tilde{\chi }}_{\min }^{2}\gt 8$ are marked in Table 3 and should be used with caution.

Table 3. TPM Results

Object	Deff (km)	pV	Γ a	Tc (K)	$\bar{\theta }(^\circ )$	a/b b	Spin c
(91) Aegina	101.4 ± 13.85	${0.052}_{-0.007}^{+0.007}$	${19}_{-19}^{+31}$	250 ± 16	17 ± 12	1.08 ± 0.19	⇑
(155) Scylla	38.41 ± 0.54	${0.054}_{-0.003}^{+0.002}$	${16}_{-5}^{+15}$	219 ± 12	34 ± 22	1.51 ± 0.27	⇑
(271) Penthesilea	65.05 ± 2.30	${0.054}_{-0.003}^{+0.003}$	${16}_{-16}^{+33}$	210 ± 1	16 ± 4	1.51 ± 0.27	⇑
(295) Theresia	30.50 ± 1.23	${0.220}_{-0.013}^{+0.013}$	${24}_{-17}^{+38}$	211 ± 3	45 ± 20	1.08 ± 0.19	⇑
(322) Phaeo	59.66 ± 1.29	${0.126}_{-0.009}^{+0.004}$	${12}_{-7}^{+11}$	210 ± 6	⋯	2.16 ± 0.39	⇓
(343) Ostara	19.65 ± 1.38	${0.113}_{-0.010}^{+0.009}$	${140}_{-60}^{+120}$	226 ± 7	⋯	1.51 ± 0.27	⇓
(444) Gyptis	162.3 ± 22.9	${0.049}_{-0.008}^{+0.007}$	${74}_{-74}^{+74}$	253 ± 25	⋯	1.08 ± 0.19	⇓

Notes.

^a Thermal inertia values are in SI units ( J m⁻² K⁻¹ s^−1/2). ^b a/b values are adjusted downward by 16% to account for the overestimation as described in MacLennan & Emery (2019). ^c Indicates either prograde ( ⇑ ) or retrograde ( ⇓ ) spin direction. ^d TPM results with ${\chi }_{\min }^{2}\gt$ 8 and thus should be used with caution.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We combine our sense of spin results with those in MacLennan & Emery (2019) and compare them to the spin poles of object shapes that are in the DAMIT ⁹ database (Database of Asteroid Models from Inversion Techniques; Ďurech 2010). In total, there are 101 objects both data sets and 77 of them have a sense of spin estimate that agrees. If we assume that the DAMIT spin axis estimates have a 100% accurate sense of spin, then the TPM has a 76.2% ± 4.3% success rate, based on binomial probability distribution. MacLennan & 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⁻² K⁻¹ 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 Figure 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 Section 4.

Zoom In Zoom Out Reset image size

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}_{\mathrm{eff}}^{\mathrm{WISE}}$) to the diameter values (one per object; ${D}_{\mathrm{eff}}^{\mathrm{TPM}}$) obtained here and plot them (colored by observing geometry) in Figure 4. There is a general agreement of within ±15% between the two data sets, with a few important notes. First, 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. Second, 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 Figure 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 (Section 3).

Zoom In Zoom Out Reset image size

A handful of asteroids with thermal inertias presented here have previous estimates from the works of Hanuš et al. (2018), Marciniak et al. (2019), and Pravec et al. (2019). We depict all these estimates in Figure 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 (Pravec et al. 2019)—our estimate is smaller by around a factor of 3, 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 388). If we were to only use the prograde solutions from our fitting, our estimate would not change, but if we consider higher roughness values that have a higher ${\tilde{\chi }}_{\min }^{2}$, 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.

Zoom In Zoom Out Reset image size

In addition to comparing the diameters and thermal inertias of individual objects from our data set to the findings from the WISE team, we compare the thermal inertia results of 220 asteroids from previous TPM works. Combined with the results from MacLennan & 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 five objects or more, notably Alí-Lagoa et al. (2020); Hanuš et al. (2015); Hanuš et al. (2016); Hanuš et al. (2018); and Marciniak et al. (2018, 2019). Alí-Lagoa et al. (2020) target some of the largest asteroids in the Main-Belt. Similar to this work and MacLennan & Emery (2019), Hanuš et al. (2015, 2018) have collectively modeled dozens of asteroids that were observed by WISE. Marciniak et al. (2018, 2019) specifically targeted asteroids with longer rotation periods (P_rot > 24 hr)—a group of objects that have lacked thermal inertia estimations. We refer to object thermal inertias presented in papers with fewer than five objects as "miscellaneous literature".

Table 4. TPM Results from the Literature

Object	HV	GV	Prot (hr)	Deff (km)	pV	Γa	Tc (K)	Source(s)
(1) Ceres	3.21	0.02	9.0742	951 ± 8	${0.100}_{-0.006}^{+0.004}$	${25}_{-10}^{+15}$	238 ± 15	[2]
(2) Pallas	4.08	0.08	7.8132 a	536 ± 5	${0.142}_{-0.005}^{+0.006}$	${30}_{-15}^{+15}$	234 ± 15	[2]
(3) Juno	5.29	0.34	7.2095	254 ± 4	${0.209}_{-0.019}^{+0.02}$	${70}_{-40}^{+30}$	245 ± 15	[2]
(4) Vesta	3.00	0.23	5.3421	530 ± 24	${0.394}_{-0.024}^{+0.011}$	${30}_{-10}^{+10}$	230 ± 10	[3]
(6) Hebe	5.71	0.27	7.2744 a	198 ± 3	${0.24}_{-0.01}^{+0.01}$	${50}_{-35}^{+40}$	205 ± 15	[19]
(8) Flora	6.35	0.24	12.865	142 ± 2	${0.252}_{-0.014}^{+0.014}$	${50}_{-30}^{+35}$	243 ± 15	[2]
(10) Hygeia	5.38	0.11	27.63	441 ± 6	${0.063}_{-0.002}^{+0.002}$	${50}_{-25}^{+20}$	210 ± 15	[2]
(16) Psyche	5.84	0.11	4.195 9 a	242.5 ± 25	${0.138}_{-0.015}^{+0.015}$	${120}_{-40}^{+40}$	212 ± 10	[21]

Note.

^a Indicates that the reported P_rot value has been rounded to four decimal places. References. [1] Alí-Lagoa et al. (2014); [2] Alí-Lagoa et al. (2020); [3] Capria et al. (2014); [4] Delbo' & Tanga (2009); [5] Devogèle et al. (2020); [6] Emery et al. (2014); [7] Hanuš et al. (2015); [8] Hanuš et al. (2016); [9] Hanuš et al. (2018); [10] Horner et al. (2012); [11] Jiang et al. (2020); [12] Keihm et al. (2012); [13] Leyrat et al. (2011); [14] Lim et al. (2011); [15] MacLennan & Emery (2019); [16] Marchis et al. (2012); [17] Marciniak et al. (2018); [18] Marciniak et al. (2019); [19] Marsset et al. (2017); [20] Matter et al. (2011); [21] Matter et al. (2013); [22] Mueller (2007); [23] Mueller et al. (2010); [24] Müller & Blommaert (2004); [25] Mueller & Lagerros (1998); [26] Müller et al. (2004); [27] Müller et al. (2011); [28] Müller et al. (2013); [29] Müller et al. (2014); [30] Müller et al. (2014); [31] Pravec et al. (2019); [32] Rozitis et al. (2013); [33] Rozitis & Green (2014); [34] Rozitis et al. (2014); [35] Rozitis et al. (2018); [36] Wolters et al. (2011); [37] Yu et al. (2014); [38] Yu et al. (2017).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

3.1. Multivariate Regression Model

We implement a forward stepwise multivariate regression model (Draper & 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; Section 2.4—object diameter (D_eff) and rotation period (P_rot). All variables, including Γ, are transformed into log₁₀ 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

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\rm{\Gamma }}=W+X({\mathrm{log}}_{10}{T}_{c})+Y({\mathrm{log}}_{10}{D}_{\mathrm{eff}})+Z({\mathrm{log}}_{10}{P}_{\mathrm{rot}}),\end{eqnarray} \tag{ 10 }$$

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 Figure 6. Previous studies quantified the thermal inertia dependence on diameter (e.g., Delbo' & Tanga 2009), rotation period (Harris & Drube 2016), and heliocentric distance (as a proxy for temperature; Rozitis et al. 2018) separately, 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.

Zoom In Zoom Out Reset image size

3.2. Noteworthy Objects

Our results show a slight excess of prograde spins at larger sizes (Figure 3), which is generally consistent with previous findings of spin vector distributions estimated from light-curve inversion methods (Kryszczynska et al. 2007; Hanuš et al. 2011; Ďurech 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 a prograde/retrograde difference in the >60 km size bin, which is consistent with the findings of Hanuš et al. (2011) and Ďurech 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 & 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 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 Equation (10) can be written as the proportionality ${\rm{\Gamma }}\propto {T}_{c}^{2.74\pm 0.29}$. Because surface temperatures generally scale with the inverse square of the heliocentric distance (Equation (7)), it follows that Γ ∝ r^α . Our results can be expressed in terms of this proportionality by using ${T}_{c}^{X}={R}_{\,\mathrm{au}}^{-2\alpha }$, 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 Figure 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 & Drube (2016); Rozitis et al. (2018) suggested that increased solar heating would allow the thermal wave to sample higher thermal inertia material in the subsurface 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 & 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². When combining the temperature dependence 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' & 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' & 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). Harris & Drube (2016), who used a NEATM-based thermal inertia estimator, found a significant correlation between thermal inertia and rotation period for asteroids with rotation periods spanning 2–200 hr. However, the works of Marciniak et al. (2018, 2019) found an abundance of low-Γ slow rotators (P_rot > 12 hr) 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 its relationship (or lack thereof) with the 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 (Rozitis et al. 2018), can be used to investigate possible changes in regolith properties as a function of depth(i.e., Harris & Drube 2016).

In this work, we applied the method of MacLennan & 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 Section 3. Our thermal inertia estimates are consistent with previous values from the literature for individual objects (Figure 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 the 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 (Section 3.1). The temperature dependence (${\rm{\Gamma }}\propto {T}_{c}^{2.74\pm 0.29}$) 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² if the temperature dependence of heat capacity is additionally considered. Instead, this relationship between thermal inertia and temperature is consistent with the temperature dependency of both the heat capacity and thermal conductivity (Section 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 data set. 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 multivariate regression model similar to that performed here in order to explore the controlling factors of regolith evolution on asteroids.

E.M.M. is supported by the NASA Earth and Space Science Fellowship #NNX14AP21H. We thank two anonymous reviewers for their thoughtful critiques, which improved the presentation of this paper. This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a joint project of the Jet Propulsion Laboratory/California Institute of Technology and the University of Arizona. NEOWISE is funded by the National Aeronautics and Space Administration.