An Investigation of the Ranges of Validity of Asteroid Thermal Models for Near-Earth Asteroid Observations

We simulate thermal flux densities using a TPM. TPMs allow for simulating the asteroid's surface temperature distribution by solving the heat transfer equation numerically for a large number of plane surface facets. The TPM we use is mostly identical to, and has been tested extensively against, the one discussed by Mueller (2007) and has been used in other works, including Mommert et al. (2014a, 2014b). Our TPM implementation accounts for the spin axis orientation (longitude and latitude in ecliptic coordinates), rotational period, thermal inertia, and surface roughness in addition to observing geometry. TPMs can provide more accurate diameter and albedo estimates than simple thermal models, but they also require additional constraints on the aforementioned parameters. In this work, we use the TPM to simulate flux densities for synthetic asteroids.

The TPM solves the heat transfer equation numerically for each surface facet until the sub-surface temperature at a certain depth varies less than a certain fractional threshold, ε_temp (usually ε_temp = 10⁻⁴), between two full rotations. For the sake of simplicity, we use a spherical shape for our model asteroids, each consisting of 496 triangular surface facets. Surface roughness is taken into account and modeled as emission from spherical craters using the approach taken by Mueller (2007); different models (no roughness, low roughness, default roughness, and high roughness) as defined in Müller et al. (2004) are utilized here. In order to verify the conservation of energy in the model, the amount of incoming and outgoing energy are compared against each other. If there is an energy discrepancy of 1% or greater, ε_temp is reduced by an order of magnitude. This adjustment improves the conservation of energy and leads to more accurate flux densities—but also requires more computation time.

The wavelengths at which the flux densities are simulated sample the peaks of the asteroids' spectral energy distributions well: 5–20 μm at intervals of 5 μm. Note that reflected solar light is not considered as part of this work. Flux density uncertainties, which are necessary for the proper interpretation of the flux densities as part of the χ² minimization in the thermal models (see Section 2.3), are assumed to follow a Poisson distribution and hence are the square-root of the derived flux at each wavelength (in order to mimic the noise properties of real observations).

We create a synthetic representation of the NEA population consisting of 1 million objects sampling the parameters of our TPM (see Section 2.1). The ranges we consider for the physical properties are typical for NEAs, and we weight them according to their measured or modeled distributions. Orbital parameters are chosen from the population of known NEAs, as detailed below.

In order to simulate flux densities, we have to assume diameters and geometric albedos for our synthetic NEAs. Geometric albedos (p_V) are drawn from the albedo distribution by Wright et al. (2016) measured for the NEA population during the cryogenic part of the WISE mission. Because the surface temperature distribution of an asteroid is independent of its size, we assume a constant absolute magnitude H_V = 18 mag (G = 0.15) for all objects; the target diameter follows from H_V and the respective geometric albedo value.

We draw ecliptic spin axis latitudes from a second-order polynomial distribution fitted to the obliquity distribution derived by Farnocchia et al. (2013, their Figure 6); this distribution is in agreement with later findings (Tardioli et al. 2017). Spin axis longitudes are drawn from a uniform distribution spanning 360° as there is no preference for spin axis longitudes. Other parameters are also drawn from uniform distributions due to a lack of additional information on the underlying distributions. We sample thermal inertias in the range ${\rm{\Gamma }}\in [20,2000]$ SI units (1 SI Unit = 1 J m⁻² s^−0.5 K⁻¹), which is motivated by Mommert et al. (2014a) and the fact that little is known about the unbiased distribution of NEA thermal inertias. Rotational periods are uniformly sampled from the interval $P\in [10\,{\rm{s}},100\,\mathrm{hr}];$ this range is consistent with the range of reliably measured rotational periods of NEAs compiled by Warner et al. (2009). We note that the distribution of rotational periods of NEAs and other asteroids is strongly biased (e.g., Masiero et al. 2009). Hence, we deliberately use a uniform distribution of rotational periods in order to probe the impact of rotational periods over a wide range. Finally, we uniformly sample the discrete roughness models (low, default, and high roughness; see Mueller 2007, for definitions of these models); we consider the absence of surface roughness unrealistic for this study. We refer the reader to Section 4 for a discussion of the range definitions and their implications.

