Thermophysical Modeling of Asteroid Surfaces Using Ellipsoid Shape Models

Many versions of TPMs, with varying levels of complexity, have been developed throughout the past few decades. Such models explicitly account for the effects that physical attributes of the surface have on the thermal emission (Spencer et al. 1989). Subsurface heat conduction, self-shadowing from direct insolation (incoming solar radiation), multiply scattered insolation, and re-absorbed thermal radiation (i.e., self-heating) are implemented in a multitude of ways. The energy conservation expressed by Equation (1) is amended to include additional terms that account for these effects:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{S}_{\odot }(1-A)}{{R}_{\mathrm{au}}^{2}}\cos (i)(1-s)+{E}_{\mathrm{scat}}^{\mathrm{solar}}+{E}_{\mathrm{abs}}^{\mathrm{therm}}\\ \qquad +\,{\left.k\displaystyle \frac{{dT}}{{dx}}\right|}_{\mathrm{surf}}-{\varepsilon }_{B}\sigma {T}_{\mathrm{surf}}^{4}=0.\end{array}\end{eqnarray} \tag{ 2 }$$

Here, k is the effective thermal conductivity, i is the angle between a surface facet's local zenith and Sun-direction (solar incidence angle), and s is a binary factor indicating whether an element is shadowed by another facet (s = 1) or not (s = 0). The terms ${E}_{\mathrm{scat}}^{\mathrm{solar}}$ and ${E}_{\mathrm{abs}}^{\mathrm{therm}}$ encompass the contributed energy input from scattered solar radiation and re-radiated thermal photons, respectively, from other surface elements. Rough topography has been modeled as spherical section craters (Spencer 1990; Emery et al. 1998), random Gaussian surface (Rozitis & Green 2011), and using fractal geometry (Davidsson & Rickman 2014). The mathematics and numerical implementations of these features differ somewhat from one another. For further details, we refer the reader to the primary papers and to Davidsson et al. (2015) for a comparison of the different implementations.

TPMs often use the time-dependent, one-dimensional heat diffusion equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial T(x,t)}{\partial t}=\displaystyle \frac{k}{\rho c}\displaystyle \frac{{\partial }^{2}T(x,t)}{\partial {x}^{2}},\end{eqnarray} \tag{ 3 }$$

to model the flow of heat into and out from the subsurface. This presentation of Fourier's Law in Equation (3) assumes that the effective thermophysical factors (thermal conductivity, bulk density, ρ, and specific heat capacity, c) do not vary with depth or temperature. One-directional heat flow into the subsurface is also often assumed, but is well-justified by the concept that thermal energy flow is aligned with the temperature gradient—either directly up or down beneath the surface.

In the subsections below, we review and assess how various physical parameters (diameter, albedo, thermal inertia, surface roughness, shape, and spin direction) can be constrained in various observing circumstances. In particular, we highlight the optimal observing configurations and possible biases that may arise from particular viewing geometries.

One of the primary motivations for developing thermal models was to obtain size estimates of objects from disk-integrated thermal infrared observations (e.g., Allen 1970; Morrison 1973), given that thermal flux emission is directly proportional to the object's projected area (i.e., Equation (18)). A single value of diameter is not uniquely defined for a non-spherical body, so an effective diameter (D_eff) is given: the diameter of the sphere (Mueller 2007) having the same projected area as the object. During the span of data collection for a non-spherical body, the projected area will almost never be constant. In this case, D_eff is best reported as a time-averaged value, or adjusted using visible lightcurve data obtained simultaneously or proximal in time to the thermal observations (e.g., Lebofsky & Rieke 1979; Harris & Davies 1999; Delbo' et al. 2003; Lim et al. 2011). If a detailed shape model is being used, then the rotational phase of the object can be shifted in order to match the time-varying flux (Alí-Lagoa et al. 2013).

Geometric albedo (p_V) is calculated directly from D_eff and the object's absolute magnitude (H_V) (Russell 1916; Pravec & Harris 2007):

$$\begin{eqnarray}&&\sqrt{{p}_{V}}=\displaystyle \frac{1329\ \mathrm{km}\ {10}^{-0.2{H}_{V}}}{{D}_{\mathrm{eff}}[\mathrm{km}]}.\end{eqnarray} \tag{ 4 }$$

The geometric albedo is then converted to the bolometric Bond albedo using the definition of the phase integral, q:

$$\begin{eqnarray}&&q\equiv {A}_{V}/{p}_{V}=0.290+0.684\times {G}_{V}\end{eqnarray} \tag{ 5 }$$

in which G is the slope parameter (Bowell et al. 1989) and the approximation A ≈ A_V is made in Equation (1).

The characteristic ability of a material to resist temperature changes when subject to a change in energy balance is quantified by its thermal inertia (${\rm{\Gamma }}=\sqrt{k\rho c}$). Atmosphereless surfaces with low thermal inertia conduct little thermal energy into the subsurface, resulting in dayside temperatures that are close to instantaneous equilibrium with the insolation and extremely high diurnal temperature differences. Surfaces having large thermal inertia store and/or conduct thermal energy in the subsurface during the day, so significant energy is re-radiated during the nighttime hours, which results in a comparatively smaller diurnal temperature change. This influence that thermal inertia has on the surface temperature distribution manifests in the observed SED. Thus, photometry at two wavelengths (e.g., a measurement of the color temperature) can provide a better measure of thermal inertia, compared to that of a single wavelength (Mueller 2007). Figure 1(a) shows the normalized SEDs of surfaces (at opposition) having different thermal inertias, roughnesses, and shapes. Surfaces with higher roughness and low thermal inertia exhibit flux enhancement at small wavelengths near opposition, lowering the wavelength location of the peak flux. Data collected that have large wavelength separations spanning the blackbody peak are most useful because they are sensitive to relatively warmer and cooler portions of the surface—in other words, the overall temperature distribution. Delbo' & Tanga (2009) estimated thermal inertias for 10 main-belt asteroids using multi-wavelength IRAS (Infrared Astronomical Satellite) observations, and Hanus et al. (2015, 2018) collectively present thermal inertia estimates for 131 objects using Wide-Field infrared Survey Explorer (WISE) data. Harris & Drube (2016) present population trends using over 100 thermal inertias derived from the NEATM η values of asteroids observed with WISE.

Zoom In Zoom Out Reset image size

Multiple observations at a single wavelength but various solar phase angles—i.e., a thermal phase curve—have also proven to be useful methods to estimate thermal inertia (e.g., Spencer 1990). Figure 1(b) shows examples of thermal phase curves for three different objects, each possessing five different thermal inertia values. Müller et al. (2011, 2017) demonstrated that observations of (162173) 1993 JU₃, Ryugu, taken on one side of opposition but widely spaced in solar phase angle (30°), can constrain thermal inertia because they effectively have two points along a thermal phase curve. As seen in Figure 1(b), this approach can be optimized when the thermal phase curve spans across both sides of opposition. Thus, observing at pre- and post-opposition virtually guarantees that thermal emission information from both the warmer afternoon (α > 0) and cooler morning sides is gathered (Müller et al. 2014).

Surfaces with large degrees of roughness exhibit warmer dayside temperatures, due both to more of the surface area being pointed toward the Sun and to the effects of multiple scattering of reflected and emitted light. The term "beaming effect" is used to describe enhanced thermal flux return in the direction of the Sun at low phase angles. As pointed out by Rozitis (2017), the flux enhancement near opposition is highly sensitive to surface roughness. They also demonstrated that telescopic observations at a near pole-on illumination and viewing geometry can be used to effectively constrain the degree of roughness, specifically when α < 40° and the sub-solar latitude is greater than 60° (Figures 1(b) and (c)). Figure 1(b) shows thermal phase curves of objects with smooth and rough surfaces of varying thermal inertia. Figure 1(c) shows the thermal flux emitted as a function of sub-solar latitude and recreates the findings of Rozitis (2017). In general, the surface roughness of airless bodies are more easily estimated near opposition—or better yet, with disk-resolved observations that offer a greater range of viewing geometries. The Rozitis (2017) study shows that the uncertainty in thermal inertia is slightly larger for an object observed at a large phase angle, as a consequence of the difficulty in estimating surface roughness.

Disk-integrated flux is also directly affected by the shape and spin vector of an object (D˘urech et al. 2017). If thermal observations sample the entire rotation period, a thermal light curve is useful for constraining these parameters. In particular, the amplitude of the thermal lightcurve is largely affected by the elongation and orientation of the spin vector. In principle, constraints can be placed on the spin vector—or, at least, spin direction—of objects because the temperature distribution is affected by the sub-solar latitude and thus the spin vector (Müller et al. 2011, 2012, 2014, 2017). Figure 1(c) shows the variation of the flux for an ellipsoid at different viewing geometries. Gathering multiple thermal lightcurves can offer amplitude measurements at different viewing geometries, which significantly improves determination of the global shape. Both Morrison (1977) and Hansen (1977b) correctly determine the spin direction by observing the change in diameter estimates for data acquired before and after opposition. The diameter estimate of (1) Ceres, a prograde rotator, before opposition was 5% larger than the diameter estimates from post-opposition observations. Müller et al. (2014) explicitly mentions that observations at either side of opposition are useful indicators of the spin direction. As shown in Figure 1(a), this effect is due to the flux excess emitted from the hotter afternoon side of a prograde rotator as it is would be observed before opposition (e.g., Lagerros 1996).

The text and figure above demonstrate how multi-wavelength thermal lightcurve observations that span across an object's blackbody peak region, at both pre- and post-opposition (with Δα > 40°), contain information directly dependent on an object's physical properties. Thermal inertia, surface roughness, and spin direction can all be constrained from measurements of an object's thermal phase curve, particularly when observations widely vary in solar phase angle. Additionally, thermal lightcurve amplitude measurements at more than one viewing geometry are a unique indicator of the elongation of a rotating object. The work described in the rest of this paper demonstrates and quantifies the effectiveness of using an ellipsoidal TPM to model the thermophysical properties of any given asteroid, given the data set described above.

In development of the smooth TPM, we find it particularly useful to parameterize the depth variable in Equation (3) as x' = x/l_s (Spencer et al. 1989), where l_s is the thermal skin depth, the length scale at which the amplitude of the diurnal temperature variation changes by a factor of e ≈ 2.718:

$$\begin{eqnarray}&&{l}_{s}=\sqrt{\displaystyle \frac{k}{\rho c}\displaystyle \frac{{P}_{\mathrm{rot}}}{2\pi }},\end{eqnarray} \tag{ 6 }$$

where P_rot is the rotation period of the body. This parameterization transforms the temperature gradient term in Equation (2) as follows:

$$\begin{eqnarray}&&{\left.k\displaystyle \frac{{dT}}{{dx}^{\prime} }\right|}_{\mathrm{surf}}\Rightarrow {\rm{\Gamma }}\sqrt{\displaystyle \frac{2\pi }{{P}_{\mathrm{rot}}}}{\left.\displaystyle \frac{{dT}}{{dx}^{\prime} }\right|}_{\mathrm{surf}}.\end{eqnarray} \tag{ 7 }$$

In order to further reduce the number of independent input variables, we also parameterize the temperature and time as T' = T/T_eq and t' = 2πt/P_rotΘ. The thermal parameter, Θ³ is given by:

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

Spencer et al. (1989) introduced this dimensionless parameter, which accounts for factors that affect the diurnal temperature variation. This parameter was realized by comparing the diurnal rotation of a body to the timescale in which thermal energy is stored and then re-radiated, per unit surface area (Spencer et al. 1989). Ignoring the effects of multiple-scattering and self-heating, this parameterization scheme changes Equation (2) to

$$\begin{eqnarray}&&{\left.\cos (i)+{\rm{\Theta }}\displaystyle \frac{{dT}^{\prime} }{{dx}^{\prime} }\right|}_{\mathrm{surf}}-{T}_{\mathrm{surf}}^{{\prime} 4}=0,\end{eqnarray} \tag{ 9 }$$

and Equation (3) to

$$\begin{eqnarray}&&\displaystyle \frac{\partial T^{\prime} (x^{\prime} ,t^{\prime} )}{\partial t^{\prime} }={\rm{\Theta }}\displaystyle \frac{{\partial }^{2}T^{\prime} (x^{\prime} ,t^{\prime} )}{\partial {x}^{{\prime} 2}}.\end{eqnarray} \tag{ 10 }$$

Surface temperatures are calculated across the surface of a sphere by solving Equation (10), given the upper boundary condition Equation (9). Because the amplitude of diurnal temperature changes decreases exponentially, the heat flux approaches zero, and the lower boundary becomes:

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{dT}^{\prime} }{{dx}^{\prime} }\right|}_{x^{\prime} \to \infty }=0.\end{eqnarray} \tag{ 11 }$$

Using the parameterized depth, time, and temperature, the number of input variables in this model is effectively reduced to two (the thermal parameter and sub-solar latitude). A finite-difference approach is used to numerically implement Equation (10), as detailed in Appendix A. This spherical, smooth-surface TPM consisted of 13 latitude bins and 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 logarithmic space, from 0 to 450) in order to generate the surface temperature look-up tables of T'.

Our rough, surface-cratered TPM is similar to that originally presented by Hansen (1977a), which was improved upon by both Spencer (1990) and Emery et al. (1998) to include heat conduction, multiple scattering of insolation, and re-absorption of thermal radiation within spherical craters. Following the procedure of Emery et al. (1998), craters are constructed with m = 40 planar elements contained within k = 4 rings that are radially symmetric about the crater center. The kth ring outward contains 4k elements, all of which are forced to have the same surface area (see Figures 1 and 2 in Emery et al. 1998, for crater depiction). As was mentioned by Spencer (1990), craters with more than four rings significantly increase the overall computational time and do not enhance the model resolution over the four-ringed crater in most cases. The overall geometry of these craters is characterized by the half-opening angle, γ, as measured from the center-line of the crater to the edge (e.g., a hemispherical crater has γ = 90°). The overall degree of surface roughness is characterized by the mean surface slope ($\bar{\theta };$ Hapke 1984), which is only a function of γ and the fraction of area covered by craters, f_R (Lagerros 1996):

$$\begin{eqnarray}&&\tan \bar{\theta }=\displaystyle \frac{2{f}_{R}}{\pi }\displaystyle \frac{\sin (\gamma )-\mathrm{ln}[1+\sin (\gamma )]+\mathrm{lncos}(\gamma )}{\cos (\gamma )-1}.\end{eqnarray} \tag{ 12 }$$

As shown in Emery et al. (1998), the fraction of energy transferred to one crater element from another is (conveniently) equal among all elements. The scattered solar and re-absorbed thermal radiation received by the ith crater facet from the jth facet are

$$\begin{eqnarray}&&{E}_{i,\mathrm{scat}}^{\mathrm{solar}}=\displaystyle \frac{{S}_{\odot }(1-A)}{{R}_{\mathrm{au}}^{2}}\displaystyle \frac{A}{1-A\tfrac{\gamma }{\pi }}\displaystyle \frac{1-\cos (\gamma )}{2m}\sum _{i\ne j}^{m}\cos ({i}_{j})\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{E}_{i,\mathrm{abs}}^{\mathrm{therm}}=\displaystyle \frac{1-\cos (\gamma )}{2m}(1-{A}_{\mathrm{th}})\sum _{i\ne j}^{m}{\varepsilon }_{B}\sigma {T}_{j}^{4},\end{eqnarray} \tag{ 14 }$$

respectively. The energy balance at the surface (Equation (2)), in our parameterized time, temperature, and depth environment, becomes:

$$\begin{eqnarray}&&\begin{array}{l}\cos (i)(1-s)+\displaystyle \frac{1-\cos (\gamma )}{2m}\\ \qquad \times \,\left(\displaystyle \frac{A}{1-A\tfrac{\gamma }{\pi }}\sum _{i\ne j}^{m}\cos ({i}_{j})+(1-{A}_{\mathrm{th}})\sum _{i\ne j}^{m}{T}_{j}^{{\prime} 4}\right)\\ \qquad +\,{\rm{\Theta }}{\left.\displaystyle \frac{{dT}^{\prime} }{{dx}^{\prime} }\right|}_{\mathrm{surf}}-{T}_{\mathrm{surf}}^{{\prime} 4}=0.\end{array}\end{eqnarray} \tag{ 15 }$$

Because planetary surfaces are highly absorbing at infrared wavelengths, the Bond albedo at thermal-infrared wavelengths (A_th) is assumed to be zero. Thus, only singly scattered re-absorption of thermal emission within the crater is considered, in contrast to multiple scattering and reabsorption of the insolation. Temperature lookup tables are generated for craters of the same values of sub-solar latitude and thermal parameter as the smooth surface TPM runs. This was done for three sets of craters of varying opening angle, in order to simulate changes in roughness. Unlike the parameterized version of the smooth surface energy balance equation, A is included as an independent input parameter, requiring us to explicitly account for changes in this parameter. Fortunately, it only appears in the solar scattering term, which is not a major contributor to the total energy budget. Our spherical, rough-surface TPM consists of 13 different latitude bins and was run for 3 values of γ = {45°, 68°, 90°}, 46 values of sub-solar latitude (0° to 90° in 2° increments), 116 values of the thermal parameter (spread out in logarithmic space, from 0 to 450), and 7 values of A^grid = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 1}⁴ to construct this large set of T' lookup tables.