In order to properly sample observation geometry, we use heliocentric and geocentric coordinates of 5000 randomly selected real NEAs for a random epoch. All selected NEAs have solar elongation angles 60° ≤ ε ≤ 180°, restricting the sample to objects that are observable from the Earth and from space-based platforms with state-of-the-art thermal-infrared telescopes. Hence, our set of orbital geometry parameters provides a snapshot of the NEA population that is accessible for observations. We verified that this subset of the NEA population represents the overall population well in heliocentric and geocentric distances, solar phase angle, and as a function of time (using different epochs). This verification was performed using a two-sided Kolmogorov-Smirnov test that compares the full ensemble of observable NEAs at different epochs with our random sample; in no case was there a significant probability that our sample was drawn from a different population.

For each of the 5000 NEA geometries, we create 200 randomized sets of physical properties, resulting in a total of 1 million synthetic NEAs considered in this study. We do not consider noise contributions to the sampled parameters in order to more effectively work out the effects at work.

The simulated flux densities are modeled using well-tested versions of the NEATM and the FRM, following the implementation by Delbó & Harris (2002). For each synthetic NEA, we use all four generated thermal flux densities and fit them simultaneously taking into account the observing geometry and flux density uncertainties. The best fit is derived using a χ² minimization of the residuals between "measured" (generated by the TPM) and thermal model flux densities, weighted by the "measured" flux density uncertainties. The use of realistic (Poissonian) "measured" flux density uncertainties is necessary to weight the "measured" flux densities at the different wavelengths correctly. Improper uncertainties will lead to distorted best-fit spectral energy distributions and hence systematic effects on the derived diameters and albedos. Using this approach, we fit diameters and beaming parameters (NEATM only) for each object—albedos follow through the combination with the assumed absolute magnitude of the object. We use the same solar system absolute magnitude as a measure for optical brightness for all objects, which does not affect the conclusion of this work. We refer to Mommert (2013) for a detailed discussion of the model implementations.

For each synthetic object, we can identify whether the NEATM or the FRM provides the more accurate diameter or albedo result by comparing our solutions to the input parameters utilized in the derivation of the flux densities using the TPM (see Section 2.1). For instance, the diameter derived with the NEATM is considered more accurate if $| {d}_{\mathrm{NEATM}}-d| \lt | {d}_{\mathrm{FRM}}-d|$, where d is the input diameter, d_NEATM is the diameter derived by the NEATM, and d_FRM is the diameter derived by the FRM. In order to characterize the impact of one specific physical or orbital parameter on diameter and albedo estimates, we marginalize over the other parameters and derive distributions where the NEATM (or the FRM) is more accurate. The success rate is defined as the fraction of cases where the NEATM (or FRM) provides a better fit over a marginalized distribution. The marginalized distribution is derived from the complete synthetic object distribution by summation over the success rates for all but one parameter. Hence, the success rate distribution for this one parameter is fully retained, whereas information from the other parameters has been marginalized out. For instance, the NEATM diameter success rate with respect to thermal inertia describes the fraction of cases in which the NEATM derives a more accurate diameter than the FRM as a function of thermal inertia; in this case, the success rates for all parameters except that of the thermal inertia are summed over, retaining the success rate distribution with respect to this parameter. This fraction can be interpreted as the probability of the NEATM to be superior to the FRM as a function of thermal inertia for a random asteroid. Success rates for the diameter and the albedo should be identical; we provide both quantities for the sake of completeness and discuss exceptions from this rule.

In order to quantify the agreement of thermal model results with the actual synthetic asteroid diameters and albedos, we introduce a second diagnostic quantity. For each synthetic asteroid, we calculate the quantity (d_NEATM − d)/d, which is the fractional error of the NEATM-derived diameter (d_NEATM) from the actual input diameter (d); errors for the FRM and for the geometric albedo can be defined correspondingly. By marginalizing over one parameter, we can define the median fractional error as a function of this parameter. The marginalization is done as described above. However, in this case, the median is derived instead of the sum of the other parameters. The marginalization enables the identification of systematic effects that occur as a function of individual parameters. The variation in the fractional errors is captured by the 1σ and 3σ envelopes that include 68.3% and 99.7% of all synthetic objects around the median in each bin, respectively.