Shape model facet temperatures are independent of one another in our TPM because global self-heating does not occur on the spherical and ellipsoidal shapes considered here. Because the insolation upon a facet is only dependent on the orientation of the surface normal to the Sun, it is possible to map, or transform, surface temperatures from one shape to another. We exploit this fact in order to map the surface temperatures from a spherical body to convex DAMIT shape models and ellipsoids. We make use of two kinds of coordinate systems, which are depicted in Figure 2. The body-centric (θ, ϕ) coordinate system defines the latitude and longitude of a facet relative to the center of the shape model. Shape model facets are often described using sets of vectors, r, given in body-centric coordinates. Alternatively, (ϑ, φ) describes a coordinate system that describes the tilt of a facet relative to the local surface, which we call the surface-normal coordinates. Both ϑ and φ are calculated using the surface normal vector n shown in Figure 2 (for a sphere, the body-centric and facet-normal coordinates are equivalent). For any given facet, on any shape model, with ϑ and φ, the temperatures can be equated to a facet on a sphere that has the same surface-normal coordinates. In Appendix B, we derive closed-form analytic expressions to map temperatures from a sphere to the effective coordinates of a triaxial (a ≥ b ≥ c) ellipsoid.

Zoom In Zoom Out Reset image size

Thermal flux is calculated by a summation of the individual flux contributions from smooth surface and crater elements visible to the observer (i.e., Equation (18)). To calculate the flux of an object not having exact Θ or sub-solar latitude values included in our lookup temperature tables, we calculate the fluxes for the closest grid points and perform a linear interpolation to compute the fluxes for the desired parameters. Instead of doing the same for A, we calculate the flux for a value from A^grid and then multiply by an adjustment factor, Λ. This approach saves time and computational cost by not computing an interpolation across three variables (or dimensions). A similar approach was described by Wolters et al. (2011), in which they calculated fluxes for a perfectly absorbing surface and employed a correction factor based on the desired A. Here, Λ was tested empirically and adjusts the model flux based on the blackbody T_eq curves of the desired A and the closest value of A^grid.

$$\begin{eqnarray}&&{\rm{\Lambda }}(A,{A}^{\mathrm{grid}})=\displaystyle \frac{1}{2}\left[1+{\left(\displaystyle \frac{B(\lambda ,{T}_{\mathrm{eq}}({A}^{\mathrm{grid}}))}{B(\lambda ,{T}_{\mathrm{eq}}(A))}\right)}^{2/3}\right].\end{eqnarray} \tag{ 16 }$$

The χ² goodness-of-fit statistic is evaluated using the modeled flux points, F_m(λ), and the observed flux measurements, F_o(λ), with the associated 1σ uncertainty, σ_o, and in general:

$$\begin{eqnarray}&&{\chi }^{2}=\sum \displaystyle \frac{{({F}_{m}(\lambda )-{F}_{o}(\lambda ))}^{2}}{{\sigma }_{o}^{2}}.\end{eqnarray} \tag{ 17 }$$

We note here that, instead of individual flux measurements, the modeled and observed fluxes used in Equation (17) represent only the extracted mean and peak-to-trough range of thermal lightcurve fluxes, along with their uncertainties as calculated via error propagation. As part of the data-fitting routine, the free parameters of the model are varied in order to minimize the goodness-of-fit as described in the following subsection.

A weighted sum of fluxes to simulate a surface that is comprised of smooth and rough surface patches:

$$\begin{eqnarray}\begin{array}{rcl}F(\lambda ) & = & \displaystyle \frac{{\varepsilon }_{\lambda }}{{{\rm{\Delta }}}_{\mathrm{au}}^{2}}{\int }_{S}\int \{(1-{f}_{R}){B}_{\mathrm{smooth}}(\lambda ,T(\theta ,\phi ))\\ & & +\,{f}_{R}(1-v){\rm{\Lambda }}{B}_{\mathrm{rough}}(\lambda ,T(\theta ,\phi ,i))\}\cos (e){dA}.\end{array}\end{eqnarray} \tag{ 18 }$$

Here, Δ_au is the observer-centric distance in astronomical units. The visibility factor, v, is analogous to the shadowing factor; v = 1 if a facet is hidden from view, otherwise v = 0. Each point on the surface is treated as a blackbody emitter with wavelength-dependent emissivity, ε_λ, and temperature, T = T_eqT':

$$\begin{eqnarray}&&B(\lambda ,T)=\displaystyle \frac{2{{hc}}^{2}}{{\lambda }^{5}}\displaystyle \frac{1}{\exp ({hc}/\lambda {k}_{b}T)-1}.\end{eqnarray} \tag{ 19 }$$

At each epoch, both the mean flux and peak-to-trough range of the thermal light curve are extracted and used in the model fitting procedure. As alluded to in Section 2, these two parameters contain diagnostic information about the object's shape, spin direction, and thermophysical properties; they can be easily extracted from non-dense thermal lightcurve data. In principle, it is just as feasible to fit models to each independent flux point that contributes to an object's infrared lightcurve. Doing so could offer insight into the shape, as departures from a sinusoidal lightcurve can indicate relative topographic lows or highs. However, such an approach works best in cases in which the thermal lightcurve is densely sampled, which is often not the case in untargeted astronomical surveys such as IRAS, Akari, and WISE. In this work, we focus on estimating only the elongation of a body by incorporating the photometric range (peak-to-trough) of the thermal lightcurve. TPM fitting to sparse thermal lightcurves is prone to systematic parameter bias from oversampling and/or heteroscedastic uncertainties across many rotational phases, which can unevenly emphasize certain rotational phases over others, thus skewing the best-fit parameters. An example of this would be a scenario in which peaks of the modeled lightcurve were fit better with the observed lightcurve minima, thus resulting in an overall lower best-fit modeled fluxes. However, the capacity and capability of incorporating thermal lightcurve data into the shape model inversion process, when it is combined with visible lightcurve data, should be explored further (D˘urech et al. 2014, 2015).

In our data-fitting approach, the shape, spin vector (${\lambda }_{l}^{\mathrm{eclip}}$, ${\beta }_{l}^{\mathrm{eclip}}$), roughness, and thermal inertia are left as free parameters that we select from a predefined sample space, and we search for the best-fit D_eff. A sphere and ellipsoids with b/c = 1 (i.e., prolate) with a/b axis ratios of 1.25, 1.75, 2.5, and 3.5 are used. For each of these shapes, we sample through 25 predefined thermal inertia values, 3 default roughness ($\bar{\theta }$) values, and 235 spin vectors. Following Table 1 in Delbo' & Tanga (2009), each value of γ is paired with a corresponding f_R = {0.5, 0.8, 1.0} value in order to produce default mean surface slopes of $\bar{\theta }=\{10^\circ ,29^\circ ,\ \mathrm{and}\ 58^\circ \}$. The thermal inertia points are evenly spread in log space from 0 to 3000 Jm⁻² K⁻¹ s^−1/2, and the spin vectors are spread evenly throughout the celestial sphere. For each shape/spin vector/Γ combination, we use a routine to find the D_eff value which minimizes χ².

The grid of spin vectors are formed by constructing a Fibonacci lattice in spherical coordinates (Swinbank & Purser 2006). Here, a Fermat spiral is traced along the surface of the celestial sphere, using the golden ratio (Φ ≈ 1.618) to determine the turn angle between consecutive points along the spiral. The result is a set of points with near-perfect homogeneous areal coverage across the celestial sphere. This Fibonacci lattice constructed here uses N = 235 points, with the lth point ($i\in [0,N-1]$) having an ecliptic latitude and longitude of

$$\begin{eqnarray}&&{\beta }_{l}^{\mathrm{eclip}}=\arcsin \left(\displaystyle \frac{2l-N+1}{N-1}\right)\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{\lambda }_{l}^{\mathrm{eclip}}=2\pi l{{\rm{\Phi }}}^{-1}\mathrm{mod}2\pi ,\end{eqnarray} \tag{ 21 }$$

respectively. The mean flux value and the peak-to-trough range at each wavelength and at each epoch are taken as input to Equation (17). The Van Wijngaarden–Dekker–Brent minimization algorithm, commonly known as Brent's Method (Brent 1973), combined with the golden section search routine (Section 10.2 Press et al. 2007), is used because it does not require any derivatives of the χ² function to be known. The algorithm uses three values of χ² evaluated at different input diameters to uniquely define a parabola. The algorithm used the minima of these parabolas to iteratively converge on the D_eff/A combination that corresponds to the global minimum of χ² value to within 0.01%.

The values of D_eff and A are linked through Equations (4) and (5) for a smooth surface, but the effect of multiple scattering alters the energy balance, and thus effective albedo (Mueller 2007),⁵ for a cratered surface:

$$\begin{eqnarray}&&{A}^{\mathrm{crater}}(\gamma )=A\displaystyle \frac{1-{\sin }^{2}(\gamma /2)}{1-A{\sin }^{2}(\gamma /2)}.\end{eqnarray} \tag{ 22 }$$

An object having an areal mixture of both smooth and rough topography has an effective bond albedo, A_eff, which is a weighted average of the albedo of a smooth surface and the bond albedo of a crater A^crater (Wolters et al. 2011):

$$\begin{eqnarray}&&{A}_{\mathrm{eff}}({f}_{R},\gamma )=(1-{f}_{R})A+{f}_{R}{A}^{\mathrm{crater}}.\end{eqnarray} \tag{ 23 }$$

To place confidence limits on each of the fitted parameters, we use the χ² values calculated during the fitting procedure. In general, the χ² distribution depends on ν degrees of freedom (three in the present work), which is equal to the number of data points (constraints; eight in this work, as described in Section 4.1) minus the number of input (free; five in this work) parameters. Because the χ² distribution has an expectation value of ν and a standard deviation of $\sqrt{2\nu }$, the best-fit solutions cluster around ${\chi }_{\min }^{2}=\nu$, and those with ${\chi }^{2}\lt (\nu \,+\sqrt{2\nu })$ represent 1σ confidence estimates.⁶ However, because of the TPM assumptions (e.g., no global self-heating, homogeneous thermophysical properties), it is often possible for the TPM to not perfectly agree with the data, in which case ${\chi }_{\min }^{2}\gt \nu$. In this case, we use ${\chi }^{2}/{\chi }_{\min }^{2}\lt \nu +\sqrt{2\nu }$ to place 1σ confidence limits. This modification effectively scales the cutoff bounds by a factor of ${\chi }_{\min }^{2}$, instead of adopting the traditional approach of using a constant χ² distance cutoff. This scaling of the χ² cutoff bounds is a more conservative approach in quoting parameter uncertainties, as it includes the systematic uncertainties that lead to the larger ${\chi }_{\min }^{2}$ in the reported parameter uncertainties. Using the reduced χ² statistic (${\tilde{\chi }}^{2}={\chi }^{2}/\nu$), we can express the solutions within a 1σ range as ${\tilde{\chi }}^{2}\lt {\tilde{\chi }}_{\min }^{2}(1+\sqrt{2\nu }/\nu )$.

In Section 5, we report TPM results for all of the objects that we analyzed, even those with a high ${\chi }_{\min }^{2}$ (indicating a poor fit to the data), which are considered unreliable and should only be used with caution.

Table 1 lists the objects, DAMIT shape models used, the associated rotation periods, and the effective diameter, as computed beforehand via a NEATM fit to the real WISE observations. The observing geometries used in the synthetic model runs are the same as the WISE observations listed in Table 3. We calculate the artificial thermal emission for WISE photometric filters W3 and W4 (12 and 24 μm) for a full rotation of the shape. From these infrared lightcurves, we extract the mean and peak-to-trough flux range for each filter as input into our fitting routine (Equation (17)). Because we now have these two quantities, observed twice for two wavelengths, there are eight data points to constrain five free parameters: diameter, thermal inertia, surface roughness, shape elongation, and sense of spin. Fluxes for each shape model are calculated using the TPM described here by using the actual WISE observing circumstances detailed in Table 3. If an object has two shape models available—a result of ambiguous solutions of the lightcurve inversion algorithm—then both were included in the analysis.

Table 1.  Shape Models and Physical Properties for Synthetic Data Set

Object	Shape Model	a/bsynth	Spin Vector	Prot (hr)	${D}_{\mathrm{eff}}^{\mathrm{synth}}$ (km)	HV	GV
	M171	1.31	−69°, 107°
(167) Urda				13.06133	43.00	9.131	0.283
	M172	1.30	−68°, 249°
(183) Istria	M669	1.39	20°, 8°	11.76897	35.48	9.481	0.221
	M182	1.22	−75°, 20°
(208) Lacrimosa				14.0769	44.35	9.076	0.232
	M183	1.23	−68°, 176°
(413) Edburga	M354	1.37	−45°, 202°	15.77149	35.69	9.925	0.296
	M521	1.41	24°, 90°
(509) Iolanda				12.29088	60.07	8.476	0.382
	M522	1.32	54°, 248°
(771) Libera	M250	1.50	−78°, 64°	5.890423	29.52	10.28	0.323
	M609	2.30	34°, 38°
(857) Glasenappia				8.20756	15.63	11.29	0.246
	M610	2.32	48°, 227°
(984) Gretia	M256	1.57	52°, 245°	5.778025	35.76	9.526	0.379
(1036) Ganymede	M261	1.05	−78°, 190°	10.313	37.42	9.236	0.311
	M403	1.76	−73°, 12°
(1140) Crimea				9.78693	31.77	9.621	0.207
	M404	1.61	−22°, 175°
(1188) Gothlandia	M479	1.71	−84°, 334°	3.491820	13.29	11.52	0.254
	M409	2.03	35°, 106°
(1291) Phryne				5.584137	27.87	10.29	0.251
	M410	2.30	59°, 277°
	M657	1.27	54°, 225°
(1432) Ethiopia				9.84425	7.510	12.02	0.282
	M658	1.32	44°, 41°
(1495) Helsinki	M656	1.74	−39°, 355°	5.33131	13.54	11.41	0.359
(1568) Aisleen	M422	2.30	−68°, 109°	6.67597	13.60	11.49	0.131
	M477	1.56	70°, 222°
(1607) Mavis				6.14775	15.10	11.32	0.256
	M478	1.93	59°, 0°
(1980) Tezcatlipoca	M274	1.72	−69°, 324°	7.25226	5.333	13.57	0.186
(2156) Kate	M438	2.15	74°, 49°	5.622153	8.678	4.339	0.186
	M704	1.68	−50°, 197°
(4611) Vulkaneifel				3.756356	12.38	11.87	0.268
	M705	1.64	−86°, 5°
	M682	3.48	−78°, 97°
(5625) 1991 AO2				6.67412	7 14.25	12.83	0.165
	M683	3.18	−52°, 265°
	M757	1.55	67°, 62°
(6159) 1991 YH				10.65893	5.148	13.38	0.175
	M758	1.65	67°, 266°

Download table as:  ASCIITypeset image

The best-fit TPM parameter results for the synthetic data set are detailed below and depicted in Figure 3. Overall, using ellipsoid shapes results in more accurate and precise estimates for diameter, thermal inertia, and shape.

Zoom In Zoom Out Reset image size

Diameter Constraints. Figure 3 includes diameter uncertainties, shown by gray bars, based off the noise of WISE data. The uncertainty overlap with the expected diameter values indicates that the uncertainty in the data is equal to or larger than the uncertainty introduced by model assumptions. The assumption of a spherical shape results in overestimation of diameter of up to 20%, as shown in Figure 3(a). The TPM performs better when using ellipsoid shapes; at most, diameters are overestimated by 10%. These offsets are particularly pronounced for highly elongated objects with high Θ values and when observed at large sub-solar latitude values. In these cases, the average observed cross-sectional area is particularly large and surface temperatures are, on average, warmer (Figure 1(c)), because a larger fraction of the surface experiences perpetual daylight. Diameter uncertainties when using ellipsoid shapes are consistent with the ±10% value seen in other thermophysical modeling papers, yet seem to imply a drop in accuracy from the NEATM (Harris 2006). However, care must be taken when comparing the performance results presented here, which are based on synthetic data generated from shape models with various spin vectors, and those presented by Harris (2006), which uses an idealized synthetic data set generated from spherical shapes with no obliquity. The work of Wright (2007), who included rough spheres in their performance test of NEATM, shows diameter accuracy of the NEATM to be on par with TPM works, such as this one.