Figure 2 shows success rates and fractional errors for rotational period, thermal inertia, surface roughness, and geometric albedo. NEATM success rates are greater than 80%, and the NEATM median fractional errors are zero within the 1σ intervals for all parameters, whereas the FRM fractional errors are systematically and significantly offset. Conditions that resemble the assumptions of the FRM (high thermal inertia, short rotation periods) lead to slightly higher success rates and lower fractional errors for the FRM and vice versa for the NEATM. However, this effect is not significant. Increasing degrees of surface roughness lead to increasing fractional errors for both the NEATM and the FRM, but with different signs. This is expected, as both the NEATM and the FRM assume Lambertian surfaces that are perfectly smooth. Finally, diameters and geometric albedos derived with both models are unaffected by the input geometric albedos of the synthetic objects.

Zoom In Zoom Out Reset image size

The physical properties have only minor effects on the applicability of the NEATM, which is superior to the FRM in the vast majority of cases.

The impact of the spin axis orientation on thermal modeling is shown in Figure 3. The NEATM success rates are greater than 90%, and NEATM fractional errors are negligible throughout the entire range in both parameters. While the fractional errors for the spin axis latitude show a dip/bump, the fractional errors for the spin axis longitude show a wavy pattern. The former effect, a dip in the fractional errors of the diameters (albedos) derived with the NEATM (FRM) and a bump in the fractional errors of the diameters (albedos) derived with the FRM (NEATM) at spin axis latitudes close to zero, is a result of the different surface temperature distributions assumed by the two thermal models that are more or less compatible with equator-on views. Similarly, the NEATM has higher success rates for spin axis latitudes around zero and the FRM has higher success rates for extreme spin axis latitudes (−90° and +90°). At spin axis latitudes close to zero (and a random spin axis longitude), it is somewhat likely to have insolation and hence higher temperatures at the poles of the object, mimicking the NEATM surface temperature distribution and preferring the NEATM. In this case, the FRM overestimates diameters and underestimates albedos. At spin axis latitudes close to the extremes (−90° and +90°), the largest fraction of the insolation falls near the equator, leading to a surface temperature distribution very similar to that of the FRM, penalizing the NEATM performance. The wavy pattern in the case of the spin axis longitude is caused by our selection effects for the orbital geometries: all geometries used in this simulation have the same epoch and solar elongations between 60° and 180° (see Section 2.2). For those synthetic objects with spin axis latitudes close to zero, the spin axis longitude decides whether insolation is focused on the poles or the equator. Hence, NEATM-preferred and FRM-preferred surface temperature distributions are spaced by 90° in spin axis longitude, causing a wave-like pattern in the fractional errors of this parameter.

Zoom In Zoom Out Reset image size

Despite the complexity caused by the different spin axis configurations, the NEATM is more likely to provide accurate diameters and albedos for a random spin axis orientation of a random NEA.

We investigate the impact of the synthetic objects' observing geometry on the thermal modeling results. As heliocentric distance, distance from the observer, and solar phase angle are closely related, we focus on the latter, which shows the effects most clearly (Figure 4). At solar phase angles up to 65°, the NEATM has a success rate close to unity and its fractional errors in diameter and albedo are close to zero. For phase angles greater than 65°, the FRM diameter success rate is greater than that of the NEATM. The situation is more complicated for the albedo success rates of the NEATM and the FRM. This behavior is caused by the fact that the absolute diameter fractional errors of the NEATM grow faster than those of the FRM; FRM diameters are closer to the real diameters for solar phase angles greater 65°. On the other hand, the absolute albedo fractional errors of the NEATM grow slower than those of the FRM; for solar phase angles 65° ≤ α ≤ 80°,  the FRM derives on average slightly more accurate albedos than the NEATM, but for α > 80°, the NEATM performs slightly better again, causing the fluctuations in the albedo success rates.