Thermal Inertia/Parameter Constraints. We transform thermal inertia estimates into thermal parameter space (Equation (8)) because the rotation periods of synthetic objects affect the temperature distribution and must be accounted for. Panels (c) and (d) in Figure 3 show the estimated Θ values for the spherical and best-fit ellipsoidal shape, respectively, against the synthetic values. When a sphere is used for a TPM, thermal inertia is overestimated at low values and underestimated at large Θ values. At either extreme of Θ, this offset is systematically different by a factor of ∼4. Ellipsoidal shapes do a much better job matching the input Θ, with systematic offsets of less than 25% (red dashed line in Figure 3(d)). We investigate this discrepancy further in Section 4.3, but note here that the estimated uncertainties (gray bars), calculated from WISE signal-to-noise ratios for each object, are far larger than the systematic offset, indicating that the assumptions in our TPM fitting do not contribute much to the overall uncertainty in thermal inertia and the χ²-based errors values reasonably reflect the overall precision in our TPM.

Shape Constraints. To assess how well our method is able to constrain the shape (elongation) of an object, we compare the best-fit ellipsoid a/b to that of the area-equivalent a/b ratio of the DAMIT shape model. Because we use a sparse grid of possible a/b values, placing meaningful confidence bars based on χ² values is not practical. Instead, we plot the frequency distribution of %a/b (=$\tfrac{a/{b}_{\mathrm{ellip}}}{a/{b}_{\mathrm{DAMIT}}}-1)$ difference in Figure 3. Fitting a Gaussian function (16%, with standard deviation of 18%) to the distribution shows that the assumptions within the model, most likely the use of a prolate ellipsoid (a/c = 1), result in a slight overestimation of the a/b axis. The standard deviation of this distribution is not reflective of the uncertainty in the mean offset in the shape accuracy, but rather of the inherent model uncertainties that arise from assuming an ellipsoid shape. Accounting for this discrepancy, we shift the TPM best-fit shape result downward by 16% and assign an uncertainty of 18%.

Spin Constraints. We also investigated the ability to constrain the spin direction of an object (i.e., retrograde or prograde rotation). In total, the ellipsoid TPM correctly identified the spin direction 58.3 ± 6.5% of the time, compared to 67.1 ± 5.7% for the spherical TPM. In Section 2.4, we pointed out that thermal inertia is a significant factor when constraining spin direction. Thus, we break down our results into the thermal inertia bins assigned for the artificial data set. Figure 3(f) shows the percentage of best-fit solutions, broken down by sphere/ellipsoid shape assumption, in which spin direction was correctly identified for a given thermal inertia. It is clear that intermediate values of thermal inertia provide more reliable constraints on the spin direction, which can be explained by the increased asymmetry of thermal phase curves at intermediate values of thermal inertia. More extreme (low and high) values of thermal inertia result in more symmetric thermal phase curves, which make it difficult to distinguish morning and afternoon hemispheres—and thus spin direction (Figure 1(a)).

We performed a linear multiple regression analysis on the percent diameter and Θ difference (%ΔD_eff and %ΔΘ_eff, between the synthetic and TPM estimates) to distinguish which factors, if any, bias the estimates. The predictor variables chosen were log₁₀(Θ), the sub-solar latitude (s-s lat.) or sub-observer latitude (s-o lat.), and the elongation of the shape model (a/b^synth), because they are the most likely to affect the response variable (Section 2). The regression models the response variable—the percent difference in diameter or Θ—as being linearly dependent on the sum of any number of independent predictor variables—in this case, four. For each predictor, the regression model computes a slope coefficient that represents the change in that variable when all others are held constant. An intercept term, quantifying the value of the response when all predictors are zero, is also computed. Table 2 shows the slope coefficients, intercepts, and the associated p-values that indicate the statistical significance of each predictor variable in the multiple regression model. We use p < 0.05, representing a 95% confidence level, to identify predictor variables that affect the response variable (p-values that are larger than this cutoff indicate that the predictor variable has little-to-no statistically relevant effect on the response variable).

Table 2.  Multiple Linear Regression Results

	%ΔDeff		%ΔΘ
	coefficient	p-value	coefficient	p-value
${D}_{\mathrm{eff}}^{\mathrm{synth}}$	−0.0009 ± 0.0004	0.02	0.0012 ± 0.0017	0.48
log10(Θsynth)	−0.0037 ± 0.0067	0.59	−0.1390 ± 0.0292	<0.01
s-s lat.synth	⋯	⋯	−0.0036 ± 0.0024	0.12
s-o lat.synth	0.0024 ± 0.0005	<0.01	⋯	⋯
a/bsynth	0.0085 ± 0.0119	0.48	0.0573 ± 0.0507	0.26
intercept	−0.0264 ± 0.0319	0.41	−0.0621 ± 0.1357	0.65

Download table as:  ASCIITypeset image

Table 3.  WISE Observation Circumstances and Fluxes

Object	UT Datea	Δtobsb	Nobsc	Raud	Δaue	α (°)f	$\overline{W3}$ g	$\diamond W3$ h	$\overline{W4}$ g	$\diamond W4$ h
(167) Urda	2010 Feb 9	1.257	13	2.840	2.647	20.33	451.8 ± 5.2	141.6 ± 10.3	1224 ± 26	350.2 ± 65.8
	2010 Aug 2	1.257	12	2.787	2.510	−21.27	601.1 ± 7.2	185.9 ± 13.8	1559 ± 30	480.7 ± 51.1
(183) Istria	2010 Feb 4	3.772	22	3.718	3.580	15.38	87.65 ± 1.27	36.88 ± 2.57	322.7 ± 9.0	119.9 ± 16.1
	2010 Jul 21	1.257	13	3.765	3.562	−15.62	70.84 ± 1.19	43.98 ± 2.23	274.7 ± 8.0	160.0 ± 14.9
(208) Lacrimosa	2010 Feb 10	0.595	9	2.888	2.700	19.98	508.2 ± 5.7	119.8 ± 10.6	1358 ± 28	344.8 ± 57.3
	2010 Aug 8	1.257	10	2.913	2.645	−20.30	442.4 ± 4.9	97.78 ± 9.89	1269 ± 25	255.1 ± 48.3
(413) Edburga	2010 Feb 10	1.257	14	3.180	3.009	18.08	164.9 ± 2.2	72.85 ± 3.98	534.7 ± 11.5	208.9 ± 24.1
	2010 Jul 26	1.257	15	2.701	2.427	−22.01	478.4 ± 5.8	180.3 ± 11.6	1219 ± 25	410.4 ± 42.0
(509) Iolanda	2010 Jan 18	0.596	8	3.342	3.164	17.11	440.3 ± 5.1	230.6 ± 9.9	1405 ± 30	594.3 ± 67.6
	2010 Jul 3	1.257	15	3.327	3.080	−17.72	470.5 ± 6.0	216.4 ± 12.2	1461 ± 24	649.9 ± 40.7
(771) Libera	2010 Jan 29	1.257	10	2.764	2.593	20.88	223.7 ± 2.8	74.50 ± 5.24	618.4 ± 15.4	223.4 ± 28.7
	2010 Jul 16	3.903	22	3.111	2.847	−18.97	146.2 ± 2.1	75.75 ± 4.31	474.6 ± 11.2	223.2 ± 19.7
(857) Glasenappia	2010 Jan 19	1.125	14	2.329	2.085	24.99	180.8 ± 2.4	70.40 ± 5.26	373.4 ± 9.3	127.4 ± 23.6
	2010 Jul 11	1.389	17	2.172	1.823	−27.75	279.2 ± 3.5	80.68 ± 6.9	527.0 ± 12.7	164.6 ± 21.7
(984) Gretia	2010 Jan 31	1.125	9	3.324	3.188	17.24	161.3 ± 2.2	76.53 ± 4.46	531.1 ± 10.99	214.3 ± 19.85
	2010 Jul 21	1.125	10	3.159	2.920	−18.71	192.6 ± 2.5	81.18 ± 5.08	604.3 ± 14.8	258.7 ± 32.8
(1036) Ganymede	2010 Jan 15	0.860	9	3.897	3.729	14.61	72.50 ± 1.10	12.50 ± 2.16	294.7 ± 7.0	39.25 ± 14.02
	2010 Jun 22	1.125	14	3.463	3.241	−17.02	141.3 ± 2.1	27.94 ± 4.14	494.8 ± 12.1	81.66 ± 24.09