Zoom In Zoom Out Reset image size

The overall trend in the diameter and albedo fractional errors is a consequence of the oversimplified assumptions going into the thermal models. NEATM's failing at phase angles greater than 65° is a result of its lack of night-side emission. As the solar phase angle increases, a larger fraction of the asteroid's night side, which is assumed to not contribute to its thermal emission, is observed. Observed flux densities are greater than modeled flux densities, leading to an overestimation of the object's diameter. The fact that the FRM overestimates the diameter for solar phase angles less than 65° and underestimates it for greater phase angles stems from its limited capabilities to adapt more complex surface temperature distributions. At low solar phase angles, the observer sees the subsolar point, which is usually hotter than any temperature the FRM assumes, leading to an overestimation of the diameter. For high solar phase angles, the subsolar point is below the horizon and a large fraction of the night side is observed, underestimating the diameter. For both models, albedos are overestimated if diameters are underestimated and vice versa.

We take advantage of the systematic effect of the solar phase angle on the median fractional errors to establish statistical correction functions for diameter and albedo estimates derived with the NEATM and the FRM for individual objects. The correction functions are third-order polynomials that are fitted to the median fractional errors shown in Figure 5 (top panel) using a χ² minimization technique; the best-fit polynomial coefficients are listed in Table 1, and the fits are plotted in the top panel of Figure 5. The correction functions are of the form $f(\alpha )={\sum }_{i=0}^{3}{c}_{i}{\alpha }^{i}$ with solar phase angle α measured in degrees. In order to derive corrected diameters (d_corr) from NEATM diameters (d_NEATM), we used the following relation:

$$\begin{eqnarray}&&{d}_{\mathrm{corr}}=\displaystyle \frac{{d}_{\mathrm{NEATM}}}{1+f(\alpha )},\end{eqnarray} \tag{ 1 }$$

where f(α) is the corresponding NEATM diameter function and α the solar phase angle in degrees. FRM and albedo corrections are performed correspondingly. The corrections can be applied to individual objects. The limitations of this approach are discussed in Section 4.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 5 shows the median fractional errors for the NEATM and the FRM after the statistical corrections have been applied to each individual synthetic object on a per-object basis. Over the entire phase angle range, the median corrected diameters and albedos agree with the actual values (fractional error equals zero) within their 1σ envelopes.