(1140) Crimea	2010 Feb 13	1.125	11	3.075	2.900	18.73	184.4 ± 2.4	77.33 ± 4.85	524.7 ± 12.0	201.6 ± 26.8
	2010 Aug 1	1.124	11	2.990	2.725	−19.76	220.7 ± 2.9	110.0 ± 5.1	624.5 ± 12.6	290.7 ± 24.4
(1188) Gothlandia	2010 Jan 18	0.992	10	2.580	2.354	22.41	59.19 ± 1.11	41.00 ± 2.05	158.7 ± 6.3	99.88 ± 12.78
	2010 Jul 2	1.389	16	2.521	2.222	−23.68	88.37 ± 1.40	58.34 ± 2.47	214.4 ± 6.8	129.9 ± 11.8
(1291) Phryne	2010 Jan 21	0.992	8	3.210	3.038	17.85	117.5 ± 1.8	122.0 ± 3.6	367.8 ± 9.7	350.8 ± 18.2
	2010 Jul 8	1.257	14	3.090	2.825	−19.12	131.0 ± 2.1	142.4 ± 4.3	393.8 ± 11.4	386.3 ± 18.9
(1432) Ethiopia	2010 Feb 4	3.771	23	2.699	2.511	21.43	15.44 ± 0.60	4.102 ± 1.217	43.91 ± 3.48	17.13 ± 8.23
	2010 Jul 28	1.389	15	2.307	1.992	−26.02	32.16 ± 0.76	7.282 ± 1.496	78.34 ± 3.84	14.90 ± 7.49
(1495) Helsinki	2010 Feb 8	4.301	19	2.490	2.284	23.33	90.81 ± 1.46	42.36 ± 3.0	208.3 ± 6.6	81.53 ± 12.85
	2010 Aug 5	1.125	10	2.267	1.940	−26.47	140.6 ± 2.0	87.26 ± 4.38	299.8 ± 8.3	154.2 ± 16.5
(1568) Aisleen	2010 Jan 24	1.125	12	2.950	2.776	19.50	27.39 ± 0.68	14.91 ± 1.37	82.59 ± 3.85	51.08 ± 7.46
	2010 Jul 7	1.257	14	3.100	2.550	−21.04	52.16 ± 0.90	21.74 ± 1.84	135.1 ± 4.8	63.58 ± 9.57
(1607) Mavis	2010 Jan 21	0.993	10	3.285	3.116	17.43	28.51 ± 0.71	20.45 ± 1.44	94.41 ± 4.53	63.07 ± 8.75
	2010 Jul 4	1.125	12	3.039	2.778	−19.47	35.16 ± 0.73	14.46 ± 1.44	107.9 ± 3.9	48.81 ± 8.31
(1980) Tezcatlipoca	2010 Jan 23	0.992	11	2.333	2.109	24.96	13.79 ± 0.53	5.748 ± 1.041	34.52 ± 2.83	22.68 ± 5.61
	2010 Jun 30	1.257	8	2.064	1.716	−29.39	29.65 ± 0.73	18.61 ± 1.54	65.74 ± 3.46	42.50 ± 6.61
(2156) Kate	2010 Jan 25	0.992	9	2.688	2.502	21.49	23.82 ± 0.68	16.42 ± 1.38	64.43 ± 3.85	34.34 ± 8.27
	2010 Jul 11	1.126	11	2.625	2.322	−22.67	24.93 ± 0.66	15.41 ± 1.29	66.16 ± 3.66	37.33 ± 7.50
(4611) Vulkaneifel	2010 Jan 25	0.993	10	3.103	2.944	18.50	21.82 ± 0.66	10.80 ± 1.36	72.11 ± 3.44	31.03 ± 7.05
	2010 Jul 11	1.257	13	3.100	2.834	−19.05	25.08 ± 0.65	13.63 ± 1.25	81.09 ± 3.86	40.84 ± 7.26
(5625) 1991 AO2	2010 Jan 28	0.993	11	3.180	3.032	18.04	28.80 ± 0.72	19.98 ± 1.43	91.72 ± 4.18	49.67 ± 8.86
	2010 Jul 14	1.125	13	3.168	2.900	−18.61	34.55 ± 0.83	31.91 ± 1.63	107.8 ± 4.3	100.6 ± 7.7
(6159) 1991 YH	2010 Feb 3	4.301	19	2.341	2.124	24.91	13.21 ± 0.84	11.16 ± 1.49	31.40 ± 3.64	27.42 ± 7.01
	2010 Jul 30	1.389	11	2.424	2.120	−24.67	14.90 ± 0.59	11.15 ± 1.17	36.87 ± 3.36	27.63 ± 6.68

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

^aUT date of the first observation. ^bTime spanned by observations (days). ^cNumber of observations used. ^dMean Heliocentric distance. ^eMean WISE-centric distance. ^fMean solar phase angle. ^gLightcurve-averaged mean flux. ^hPhotometric range of lightcurve.

Download table as:  ASCIITypeset image

The diameters calculated from our approach are overestimated when observations occur at high sub-observer latitudes (p < 0.01) and for small diameters (p = 0.02) as marked in gray in Table 2. From the multiple regression results performed on Θ, the only predictor that explains the variance in response variable is Θ^synth, marked in Table 2 with gray. The negative sign represents an overestimation for small Θ and/or underestimation at large Θ. A re-examination of panel (d) in Figure 3 shows that this offset is due to points with Θ > 6. For lower values of Θ, the best-fit line is likely skewed upward due to the underestimation at the larger end. Thus, we employ a formulation to correct objects for which Θ > 6:

$$\begin{eqnarray}&&{{\rm{\Theta }}}^{\mathrm{corr}}={10}^{\left(\tfrac{b}{1+m}\right)}{{\rm{\Theta }}}^{\left(\tfrac{1}{1+m}\right)}.\end{eqnarray} \tag{ 24 }$$

To determine the correction factor, we fit a line to objects with Θ > 4 (shown by the purple-dotted line in Figure 4) and compute m = −0.136 and b = 0.026. This fit provides a means of removing the systematic underestimation in thermal parameter, and by proxy, thermal inertia. The bias seen in panel (d) of Figure 3 vanishes, as shown by the red best-line in Figure 4 re-computed after using Equation (24) to adjust the Θ^synth > 6 values. Thermal parameters this large are more likely to be found for objects with unusually high thermal inertia and/or among slow rotators—or for icy bodies with low surface temperatures in the outer solar system (e.g., Table 1 in Spencer et al. 1989).

Zoom In Zoom Out Reset image size

Upon closer inspection of the multiple regression results, we see that the skew at small sizes is due to particular objects being observed at high latitudes rather than by the intrinsic diameter of the object. The effect that the viewing geometry has on the size estimation can be explained by the fact that the thermal lightcurve-averaged flux does not properly represent the effective diameter. For the extreme case of observations taken from a "pole-on" geometry, the thermal flux represents the largest feasible cross-sectional area, which would result in the overestimation of diameters, regardless of the shape used in the TPM.