Diameters and albedos for about 2500 NEAs have been derived from thermal-infrared observations. The majority of these observations have been performed during the previous decade, namely by the Spitzer Space Telescope (Trilling et al. 2010, 2016) and WISE (e.g., Mainzer et al. 2011a), and diameters and albedos have been derived using the NEATM. Due to spacecraft design constraints, most observations were performed near quadrature, leading to similar constraints on the solar phase angle for NEA observations from both observatories. Both the Spitzer programs (ExploreNEOs, NEOSurvey, and NEOLegacy; see the SpitzerNEOs database http://nearearthobjects.nau.edu/spitzerneos.html) and the WISE observations (Mainzer et al. 2011a) had a significant fraction of their targets observed at solar phase angles greater than 65° but generally less than 90°. According to Figure 4, this implies that for these observations, diameters are likely to be overestimated by up to 25% and albedos underestimated by up to 40% (relative). These effects are of the same order of magnitude as the NEATM accuracies derived by Harris et al. (2011) and hence represent a systematic effect that is not negligible.

We estimate the magnitude of this effect on the ensemble of radiometric NEA diameters and albedos available to date by querying the SpitzerNEOs database (2017 October) and find that about 25% of all Spitzer-observed NEAs were observed at solar phase angles greater than 65°; ∼5% were observed at α > 80° and <0.5% were observed at α > 90°. A similar analysis for WISE-observed NEAs is not possible, as observing geometries for these targets have not been published. However, according to Figure 7 in Mainzer et al. (2011a), the ratios seem to be similar. Hence, we conclude that a small fraction (<5%, the fraction of Spitzer-observed NEAs with α > 80°) of existing NEA diameter and albedo estimates are likely to be affected by systematic effects that are of the same order of magnitude as the expected NEATM accuracy (Harris et al. 2011). We further strongly encourage researchers to include information on the observing geometry with future publications using thermal models.

Our study is limited by the fact that all synthetic asteroids are assumed to be spherical and homogeneous in their properties. Furthermore, we assume no noise model in the derivation of the synthetic asteroid flux densities (flux density uncertainties following Poisson noise have been derived for the χ² minimization in the thermal models (Section 2.3); note that the generated fluxes have not been modified based on these uncertainties). While these simplifications allow us to better probe the effects we wanted to study, it is a gross simplification of real asteroid observations. Hence, we deliberately do not suggest that the fractional error confidence intervals shown in Figures 2–5 provide a measure of the accuracies that the NEATM and the FRM can provide for real asteroid observations. The result that the NEATM fractional error 1σ confidence intervals are of the order of 10% in most cases provides merely a lower limit to the diameter accuracy that can be reached by the NEATM under idealized conditions that minimize uncertainties in the measurement process and the object's light curve. In light of this lower limit, a 20% uncertainty on diameters derived with the NEATM, as derived by Harris et al. (2011), seems to be realistic. A similar argument for the albedo uncertainties is not useful, as albedo uncertainties are mainly driven by the accuracy of the target's absolute magnitude in the optical (Harris & Lagerros 2002) and light-curve effects.

The results of this study are also limited by the applicability of the TPM, which does not consider lateral heat conduction and is only valid in cases in which the diameter of the body is significantly larger than the thermal skin depth, which is usually of the order of a few centimeters (Spencer et al. 1989; Mueller 2007). Hence, the results of this work are valid for meter-sized and larger asteroids. However, future observations of decimeter-sized bodies in close encounters with Earth or so-called minimoons might require additional modeling for a full interpretation of the results that can be obtained from thermal-infrared observations.

Figures 6 and 7 illustrate two cases of synthetic asteroids observed at large and small solar phase angles, respectively, and explain the superiority of the FRM in the former and the NEATM in the latter case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 demonstrates why the FRM provides more accurate model diameters than the NEATM at high phase angles. The NEATM's assumption of zero night-side emission plays an increasing role with increasing solar phase angle. At large phase angles, the NEATM assumes that a large fraction of the observed surface area is on the object's night side and hence cold—for solar phase angles greater 90°, the subsolar point is not visible to the observer. For the synthetic body depicted in Figure 6, the night-side surface temperature is much higher than assumed by the NEATM. Hence, the NEATM tries to compensate for the lack of night-side emission by increasing the model diameter, yielding an overestimation of the object's diameter. In this case, the FRM provides a better approximation of the surface temperature distribution and the object's diameter.

Interestingly, the NEATM provides a much better fit to the object's SED than the FRM. This is the result of the NEATM's ability to modulate the subsolar temperature to fit the observations. However, in this case, a good fit does not equal a good approximation of the object's diameter, as the assumptions on which the model is built describe this configuration poorly. Hence, the simpler FRM provides a better approximation of this situation. The ideal solution to this problem would be to utilize a TPM to solve for the target's diameter and albedo. However, this approach requires additional information on the target's spin axis orientation, rotation period, and other properties. Because this wealth of information is usually not available, we suggest the comparison of NEATM and FRM results, as these will most likely bracket the real diameter and albedo.

In Figure 7, we present the surface temperature distribution of a random synthetic asteroid observed at low solar phase angles. Although the object's surface temperature distribution is nearly iso-latitudinal, the NEATM better describes the object's diameter, which is at least in part due to the fact that the subsolar point (the surface's hottest point) lies within the observed hemisphere, as is assumed for the NEATM for small phase angles. Hence, low phase angle observations are usually better described with the NEATM, even though the temperature distribution is somewhat close to that of the FRM. This observation can be attributed to the fact that real asteroids have non-zero thermal inertias, finite rotational periods, and spin axis orientations that are different from the equator-on assumption on which the NEATM is built, leading to a maximum surface temperature that is lower than predicted by the subsolar temperature.