In some of the largest publications of independently derived thermal inertias (Delbo' & Tanga 2009; Hanus et al. 2015, 2018), convex shape models were used in the TPM approach. One critique of using convex shape models is that actual asteroids can potentially harbor shape concavities, which raises potential concern as to whether or not large deviations from a spherical shape will bias the temperature distribution in any significant way (Rozitis & Green 2013). Radar observations of the NEA (341843) 2008 EV₅ showed a large concavitiy, and thermal observations were analyzed in detail by Alí-Lagoa et al. (2013). They found no evidence that the concavity influenced the results, as its effects on the thermal emission (i.e., shadowing and global self-heating) were likely below the signal-to-noise of the observations. Using a larger set of shape models with concavities, Rozitis & Green (2013) demonstrated that the effects of global self-heating and shadowing have negligible effects on the temperature distribution, compared to those of thermal inertia and surface roughness; this is consistent with the findings of Lagerros (1997), who investigated the effect using objects random Gaussian shapes. All of these studies justify the use of convex, prolate ellipsoid studies in our TPM approach.

Although not the case here, thermal inertia bias could arise if the size is fixed during the fitting procedure. For example, diameters that are estimated from radar observations can suffer from an overestimated z-axis if the object is observed at near-equatorial view, yielding a higher thermal inertia (Rozitis & Green 2014). Some studies have investigated how changes in shape and spin vector can affect the value of thermal inertia. Using both a spherical and radar shape model for the asteroid (101955) Bennu, Emery et al. (2014) found that the radar shape model gave a lower thermal inertia result compared to a spherical shape with the same spin vector. The lower thermal inertia is explained by the oblateness of Bennu, as surface facets are systematically tilted away from the Sun-direction, cooling the surface temperatures relative to that of a sphere and requiring a lower model thermal inertia to compensate. As briefly mentioned before, our use of prolate ellipsoids may likely be the cause of overestimation of the elongation of objects. For example: if we were to make b/c = 1.1, the average orientation of the shape facets would veer away from the Sun and effectively lower the modeled surface temperatures. By lengthening the b/c axis, there is less need to lengthen the a/b axis in order to match the surface temperature.

It is becoming increasingly common for thermal infrared observations to be used to refine the detailed shape models created from the delay-Doppler radar and/or the visible lightcurve inversion techniques (D˘urech et al. 2017). Shape models from these methodologies can sometimes yield incorrect z-axes if the object is observed nearly equator-on (Rozitis et al. 2013; Rozitis & Green 2014). In short, thermal data constrains the effective diameter—and by extension, the z-axis dimension—if the x- and y-axis dimensions were known to great accuracy from the observations. Additionally, Hanus et al. (2015) utilized thermal infrared data to refine convex shape models derived from lightcurve inversion. They varied a convex shape model within the uncertainty of the photometric errors to generate a set of shapes that were used to generate thermal fluxes from a TPM, which were then fit to WISE observations. Although not much work has been done to derive new shapes from only disk-integrated thermal emission data, approaches that combine thermal observations with optical lightcurves (D˘urech et al. 2012, 2014) and other telescopic observations (e.g., stellar occultations, optical interferometry, and delay-Doppler radar; see D˘urech et al. (2015, 2017)) have been employed in order to refine pre-existing shape models. Generally, the refined shape models appear smoother (D˘urech et al. 2012; Hanus et al. 2015), which is likely in part due to the thermal emission being sensitive to large-scale curvature, especially at lower thermal inertia values (Lagerros 1996).

The WISE mission, an astrophysics mission designed to map the entire sky, operated in its fully cryogenic mode from 2010 January 14 to August 5 at wavelengths centered near 3.4, 4.6, 12, and 22 μm, denoted W1, W2, W3, and W4 (Wright et al. 2010). A data-processing enhancement called NEOWISE (Mainzer et al. 2011a) detected moving solar system objects, most of which were asteroids in the main belt and in near-Earth orbits. During each grouping of observations (epoch), a moving object is typically detected around 10 to 20 times, in ∼1.6 hr multiples—the orbital period of the spacecraft. This means that flux measurements are separated in time by more than ∼1.6 hr. Therefore, depending on the object's range of motion on the sky, each epoch of observations can potentially span up to 36 hr. NEOWISE reports each moving object detection to the Minor Planet Center (MPC⁷ ), where the start time, R.A., and decl. of each observation can be retrieved. The set of times and locations are used to parse the WISE All-Sky Single Exposure (L1b) catalog on the Infrared Science Archive (IRSA) maintained by the Infrared Science and Analysis Center (IPAC⁸ ). We select detections reported to within 10'' and 10 s of those reported to the MPC. These constraints are "relaxed" relative to the accuracy of the telescope's astrometric precision, to guarantee that IPAC returns flux information for each reported MPC observation. These criteria also return many spurious detections. In the following, we discuss our method of rejecting spurious sources returned from these generous search criteria.

According to the WISE Explanatory Supplement (Cutri et al. 2012), photometric profile fits are unreliable for W3 < −3 mag or W4 < −4 mag, due to saturation of the detector. Because nonlinearity is present for sources with W3 < 3.6 mag and W4 < −0.6 mag, we increase the magnitude uncertainty to 0.2 mag for objects brighter than these values (Mainzer et al. 2011b). We shift the isophotal wavelengths and zero magnitude point of W3 and W4 to account for the red-blue calibrator discrepancy described in the Explanatory Supplement. The raw magnitudes that are reported were calibrated assuming that the flux across each filter was that of Vega's spectrum. Thus, a color-correction must be made to account for the discrepancy between the spectrum of the object and that of Vega (Wright et al. 2010; Cutri et al. 2012). For each individual observation, NEATM was used to calculate the flux spectrum across the full bandpass for each WISE band for use in the color correction, which ranged from 0.87 to 1.0 for W3 and was nearly constant at 0.98 for W4. Based on an analysis of asteroid flux uncertainties in consecutive frames by Hanus et al. (2015), we increase the flux error in W3 and W4 by factors of 1.4 and 1.3, respectively.

In order to filter out bad observations of an asteroid—for example, in a situation where it passes near a background star or when the query returns a detection of an unwanted object—we employ "Peirce's criterion"⁹ (Peirce 1852; Gould 1855) as outlined and demonstrated by Ross (2003). We flag spurious observations using this algorithm based on W4–W3, because it will identify and remove sources of an anomalous color temperature. Using color, rather than raw flux, avoids the possibility of removing seemingly anomalous observations of the minimum or maximum flux of a highly elongated object. This rejection method is a simpler alternative to that implemented by Alí-Lagoa et al. (2013) and Hanus et al. (2015, 2018) on WISE data, in which the WISE inertial source catalog was checked explicitly for possible flux contamination of stars.

Mean fluxes, denoted as $\overline{W3}$ and $\overline{W4}$, were calculated by simply taking the error-weighted mean of all observations. The photometric range for each band, denoted as $\diamond W3$ and $\diamond W4$, was calculated by subtracting the minimum and maximum fluxes. Standard error propagation was used to estimate 1σ uncertainties for each of these parameters, and then input as σ_o in Equation (17). These bright, slow-moving objects with P_rot < 36 hr have sufficient coverage in the rotational phase to provide estimates of these parameters. However, some objects may have only been detected a handful of times, which could only sparsely sample the rotational phases. In a follow-up paper, we will develop and present a more rigorous method for accurately estimating the flux mean and range, with associated errors, from sparse lightcurve coverage.

A summary of our TPM fits and corresponding 1σ uncertainties is given in Table 4. For diameter, albedo, and thermal inertia, uncertainties are based on the χ² values calculated during the fitting routine. In some cases, two values of roughness could not be distinguished with respect to one being better than the other, so both values are reported. In the case of (6159) Andreseloy, all roughness values tried provided statistically indistinguishable fits to the data. The spin direction reported for each object here reflects the preferred spin vector. Often, many spin vector solutions lie within the 1σ uncertainty bounds, so reporting a single solution would not be meaningful. The spin direction of (1036) Ganymede could not be independently determined here, which is likely due to the nearly spherical shape of the object combined with its low thermal inertia reducing the asymmetry of the thermal phase curve, as demonstrated in Figure 1(b).

Table 4.  WISE Data TPM Results

Object	Deff (km)	pV	Γa	$\bar{\theta }(^\circ )$	a/bb	Spinc	${\tilde{\chi }}_{\min }^{2}$	MBA/NEA
(167) Urda	39.48 ± 0.89	${0.252}_{-0.019}^{+0.010}$	${51}_{-16}^{+20}$	38 ± 13	1.51 ± 0.27	⇓	1.27	MBA
(183) Istria	31.43 ± 2.92	${0.288}_{-0.033}^{+0.029}$	${21}_{-10}^{+12}$	47 ± 13	1.51 ± 0.27	⇑	5.16	MBA
(208) Lacrimosa	40.44 ± 1.37	${0.253}_{-0.014}^{+0.012}$	${77}_{-22}^{+31}$	41 ± 13	1.51 ± 0.27	⇑	1.66	MBA
(413) Edburga	33.44 ± 1.75	${0.169}_{-0.022}^{+0.012}$	${41}_{-10}^{+19}$	9 ± 6	1.51 ± 0.27	⇓	1.98	MBA
(509) Iolanda	54.39 ± 3.86	${0.243}_{-0.027}^{+0.019}$	${8.6}_{-8.6}^{+12.2}$	18 ± 7	1.51 ± 0.27	⇑	3.01	MBA
(771) Liberad	29.23 ± 2.10	${0.160}_{-0.019}^{+0.013}$	${61}_{-26}^{+34}$	53 ± 39	1.51 ± 0.27	⇓	13.5	MBA
(857) Glasenappia	13.62 ± 0.84	${0.297}_{-0.021}^{+0.013}$	${58}_{-24}^{+30}$	<20	2.16 ± 0.39	⇑	1.52	MBA
(984) Gretia	34.72 ± 1.18	${0.227}_{-0.023}^{+0.012}$	${28}_{-7}^{+8}$	53 ± 10	1.51 ± 0.27	⇑	1.81	MBA
(1036) Ganymede	35.85 ± 1.95	${0.278}_{-0.027}^{+0.017}$	${15}_{-15}^{+22}$	40 ± 32	1.08 ± 0.19	⋯	0.55	NEA
(1140) Crimea	30.13 ± 1.18	${0.276}_{-0.028}^{+0.020}$	${23}_{-23}^{+15}$	23 ± 21	1.51 ± 0.27	⇓	1.50	MBA
(1188) Gothlandia	13.52 ± 0.84	${0.238}_{-0.022}^{+0.019}$	${38}_{-13}^{+21}$	41 ± 27	1.51 ± 0.27	⇓	0.43	MBA
(1291) Phryne	27.03 ± 1.65	${0.186}_{-0.017}^{+0.014}$	${20}_{-6}^{+16}$	45 ± 17	2.16 ± 0.39	⇑	1.78	MBA
(1432) Ethiopia	7.15 ± 0.67	${0.535}_{-0.070}^{+0.058}$	${71}_{-65}^{+180}$	⋯	1.51 ± 0.27	⇑	0.56	MBA
(1495) Helsinki	13.31 ± 0.59	${0.271}_{-0.033}^{+0.017}$	${19}_{-13}^{+13}$	12 ± 8	1.51 ± 0.27	⇑	0.52	MBA
(1568) Aisleen	11.66 ± 1.01	${0.328}_{-0.038}^{+0.034}$	${51}_{-22}^{+41}$	46 ± 38	2.16 ± 0.39	⇓	2.74	MBA
(1607) Mavisd	14.52 ± 1.72	${0.249}_{-0.037}^{+0.032}$	${37}_{-25}^{+42}$	58 ± 50	1.51 ± 0.27	⇑	21.45	MBA
(1980) Tezcatlipoca	5.68 ± 0.58	${0.205}_{-0.040}^{+0.035}$	${170}_{-110}^{+170}$	57 ± 39	1.51 ± 0.27	⇓	1.67	NEA
(2156) Kate	8.04 ± 0.45	${0.294}_{-0.025}^{+0.021}$	${56}_{-23}^{+23}$	49 ± 32	2.16 ± 0.39	⇑	0.39	MBA
(4611) Vulkaneifel	12.10 ± 1.12	${0.216}_{-0.028}^{+0.023}$	${32}_{-32}^{+23}$	<53	1.51 ± 0.27	⇓	1.01	MBA
(5625) Jamesferguson	14.46 ± 0.86	${0.062}_{-0.006}^{+0.005}$	${52}_{-15}^{+14}$	>36	2.16 ± 0.39	⇓	2.46	MBA
(6159) Andreseloy	5.65 ± 1.37	${0.247}_{-0.061}^{+0.061}$	${60}_{-60}^{+177}$	⋯	1.51 ± 0.27	⇓	0.19	MBA

Notes.

^aThermal inertia values are in SI units (Jm⁻² K⁻¹ s^−1/2). ^ba/b values are adjusted downward by 16% to account for the overestimation as described in Section 4.2. ^cIndicates either prograde (⇑) or retrograde (⇓) spin direction. ^dTPM results with ${\chi }_{\min }^{2}\gt 8$ and thus should be used with caution.

Download table as:  ASCIITypeset image

Comparing the diameter estimates from NEATM fits presented by the WISE team (i.e., Mainzer et al. 2011b; Masiero et al. 2011) to the values obtained here, we observe agreement (within ±15% of another) between the two sets of results. However, for objects ∼40 km and above, TPM diameter estimates are systematically higher. The objects in question (Urda, Lacrimosa, and Iolanda) did not saturate, nor were they bright enough to lie in the nonlinear regime of the WISE detectors. This discrepancy may be due to one the flux corrections described above or from the model differences between our TPM approach and the NEATM used by the WISE team.

When comparing to other previous works, our results are consistent with the thermal inertias reported: (771) Libera and (1980) Tezcatlipoca have thermal inertia estimates of 65+85/−35 and 220+380/−204 Jm⁻² K⁻¹ s^−1/2, respectively, made by Hanus et al. (2015), though the reader should take note that the high chi-squares indicate that our fits for Libera were relatively poor, contributing to the relatively high parameter uncertainties. (1036) Ganymede has several thermal inertia estimates: 24 ± 8 (Rozitis et al. 2018), 35+65/−29 (Hanus et al. 2015), and 214 ± 80 Jm⁻² K⁻¹ s^−1/2 (Rivkin et al. 2017).¹⁰ The higher estimates are around an order of magnitude greater than the smallest estimates and can be explained by the thermal inertia dependency on surface temperature. Our thermal inertia estimate for Ganymede was based on data collected at R_au = 3.5 and 3.8, and is thus consistent with similar estimates at the same distance.

Eleven objects in this study have been analyzed by the recent work of Hanus et al. (2018), in which the varied-shape TPM was used to derive thermophysical properties from the WISE data set. A comparison of the thermal inertias for each of these objects in shown in the right panel of Figure 5. Some objects have two estimates presented in Hanus et al. (2018) due to ambiguous shape models that produce equivalent fits to the data. For 7 out of these 11 objects, there is very good agreement (i.e., the estimates are within the 1σ uncertainties) between the estimates from their work and ours. We note that our model fit for (1607) Mavis was inaccurate, as indicated by the large ${\tilde{\chi }}_{\min }^{2}$ in Table 4, for which the error bars for each parameter are noticeably large¹¹ and consistent with Hanus et al. (2018). For two objects, (413) Edburga and (984) Gretia, the 1σ error bars just barely miss overlapping, and two others, (167) Urda and (1495) Helsinki, have very different estimates. The two estimates Hanus et al. (2018) present for (167) Urda are just over twice as large as ours and have very small reported uncertainties of ±5. Hanus et al. (2018) report a thermal inertia just over zero for (1495) Helsinki, also with a very small uncertainty. For each of these objects, our 2σ uncertainty bounds encompass the Hanus et al. (2018) estimates and both works overall appear to deliver no systematically different estimates from another.

Zoom In Zoom Out Reset image size

With the combination of our TPM results with the subset compiled by Delbo' et al. (2015) and the large data set of Hanus et al. (2018), we note that the number of known thermal inertias of objects in the 5–50 km size range increases by 20. Restricting our analysis to the compiled thermal inertias from Delbo' et al. (2015), we detect a negative correlation with diameter (Figure 6) evidenced by the Spearman's rank coefficient of r_s = −0.55; this correlation remains when including the ∼120 thermal inertias of Hanus et al. (2018). From the findings of Rozitis et al. (2018), however, thermal inertia should be adjusted to account for its dependency on heliocentric distance by using the formula: ${\rm{\Gamma }}={{\rm{\Gamma }}}_{0}{R}_{\mathrm{au}}^{a}$. Doing this allows for a better comparison of hot, small NEAs and cool, larger MBAs because the effects of temperature are partially accounted for. We use a = 3/4 here, which was suggested previously by Delbo' et al. (2015) and Mueller et al. (2010). Even when performing the adjustment on asteroid thermal inertias, the correlation remains statistically significant (with r_s = −0.41 and p < 0.001) as shown in the bottom panel of Figure 6. However, this correlation becomes insignificant if we use a = 4/3, which is large but well within the range of empirically derived values from analyzing individual objects (i.e., Ganymede and 2002 CE₂₆; Rozitis et al. 2018). While the purpose of this work is not to investigate this dependency, we plan a follow-up work to increase the number thermal inertia estimates of small main-belt asteroids, combine our results with Hanus et al. (2018), and revisit this dependency in much greater detail.

Zoom In Zoom Out Reset image size

The prograde/retrograde spin directions reported here are in agreement with those from DAMIT (Table 1), with the exception of three objects: (208) Lacrimosa, (1495) Helsinki, and (6159) Andreseloy. Incorporating the ambiguous spin direction for Ganymede, this means that 17 out of 21 (81 ± 9%) objects matched the spin vectors in the DAMIT shape models. Comparing this result to the findings of our proof-of-concept study in Section 4.2, there is a noticeable difference. Assuming the DAMIT shape models are 100% accurate, there is much improvement in the percent of synthetic objects in which the spin direction was correctly identified (∼58 and ∼67% for ellipsoids and a sphere, respectively). When broken down by thermal inertia, as in panel (f) of Figure 3, these results are comparable to the best-case scenario (spherical shape for Γ = 40 Jm⁻² K⁻¹ s^−1/2) and outperform each situation in which an ellipsoid shape is used. A reasonable explanation for this discrepancy is the inclusion of roughness in the real-world application of the TPM, which adds to the asymmetry of the thermal phase curve—particularly for asteroids in the main belt (e.g., Figure 1(b)).