Evidence of surface heterogeneity on active asteroid (3200) Phaethon

Thermal infrared emission and thermophysical modeling techniques are powerful tools in deciphering the surface properties of asteroids. The near-Earth asteroid (3200) Phaethon is an active asteroid with a very small perihelion distance and is likely the source of the Geminid meteor shower. We estimate and interpret the thermal inertia of this extraordinary asteroid using observations that span ten distinct sightings. The variation in thermal inertia over these sightings is inconsistent with the expected temperature-dependent thermal inertia theorized from radiative heat transfer within the regolith. Thus, we test whether the variation in thermal inertia can be explained by modeling a regolith layer over bedrock and two spatially heterogeneous scenarios. We find that the model in which Phaethon's north and south hemispheres have distinctly different thermophysical properties can sufficiently explain the thermal-inertias determined herein. In particular, we find that a boundary located between latitudes -30 deg and +10 deg separates fine-grained southern latitudes from a northern hemisphere that is dominated by coarse-grained regolith and/or a high coverage of porous boulders. We discuss the implications related to Phaethon's activity and potential association with 2005 UD.


Introduction
Asteroid surface properties can be inferred from the interpretation of thermal inertia, Γ, which is an indicator of resistance to temperature changes (Delbo et al., 2015, and references therein). Thermal inertia is defined as Γ = √ kρc s , where k (thermal conductivity), ρ (mass density), and c s (specific heat capacity) are effective bulk properties of the near-surface. Both the grain size and porosity of the regolith are major factors that influence heat transport (e.g., Wood, 2020), specifically the conduction through contacts between the regolith grains and the radiative transfer within the pores between them (Watson, 1964;Keihm et al., 1973;Gundlach and Blum, 2013;Sakatani et al., 2017;Ryan et al., 2020). Generally speaking, for atmosphereless bodies such as asteroids, the effective thermal conductivity can vary by several orders of magnitude with changes in the grain size (Presley and Christensen, 1997;Presley and Craddock, 2006). Investigations of asteroid surfaces using thermal inertias and so-called thermally characteristic grain sizes are particularly useful when investigating asteroid surface processes and in comparing individual objects to the broad population (MacLennan and Emery, 2021Emery, , 2022. The near-Earth asteroid (NEA) (3200) Phaethon is one of the few asteroids linked to a meteor stream (Williams and Wu, 1993;Kasuga and Jewitt, 2019). The orbital similarity to Geminid stream and observed dust activity during its perihelion passage (Jewitt et al., 2013;Hui and Li, 2016) indicates a link to the annual Geminid meteor shower (Whipple, 1983;Gustafson, 1989;Ryabova et al., 2019). In less than a decade the JAXA DESTINY +1 mission will conduct in-situ flyby observations of Phaethon's dust environment Arai and DESTINY + Team, 2019). Its reflectance spectrum indicates a link to the carbonaceous chondrite meteorites (Licandro et al., 2007;Clark et al., 2010;de León et al., 2012;Kareta et al., 2018;Kareta et al., 2021). However, there is no evidence of hydration expected from the presence of, e.g., phyllosilicates (Takir et al., 2020); which is most likely a result of intense heating during its perihelion passage .
The thermal inertia of Phaethon has been estimated to be 600 ± 200 JK −1 m −2 s −1/2 by Hanuš et al. (2016) and 880 +580 −330 JK −1 m −2 s −1/2 by Masiero et al. (2019). While Hanuš et al. (2016) used a convex shape model derived from lightcurve inversion techniques, Masiero et al. (2019) assumed a spherical shape. While these thermal inertia estimates are consistent with one another, the effective diameter (D eff -the diameter of an equal-volume sphere) estimated by Hanuš et al. (2016) and Masiero et al. (2019) (D eff = 5.1 ± 0.2 km and 4.6 +0.2 −0.3 km, respectively) are in slight disagreement when only the 1-σ uncertainties are considered. It is not clear if the differences are due to shape assumptions, or the distinct observational datasets used in each respective investigation. Radar observations of Phaethon during its 2017 close approach to Earth indicate a spinning top shape with an effective diameter of D eff ∼ 5.5 km (Taylor et al., 2019), which is consistent with the lengths of stellar occultation chords Herald et al., 2020). The discrepancy between the smaller size estimates from thermal infrared and larger size derived from radar observations is yet to be reconciled.
In this work we reassess Phaethon's size and thermal inertia by using a detailed shape model derived from optical lightcurve and radar observations. This shape model is qualitatively similar to the soon-tobe published version, but differs in some respects: slightly different dimensions, and minor variations in some small surface features. After constraining the size using a subset of the observations, we estimate the thermal inertia for each epoch individually. Our TPM results are then analyzed in the context of the heliocentric distance of each observational epoch. In particular, we attempt to discern details about the regolith structure, and/or spatial variability in Phaethon's surface and speculate on how the findings relate to Phaethon's activity and potential relationship with 2005 UD.

Infrared Observations
We collect and use all the publicly-available thermal infrared data for Phaethon. In total, there are 11 distinct observing epochs from 5 telescopes that represent a rich dataset acquired over a time span of nearly 36 years. Specifically, we use thermal observations from Infrared Astronomical Satellite (IRAS), United Kingdom Infrared Telescope (UKIRT), the Infrared Spectrograph (IRS) and Multiband Imaging Photometer for Spitzer (MIPS) instruments on the Spitzer Space Telescope (hereafter Spitzer), Akari satellite (Usui et al., 2011), and 6 sightings made by the Wide-Field Infrared Survey Explorer (WISE) mission and the follow-up NEOWISE survey (Mainzer et al., 2011(Mainzer et al., , 2014. We describe below the telescopes and circumstances of each observation. A summary of all these data and observing geometries is given in Table 1. The IRAS (Infrared Astronomical Satellite) survey observed Phaethon simultaneously at 12, 25, 60, and 100 µm on 1983-Oct-11 (Tedesco et al., 2002). In total, six data points were collected in each of the three shorter wavelength bands, and we omit a single 100 µm observation in order to avoid the uncertainty of emissivity values at longer wavelengths (e.g., Müller et al., 2014).
Photometry at visible and infrared wavelengths were obtained using UKIRT on the nights of 1984-Dec-20,21 by Green et al. (1985). The measurements acquired at < 4 µm contain a significant amount of reflected visible light. To avoid possible contamination from reflected light, we only use the flux values at wavelengths from 4.7 to 19.2 µm. Each photometric data point was acquired at an individual time and thus represent a different part of Phaethon's rotational phase. The Akari astronomical survey detected Phaethon at two wavelengths (9 and 18 µm) on 2006-Sep-23 (Usui et al., 2011). The color-corrected Akari photometry are taken from the Asteroid catalog using Akari (AcuA) archive 2 .
The Spitzer Space Telescope observed Phaethon with the Infrared Spectrograph (IRS) on 2005-Jan-14 (Hanuš et al., 2016) and with the Multiband Imaging Photometer (MIPS; Rieke et al., 2004) on 2007-Jan-01. We use the IRS observations-a high signal-to-noise spectrum that covers wavelengths in the range 5-30 µm-analyzed by Hanuš et al. (2016). This spectrum has an absolute calibration error of at least 5% for each independent datapoint. The MIPS observation consists of one data point at 24 µm (effective wavelength of 23.675 µm), which has not been previously utilized. For the MIPS observations, we download the original image files via the Spitzer Heritage Archive 3 and construct an image mosaic using MOPEX software (Makovoz and Khan, 2005). Aperture photometry is then performed using a radius of 3.5 and a sky annulus defined by inner and outer radii of 20 and 32 , respectively, for which we apply an aperture correction of 2.57 4 . Finally we used a color correction of 0.956 for a 205 K blackbody  to calculate a flux of 60.9 µJy, for which we conservatively assign a flux uncertainty of 5% .
The WISE telescope observed Phaethon first during the cryogenic phase of the mission (Mainzer et al., 2011) in 2010 at a heliocentric distance of 2.317 au (near its aphelion distance of 2.403 au). At this distance, thermal emission dominated the W3 and W4 bands, approximately 12 and 22 µm, respectively. The next five sightings occurred during the reactivated NEOWISE phase of the mission (2015-2019), during which only the W1 and W2-bands, approximately 3.4 and 4.6 µm, respectively, were operational (Mainzer et al., 2014). However, these NEOWISE observations were collected when Phaethon was closer to the Sun where the thermal contribution at these wavelengths was significant (Masiero et al., 2019). Hereafter, we refer to the first sighting during the cryogenic mission as "WISE" and the latter five sightings as "NEOWISE(X)", where the "X" is the sighting number listed in Table 1.

Shape Model Properties
The preliminary shape model used herein (Marshall et al., 2021) was derived using lightcurve observations from 16 apparitions (from 1989 to 2019), and radar data from Arecibo Observatory (from 2007 and 2017) and from the Goldstone Deep Space Communications Complex (from 2017 only). Shape modeling was done using the SHAPE software (Magri et al., 2007). A publication describing the final shape model is in preparation. The modeled rotation period (P rot = 3.603955 h) and spin axis orientation (ecliptic longitude and latitude of λ = 316 • and β = −48.7 • , respectively) are very similar to the convex model of Hanuš et al. (2018). The ratios of the Dynamically Equivalent Equal-Volume Ellipsoid (DEEVE) axes are a/b = 1.08 and b/c = 1.20. This shape model has a mean (volume-equivalent) diameter of ∼ 5.38 km, based on constraints from radar observations. The occultation chords were not directly incorporated into this stage of modeling, but they suggest an average diameter about 4% smaller than the best-fit radar size, consistent within the expected uncertainties. Note that we scale the size of the shape model in our TPM to independently derive a best-fit size from the thermal observations. Phaethon's rotation phase is known from joint lightcurve and radar observations that span over most of the dates in the infrared dataset. This is particularly useful for the NEOWISE(2) and NEOWISE(4) data that consist of only a few datapoints. The exceptions are for the IRAS and UKIRT epochs, in which Phaethon's rotation phase in the 1980s would be accurately known if its rotation period was identical to the value observed from 1994-2019. However, as noted in Hanuš et al. (2016), a single lightcurve from 1989 by Wisniewski et al. (1997) seems to have a different rotation phase 5 , and there are no other published lightcurve observations from before 1994 that might allow us to resolve this discrepancy. Phaethon's rotation period may have changed slightly (possibly due to activity near perihelion), or perhaps there is an unaccounted for error in the 1989 lightcurve. In any case, the uncertainty in rotation phase becomes too large when the information from lightcurves obtained after 1994 are extrapolated backwards by 10+ years. The respective rotation phases in the IRAS and UKIRT observations are thus left as free parameters and adjusted to minimize the χ 2 for each set of input parameters.

Thermophysical Modeling
In this work we present and implement a TPM that incorporates a preliminary non-convex shape model of Phaethon (Sec. 2.2) derived from radar and lightcurve observations. The TPM incorporates one-dimensional heat conduction, small-scale surface roughness, and global self-heating and shadowing effects; the latter two of which are relevant for non-convex shapes. Closely following the Advanced Thermophysical Model (Rozitis and Green, 2011), we modify the TPM presented in MacLennan and Emery (2019) to include features that account for the above effects and call it shapeTPM. The assumed thermophysical properties in shapeTPM are discussed below and are assumed to be the same across the surface, and constant with depth and temperature. We later explore in this work how these assumptions may not hold. However, in the context of thermal inertia determination for individual sightings, it is sufficient to assume these points.
The energy balance at the surface of an asteroid can be expressed as: where E = S /r 2 au is the insolation at 1 au, S = 1367 W m −2 , scaled to the desired heliocentric distance, r au . The first term in Eq. 1 is the amount of direct insolation absorbed by a surface with bolometric Bond albedo A (hereafter Bond albedo) and solar incidence angle, z (the angle between a surface facet's local zenith and Sun-direction). The effective thermal conductivity, k ef f , and surface temperature gradient, ∆T /∆x| surf , determine the subsurface heat conduction. Radiated energy is calculated using the Stefan-Boltzmann constant, σ 0 , bolometric emissivity ε B , and surface temperature, T . The total energy budget includes scattered solar radiation (E scat ) and re-absorbed thermal radiation (E rad , i.e., self-heating), both of which occur among the large-scale facets of the shape model and between small-scale surface roughness elements. The shadowing factor, s, indicates if the surface (shape facet or roughness element) is shadowed by another part of the surface (s = 1) or not (s = 0). For instance, shape facets that exist in topographical lows can be shadowed by other facets and, at a smaller scale, rough elements often block sunlight from reaching other rough elements.
In order to quantify the energy exchange for shape facets and rough elements, we must calculate their view factors, i.e., the amount of radiation that is exchanged between two given areas. The view factor (f j→i from facet j to i) calculation can be approached in different ways, such as the summation via sub-sampled areas of a shape facet (Rozitis and Green, 2011). We choose to perform a line integral around shape facet edges (e.g., Lagerros, 1998), which has been shown to provide improved accuracy and is computationally faster than double-area integration (Walton, 2002), particularly between facet pairs with a shared edge. We use the formulation of Mitalas and Stephenson (1966): Here, view factor subscript denotes that the radiation is received by facet i from facet j, which have areas a i and a j , respectively. The qth and pth edges of facet i and j are divided into N segments of length ∆v i and ∆v j , respectively. The distance between the midpoints of each segment pair is r ij and the angle between the edges is given by Φ pq . Each possible combination of edge pairs is considered. The derivation of Eq. 2 is detailed in Appendix A and the variables from Eq. 2 are shown in Fig. A.9. View factors are used to compute the amount of reflected solar energy prior to the first time step. When only single scattering is considered, then the total amount of reflected insolation on facet i from any other facets in view (j) is (Rozitis and Green, 2011): Similarly, the total amount of thermal radiation received from other facets is: and is dependent on facet surface temperatures and their Bond albedo at thermal wavelengths (A th ∼ 0; Rozitis and Green, 2011). The scattered solar energy received for each facet is computed before a TPM run as all the factors are independent of adjacent time steps. However, because temperatures are not known a priori, the self-heating terms must be solved for at each time step. The assumption of single scattering has been shown to be a sufficient for low Bond albedos (Rozitis and Green, 2011), and because Phaethon's Bond albedo has been independently estimated to be A ≈ 0.048  and A ∼ 0.05 (Shinnaka et al., 2018;Devogèle et al., 2018) we claim that it is an appropriate approximation for our study. As in MacLennan and Emery (2019), the modeling of small-scale surface roughness (i.e. on the order of the thermal skin depth) is approximated by using spherical section craters (Hansen, 1977;Emery et al., 1998). The degree of roughness is determined by the opening angle (γ) of the spherical crater, measured from the center-line of the crater to the edge 6 , and the fraction of surface area that is covered by craters (f R ). Independently changing these two parameters will alter the overall surface roughness, but contribute to mean surface slope,θ, via the following (Lagerros, 1996): For smallθ, different combinations of γ and f R yield similar thermal emission profiles (Emery et al., 1998), however this is not true when either roughness parameter is large and is different for other roughness 6 A hemispherical crater has a γ = 90 • , for example.
representations (e.g., fractal surface; Davidsson et al., 2015). We pair crater opening angles with crater coverage fraction to define a set of three roughness values (in addition to a smooth surface) that we use for TPM fitting. The parameters for low, medium, and high default roughness surfaces are given in Table 2.
The incident energy for crater elements will differ from that of a smooth surface, as seen in the energy budget calculation. First, we calculate the scattered solar and re-absorbed thermal radiation received by the kth crater element from the lth crater elements on the ith facet: and respectively (Emery et al., 1998). As in MacLennan and Emery (2019) we use m = 40 crater elements, which is sufficient resolution for most cases (Spencer, 1990). The above formulations represent an infinite number of reflections within a crater, as opposed to the single reflection that is assumed between any two facets. If the cratered facet is globally shadowed, then none of the crater elements receive any insolation and E crater scat = 0. Crater shadowing is also accounted for wherein the scheme presented by Emery et al. (1998) is used to determine if the kth crater element is shadowed from the Sun by the walls of the crater (s k = 1) or not (s k = 0).
The view factor between a crater element, l, on facet j and the k th crater element on facet i is calculated using the view factor between the facets (f j→i ) and emission angles of each crater element (i.e., relative to a flat surface) as (Rozitis and Green, 2013): where v l,k indicates if there is line-of-sight visibility between the two crater elements. To assess v l,k we implement the same procedure as was used to determine s k , but replace the solar incidence angle with the relevant emission angles. All of the elements are forced to have the same area, because of the constructed geometry of the craters, and no area factor is necessary in Eq. 8. In order to be consistent with the single-scattering approximation, Eq. 8 is calculated only one time: as opposed to an iterative approach (i.e. Rozitis and Green, 2013). In sum, the global scattering and self-heating terms in Eq. 1 are calculated as E scat = E rough ik,scat and E rad = E smooth i,rad + E crater k,rad .

Reflected Light Removal
We must estimate and account for any unwanted contribution of reflected light in the observed fluxes before fitting the TPM fluxes to the observations. Observations from IRAS, UKIRT, Spitzer, Akari and the WISE W3 and W4 bands are dominated by infrared emission from Phaethon's surface, and are therefore not contaminated by any reflected light. However, the short wavelength W1 and W2 of the NEOWISE sightings may have a significant reflected light component. The exact amount depends on the light-scattering behavior and the heliocentric distance, as lower surface temperatures will emit less thermal radiation at shorter wavelengths, allowing reflected light to dominate over any thermal emission. We model the reflected light in W1 and W2 in a similar manner to Alí-Lagoa et al. (2013) as detailed below.
The amount of reflected light at a given wavelength received at Phaethon's surface is calculated from, where F λ 1au is the solar spectrum at 1 au, H V is V -band (0.55 µm) absolute magnitude, R λ is the reflectance relative to 0.55 µm, and Φ λ (α ) is the photometric phase function. In Eq. 10, H V is converted to flux units and is scaled to the desired observer-centric distance. The two-parameter H, G phase function (Bowell et al., 1989) , is used to calculate the flux at a given solar phase angle. We use H V = 14.2 mag (Devogèle et al., 2020) and a slope parameter of G V = 0.06 (Ansdell et al., 2014), the former of which is revised from Tabeshian et al. (2019). Because the phase function is reported for 0.55 µm, we assume that it is the same for W1 and W2. All observed fluxes are scaled to r au = ∆ au = 1 in order to compare across all of the sightings, which represent many different observing geometries. A large source of uncertainty in Eq. 10 is the reflectance value (R λ ). As a B-type, Phaethon's nearinfrared spectrum exhibits smaller reflectance values at larger wavelengths over a range of 0.5−4 µm (Licandro et al., 2007;Kareta et al., 2018;Takir et al., 2020). Specifically, the reflectance value at 3.5 µm is approximately 70% of the value at 0.55 µm (Takir et al., 2020). The reflectance spectrum of Phaethon's throughout the visible and near-infrared bears resemblance to carbonaceous chondrites: CK meteorite spectra (Clark et al., 2010;de León et al., 2012) and heated CI or CM material (Licandro et al., 2007;Kareta et al., 2021). Because Phaethon doesn't exhibit any absorption features in the 3 µm region (Takir et al., 2020) we considered reflectance spectra of dehydrated carbonaceous chondrites presented in (Trigo-Rodríguez et al., 2013). The reflectances of these meteorites in the 3−6 µm range do not significantly (∼2%) vary from a neutral (flat) slope, therefore we claim that an assumption of equal reflectance for W1 and W2 is reasonably justified.
The observed fluxes in the W1 and W2-bands, along with the modeled sunlight flux are shown in Fig. 1 where our assumption of R W 1,W 2 = 0.7 with H V = 14.2 mag is shown as a solid black line. We also find that an error of 0.1 in R W 1,W 2 is equivalent to an uncertainty of 0.15 mag, as shown by the dotted and dash-dotted lines. This translates to a relative uncertainty of ±15% in the fraction of reflected light; meaning that a 10% reflected light component has an uncertainty of ±1.5%. Our nominal reflected light model is very close to that found by Masiero et al. (2019) (shown as a blue dashed line in Fig. 1), who used H V = 14.31 mag and found R λ = 0.8 when using a model that fit reflected and thermal components simultaneously.
Due to the uncertainty with modeling reflected light, and associated model sensitivity, we only use observations that have more than an 85% contribution from thermal emission. In the case of the NEOWISE(4) sighting, both the W1 and W2 observations were dominated by thermal emission, but the W2-band was saturated (Masiero et al., 2019) and is not used in thermal model fits. The red cells in Table 3 that are highlighted in red indicate which W1 and W2 observations that are used in thermal modeling fitting described below.

Data Fitting
In the data fitting procedure, several values of the input parameters are sampled: effective diameter, thermal inertia, and surface roughness. We determine the best-fit effective diameter by searching across all observations, while the thermal inertia is estimated independently for each epoch. For the sightings with  a reflected light component, we multiply the observed fluxes by the modeled relative fraction of thermal emission. For each combination of input parameters the chi-squared statistic is calculated, where σ obs is the flux uncertainty. For comparison between datasets, we use the reduced chi-squared:χ 2 = χ 2 /ν, where ν =(number of datapoints)−(number of model parameters) is the degrees of freedom for the relevant subset of data.

Size and Surface Roughness
We begin our analysis by first placing constraints on Phaethon's size and surface roughness, as these parameters are the same for all observations. To begin, we individually step through each size/roughness combination for each epoch and select the thermal inertia that yields the best-fit solution (i.e., minimizes chisquared value). Only the IRAS, UKIRT, Spitzer-IRS, Akari, Spitzer-MIPS and WISE observations provided meaningful, independent constraints on Phaethon's size. These data represent observations collected at different orbital positions (in particular, relative to perihelion passage), that modeling work of MacLennan et al. (2021) showed may have some effect on the surface temperatures of a hypothetical idealized asteroid. However, orbital heating effects are not a significant factor for this dataset because of Phaethon's oblique spin axis and because we are using both pre-q and post-q observations (see also discussion in Sec. 5). Both the size and roughness are constrained when accumulating the χ 2 values across all the aforementioned datasets.
Next, we sum theχ 2 min values for each size and roughness combination taken across all sightings to generate a set of best-fit solutions (lower panel in Table 4). Aχ 2 threshold criterion is used to establish the upper and lower 1-σ uncertainties MacLennan and Emery, 2019),χ 2 <χ 2 min (1+ √ 2ν/ν). The size and roughness determinations are listed in Table 4 for each individual sighting, for all the sightings in the table, and when grouped by orbital position (pre-q and post-q as defined in Table 1). For example, we find that the best-fit size and roughness combination for the post-q datatsets (IRAS and Akari) is 5.5 ± 0.2 km for a smooth surface. In our discussion (Sec. 5) we explore these results in more detail. For now, we determine that the global best-fit solution is D eff = 5.4 ± 0.1 km for low roughness. With these size and roughness constraints, Phaethon's thermal inertia can be independently estimated for each sighting.

Thermal Inertia
To begin our thermal inertia analysis, we sum theχ 2 values for each size and roughness combination across all sightings. By doing this, the size uncertainty in each subset of observations is accounted for. The modeled thermal inertia value that corresponds to theχ 2 min is taken as the best-fit value. Standard uncertainties are then determined using theχ 2 threshold criterion as was done for the size and roughness procedure. Although the global best-fit corresponds to low roughness (θ = 11 • ), for comparison we estimate thermal inertia across all default roughness values. All sightings except for the IRAS dataset yielded meaningful thermal inertia constraints. Thermal inertia constraints are shown in Table 5 for each roughness value along with the estimated blackbody temperature (T bb ) that is determined by fitting a Planck function to the thermal fluxes.
For airless bodies, the effective thermal conductivity can be expressed as a sum of the radiative (r r ) and solid (k s ) components (Watson, 1964;Gundlach and Blum, 2013;Sakatani et al., 2017;Ryan et al., 2020). The radiative component of conductivity is temperature dependent (k r ∝ T 3 ), which implies that thermal inertia increases with temperature as Γ ∝ T 3/2 , and that thermal inertia is dependent on the heliocentric distance (Γ ∝ r −3/4 au Delbo et al., 2015). Therefore, the surfaces of high-e NEAs such as Phaethon (e = 0.889) are well-suited natural laboratories for observing the temperature dependence of thermal inertia. Rozitis et al. (2018) calculated thermal inertias at different heliocentric distances for three NEAs and found that the thermal inertia of each asteroid increased with decreasing heliocentric distance. They found that the heliocentric dependence of thermal inertia depended on spectral emissivity (or, equivalently, the spectral reflectivity) and that, notably, (1036) Ganymed exhibits a variation that is potentially stronger than expected from theoretical predictions. Using TPM results for over 500 asteroids, MacLennan and Emery (2021) estimated the temperature dependence of thermal inertia of Γ ∝ T 2.74±0.29 , which is also much stronger than the theoretical expectation.
Phaethon's thermal inertia increases with temperature across all sightings and for each default roughness assumed in the TPM. The variation in thermal inertia with heliocentric distance for low roughness can be captured by the power-law (Rozitis et al., 2018): Γ = Γ 0 r ζ au where Γ 0 is the thermal inertia at 1 au and ζ is the variation exponent. The best-fit power law parameters for the entire dataset are Γ 0 = 630 +80 −70 JK −1 m −2 s −1/2 and ζ = −2.45 +0.04 −0.05 . We do not see this as a particularly useful encapsulation of the results because heliocentric distance is only a proxy for temperature, therefore its association to thermal inertia is not direct. Phaethon's thermal inertia variation with temperature (in the range 203 K < T bb < 313 K) across all datasets can be expressed using the power-law formula: Γ = 1.74×10 −9 T 4.65 bb . Such an extreme temperature dependency is highly physically unlikely, yet we find that it is possible to match the datatset using two curves of T 3/2 that are scaled to different intercepts (top panel of Fig. 6). This could be an indicator that there are two distinct regions on Phaethon's surface of differing thermal inertia. We consider such a scenario for Phaethon's surface in the following section. In any case, these results indicate that the theoretical relationship between thermal inertia and temperature is supported. The results of Rozitis et al. (2018) and MacLennan and Emery (2021), should be revisited in future investigations.

Surface Models
We now construct different surface models aimed at interpreting Phaethon's variation in thermal inertia in a more geologic context. The surface models used herein are an oversimplification but useful to model the grain size and porosity using a temperature dependent heat capacity and accounting for radiative heat transfer within the regolith. This is crucial for comparing thermal inertias estimated at different heliocentric distances. However, a major caveat is that these grain sizes and porosities are likely not accurate in an absolute sense, although they are useful for model comparisons.

Homogeneous Surface
We now interpret the thermal inertia dataset in order to gain insight into Phaethon's surface, specifically the size of regolith particles/rocks and possible heterogeneity. First, we attempt to explain Phaethon's thermal inertia variations with a homogeneous surface consisting of grains of a single size that do not vary with depth. Because heat transport within the regolith occurs via solid transfer (through grains and grain contacts) and radiation (in the voids between grains) we use the thermal conductivity model of (Gundlach and Blum, 2013). The radiative and solid thermal conductivity components are summed to model the effective thermal conductivity: where d g is the grain diameter, which we will refer to simply as the "grain size". Here, a regolith grain is a particle with no porosity and therefore only accounts for solid thermal conductivity throughout the material. The grain density (ρ grain ) and grain thermal conductivity (k grain are respectively taken as 2800 kg m −3 (Macke et al., 2019), and 1.5 W m −1 K −1 (Opeil et al., 2010). Poisson's ratio (v) and Young's modulus (Y ) are respectively taken from (Flynn et al., 2018) v = 0.14 and E = 18.9 × 10 8 Pa. The volume-filling factor, ψ, describes the fraction of space occupied by the grains. The non-isothermal factor, f k , is expressed as a function of a dimensionless solid thermal conductivity, Λ s , and accounts for thermal gradients across a single particle (van Antwerpen et al., 2012;Ryan et al., 2020): f k = a 1 tan −1 (a 2 Λ −a3 s ) + a 4 , where a 1 = −0.568, a 2 = 0.912, a 3 = 0.765, a 4 = 1.035, and Λ s = kgrain 4dgσ0T 3 (Ryan et al., 2020). When the particle is particularly large or if the temperatures are high then f k < 1, and if Λ s > 25 then f k is set to unity. The thermal inertia is computed from the effective thermal conductivity, temperature-dependent heat capacity (c s (T )), grain density, and volume-filling factor: Γ = k eff ρ grain c s ψ where, the heat capacity is modeled as: c s (T ) = 2.168 × 10 −1 + 42.581 × 10 −2 T + 4.425 × 10 −2 T 2 − 2.060 × 10 −4 T 3 + 2.853 × 10 −7 T 4 (Opeil et al., 2020). We also follow the procedure in MacLennan and Emery (2022) to estimate the thermal skin depth, l s : for each thermal inertia. The thermal skin depth is the length scale over which the diurnal temperature range changes by a factor of e ≈ 2.718 and the thermal inertia represents the thermophysical properties of the regolith over 1-2l s . As discussed in (MacLennan and Emery, 2022), the grain sizes estimated from thermal inertia may not be physically representative of the surface. Instead they represent a thermally-characteristic grain size that is most related to the thermal conductivity. The grain size and skin depths for each epoch are depicted in the bottom panel of Fig. 2. The estimated grain sizes span an order of magnitude and are positively correlated with the skin depth, which varies by over a factor of two. This correlation is most likely due to the fact that both parameters are strongly influenced by the thermal conductivity. The grain-size estimates are statistically equal to the skin depth for largest values and only the smaller grain-size estimates are smaller than the skin depth. These results can be interpreted in two ways: 1) that the thermophysical properties (e.g., thermal conductivity) vary with depth, or 2) the surface is comprised of a wide distribution of grain sizes and boulders.
In the interpretation of thermal inertia, two assumptions about asteroid surfaces are often made: 1) the thermophyscial properties are constant with depth, and 2) the surface coverage is homogeneous. In what follows, we change these assumptions by first considering a layered two-component regolith model, and then considering a latitude-dependent two-component regolith model. By systematically changing these assumptions, we aim to develop a realistic model that is consistent with Phaethon's thermal inertia estimates.

Layered Model
The layered surface model consists of a regolith over a solid bedrock. The temperature-dependent radiative heat conduction within regolith will cause an increase in the thermal skin depth as the temperature increases. As a result, the thermal inertia (which represents the bulk thermophyscial properties over ∼ 1−2l s ) at smaller heliocentric distances (i.e., higher temperatures) will be larger than a regolith with no bedrock. We test various regolith layer thickness values in order to assess its influence on the effective thermal inertia.
The thermal conductivity of the regolith layer is modeled from Eq. 11. For the underlying bedrock we assign the laboratory-measured values of the thermophysical properties of CM carbonaceous chondrite meteorites. The thermal conductivity and heat capacity values are the same as in the homogeneous model. The porosity, particle size, and thickness of the regolith layer are varied. We express the thickness of the layer in relation to the skin depth of a pure regolith at 300 K (l s300 K ). Surface temperatures are computed for a spherical object at heliocentric distances from 1 to 2.5 au in 0.1 au increments, as well as for the heliocentric distances of the observations themselves (Table 1).
To compare this model to the set of estimated thermal inertias we first calculate the bulk (effective) thermal inertia. This is done by comparing the diurnal temperature variation of the layered model to that of a model that assumes a homogeneous regolith structure. This approach was one of the approaches taken by Biele et al. (2019) to model the effective thermal inertia of a bedrock covered by a layer of regolith. We then fit a power law function, Γ = Γ 0 r ζ au , used in (Rozitis et al., 2018) to characterize the variation in thermal inertia with heliocentric distance for all scenarios of the layered surface model. We compute Γ 0 and ζ for each combination of the three input parameters: layer thickness, grain size and porosity of the regolith layer, and are depicted in Fig. 3.

Quasi-Hemispherical Model
We now construct a model that incorporates latitude-dependent grain size, while also accounting for the temperature-dependence of thermal conductivity and heat capacity. In principle, surface heterogeneity can be studied when an asteroid is observed from many viewing aspects. The sub-solar and sub-observer positions indicate the relative amount of flux that is emitted and observed from an asteroid (e.g., Nugent et al., 2017). A greater fraction of thermal flux will be observed from warmer areas (determined by the sub-solar latitude) with smaller emission angles (a function of the sub-observer point). The plane of sky views in Fig. 4 show Phaethon illuminated by reflected sunlight at each observing epoch. Coincidentally, the heliocentric distance of Phaethon correlates with the viewing aspect and therefore the parts of the surface from which the observed thermal flux is emitted. The thermal fluxes acquired at larger heliocentric distances have a significant portion of flux that arises from south of Phaethon's equator, whereas observations in which Phaethon is closer to 1 au mainly sample the northern hemisphere. In addition to the viewing aspect, the wavelength of observation determined the portions of the surface for which thermal fluxes are weighted (Nugent et al., 2017). Longer wavelengths are sensitive to cooler surface temperatures which-for a low spin obliquity-are found further from the equator. On the other hand, short wavelength observations sample the warmer equatorial region.
Informed by these principles, we calculate distributions of observed thermal flux emitted from Phaethon's surface for each observing epoch. The emitted flux from the fitting procedure is used to produce rotationallyaveraged flux within 10 • latitude bins. The observed thermal flux distributions for each of the ten epochs are shown in Fig. 5, and epochs on which Phaethon was observed using multiple filters or a spectrometer are presented as averages of the relevant wavelengths. For comparison purposes, each panel shows the distribution for the WISE 2010-Jan-7 epoch as a filled grey histogram. It can be seen that the NEOWISE(3) epoch differs from the WISE epoch with more flux observed at the equator because of the shorter wavelengths that were active and detected Phaethon. On the other hand, Phaethon's northern hemisphere comprise more than 90% of the emitted flux observed by Akari and the NEOWISE(4) sighting.
We showed from the thermal conductivity model that smaller grains are a better fit at larger heliocentric distances, and larger grains are more consistent with the smaller heliocentric distances. Thus, we adopt a two-component model by which a larger grain size is assigned to northerly latitudes (ranging from 5 mm to 3 cm) and a smaller grain size for southern latitudes (from 10 µm to 400 µm). Because the upper limit for a modeled grain size is around the thermal skin depth (MacLennan and Emery, 2022) we also consider a boulder-dominated surface for northern latitudes that represents any rock that is greater than ∼ 3 cm. Finally, the porosities of both regions are varied from 10% to 70%. The thermal inertia of porous boulders is therefore modeled, but the effect of radiative heat transfer are not included. The boundary latitude between the two regions is varied by 10 • increments from −40 • S to +40 • N.   Results of this two-component model are shown in Fig. 6. The best-fit grain sizes (and respective upper and lower limits) for southern and northern areas are 40 µm (< 300 µm) and 2 cm (> 5 mm). A boulder-dominated surface (or bare-rock) surface is represented by a "B", indicating widespread presence of grain sizes larger than the thermal skin depth (3.5 cm). No lower limit to the grain sizes in the south and no upper limit to the northern region can be determined. This is due to limitations when attempting to model low thermal inertias using Eq. 11 (for southern latitudes), and interpreation issues when the grain size is comparable to the skin depth (for northern latitudes). These grain size constraints can also be expressed in terms of the thermal inertia at 1 au using the previous model inputs with Eq. 11 and Eq. 12: Γ 0 = 55 +60 −55 JK −1 m −2 s −1/2 and Γ 0 = 800 +400 −400 JK −1 m −2 s −1/2 for the southern and northern areas, respectively.

Equatorial Ridge Model
Phaethon's spinning top shape consists of a characteristic equatorial ridge-in which the local topography is raised relative to a spherical shape. Ridges such as these can form out of the downward flow of material from higher latitudes because the local gravitational potential is lower at the equator. Therefore we test a scenario in which smaller particles are concentrated in an equatorial band and larger particles or boulders are found at higher latitudes. The width of this equatorial band is varied from 5 • to 25 • in 5 • increments and several particle sizes are tested: 50, 60 70, 80, 90 100, 200, and 300 µm. Larger particle sizes are assumed for the areas outside the equatorial band 0.8, 0.9, 1, 2, and 3 cm. Additionally, we use the thermal inertia of solid CM material to approximate a boulder covered surface. We found that the best-fit combination of parameters result in a χ 2 min of ∼ 10 that does not indicate any improvement over the χ 2 min ∼ 1.9 value found for the quasi-hemispherical model.

Size Reconciliation
Using the radar shape model and combining data from five telescopes we are able to place an estimate on Phaethon's size that is consistent with the radar observations (i.e., Taylor et al., 2019). We also quickly note here that the best-fit roughness parameters determined herein (γ = 45 • , f c = 0.5) are nearly identical to that found by Hanuš et al. (2016) (γ = 50 • and f c = 0.5). The size constraints for each data subset listed in Table 4 present some interesting trends. Our best-fit effective diameter measurements that are based only on IRAS, UKIRT, and Spitzer-IRS observations are 6.4 ± 0.2 km, 5.4 ± 0.2 km, and 5.5 ± 0.1 km, respectively. These values are all systematically larger than, but partially consistent with, the Hanuš et al. (2016Hanuš et al. ( , 2018 estimates based on the same respective data sets: 6.0 +0.5 −0.3 km, 4.6 +0.4 −0.2 km, and 5.1 ± 0.2 km. Furthermore, a spherical shape assumption was used to estimate effective diameters of 4.17±0.13 km and 4.6 +0.2 −0.3 km using Akari data and WISE/NEOWISE data, respectively Usui et al. (2011);Masiero et al. (2019). Both these values are markedly inconsistent with our estimations of 5.9 ± 0.4 km and 5.7 ± 0.2 km, for Akari and WISE data, respectively. These results strongly suggests that differences in the shape models used in each work are the cause of the discrepancies in size determination.
As noted in Emery et al. (2014), one important factor in shape input for thermophysical modeling is the tilt of facets relative to the direction of the Sun. To compare across the three shapes used in TPM analyses of Phaethon, we depict the distribution of facet tilts relative to the surface normal at the equator on a sphere in Fig. 7. A sphere has most of the surface area at the equator and tilted outwards (0 • ), whereas the spinning top shapes have higher percent of the surface at ±45 • . Interestingly, the radar model used in this work has a bi-modal distribution of facet tilts because of the equatorial ridge. Although the convex shape model exhibits similar peaks, it also has a third peak at around −15 • . We suspect that this is most likely due to the lack of a pronounced equatorial ridge on the convex model.
The uncertainty in the z-axis of an asteroid shape, which can differ between convex lightcurve shapes and radar shapes, can affect the best-fit thermophysical model parameters (Rozitis and Green, 2014). To assess possible complications we compare the shape model dimensions of the convex shape model of Hanuš et al. (2016) and radar shape model in this work. Using the maximum and minimum projected areas we  Hanuš et al., 2016Hanuš et al., , 2018 and radar shape (dash-dotted red) compared to a sphere (solid black) calculate the axis ratios of an area equivalent ellipsoid (a/b and b/c): a/b radar = 1.09, a/b convex = 1.10, b/c radar = 1.22, and b/c convex = 1.30. This indicates that stretching the radar shape by 7% along the z-axis will force the same area-equivalent b/c as the convex shape. This difference between shape models is drastically smaller than the case of (1620) Geographos as studied in Rozitis and Green (2014). Using a radar model with the z-axis stretched by 7% applied to Spitzer-IRS spectrum we find best-fit TPM parameters indicating a slightly smaller size of D eff = 5.2 km and somewhat larger thermal inertia of 480−560 JK −1 m −2 s −1/2 . Unlike Geographos, the radar observations of Phaethon partly viewed the pole (Taylor et al., 2019), and the shape model used herein is partly based on lightcurve observations. Therefore, we posit that any possible uncertainty in the z-axis is insignificant.
Lastly, we note that the relative differences among the independent diameter estimates from each of the observational sightings are consistently offset from another regardless of what shape model is used. For example, the IRAS size estimates are larger than the value from Spitzer-IRS data both in this work and in Hanuš et al. (2016). We show that diameter estimates using the IRAS and Akari observations, which were acquired post-perihelion, and the WISE sighting before Phaethon reached its aphelion point are all larger than the estimates when using the UKIRT and Spitzer-IRS data, which were acquired before perihelion ( Table 4). The relative offsets in size estimates are consistent with the theoretical prediction by MacLennan et al. (2021) that post-perihelion temperatures are larger than what would be calculated from a diurnal TPM. In this scenario, a period of extreme heating during perihelion passage for a surface with a finite thermal inertia introduces a non-negligible seasonal thermal phase lag. The energy absorbed during a short perihelion passage can thus affect temperatures throughout the orbit; mainly at larger heliocentric distances.
However, a direct comparison to the model predictions of , which were made assuming a spherical object with zero spin obliquity, are not appropriate. Specifically, the fact that Phaethon's hemispheres have different heating periods that are roughly symmetric with respect to perihelion means that seasonal effects are not as extreme compared to a zero obliquity case. By using both pre-and post-perihelion observations at small heliocentric distances, our analysis inherently mitigates any possible unaccounted for seasonal temperature (and effective diameter) bias. Nonetheless, future work is necessary to investigate biases in the sizes among NEAs with high orbital eccentricities.
The bulk density of Phaethon was estimated to be 1670±470 kg m −3 by Hanuš et al. (2018) using a convex shape model with an nominal effective diameter of 5.1 km and thermal inertia of 600 ± 200 JK −1 m −2 s −1/2 . Because the diameter and bulk density affect the orbital drift (i.e., the Yarkovsky Effect Hanuš et al., 2016;Vokrouhlický et al., 2015), in a similar way, we can use the 5.4 km effective diameter to scale the previous density estimate to be 1580 ± 450 kg m −3 . This bulk density is larger than that of (101955) Bennu and (162173) Ryugu (∼ 1190 kg m −3 Watanabe et al., 2019)), which are also classified as B-types. Taking CM chondrites as a compositional analog for Phaethon, this bulk density implies a macroporosity of ∼ 56%, which is comparable to these two other carbonaceous asteroids . Along with this macroporosity estimate, Phaethon's top shape and fast spin rate have implications for its internal structure and formation history in comparison to Ryugu and Bennu (e.g., Barnouin et al., 2019) that extend far beyond this work.

Surface Heterogeneity and Implications for Activity
From our model analyses of Phaethon's thermal inertia we have shown that the northern and southern hemispheres have different thermophysical properties. Specifically, we can conclude that the thermal inertia of northern latitudes from ∼ 10 • to +60 • is different than southern latitudes from −30 • to −60 • . The difference in thermal properties could indicate differences in thermally-characteristic grain size, the abundance of boulders, or some combination of the two. The equatorial region from ∼ −30 • to 10 • may likely be a mixture of the northern and southern terrains. On the other hand, we cannot be certain about the surface properties for polar latitudes (greater than ±60 • ), in each hemisphere because very little thermal emission was observed to come from Phaethon's these latitudes. There could be a significant population of larger boulders that evade detection at southern latitudes, for example, or a possibility of fine-grained regolith at the north pole.
The constraint on the latitude boundary in our quasi-hemispherical model approximately ranges from −35 • to +15 • . It is probable that the equatorial ridge on Phaethon represents the transition from the southern fine-grained surface to the more rocky northern latitudes. Yet, it is unlikely that the transition exists as a sharp step from one hemisphere to another and may be comprised of an intermediate mixture of material from both hemispheres. It's also interesting to note that the equatorial ridge is not constant in latitude. It is conceivable that material flows that have built up the ridge are distributed unevenly. Evidence of variation in the degree of linear polarization with rotation could indicate this, but the data are not specific enough to make such a conclusion.
From our modeling in this work we cannot rule out a scenario in which Phaethon's surface is both heterogeneous and layered. For example, the southern hemisphere could consist of a finer-grained layer over a coarse-grained or boulder-like bedrock. Because we do not have observations of southern latitudes at smaller heliocentric distances, it is not possible to test such a model. Alternatively, differences in the temperature-dependence of thermal conductivity in the surface could increase or decrease the perceived thermal inertia differences between the hemispheres. For example, if the temperature dependence were weaker than that predicted by radiative heat transfer mechanism, then the difference(s) in thermophysical properties would need to increase in order to explain the thermal inertia observations presented herein. In the very unlikely case that the temperature dependence of thermal conductivity is a stronger function than theoretically predicted, then the differences would be lessened.
Several activity mechanisms have been proposed explicitly for Phaethon, or are reasonably plausible based on the extreme thermal environment (Jewitt and Li, 2010): volatilization of some elements (Springmann et al., 2019;Masiero et al., 2021), thermal fracturing (Molaro et al., 2020), meteoroid collisions (Szalay et al., 2019), radiation pressure (Bach and Ishiguro, 2021), electrostatic lofting (Kimura et al., 2022), and thermal destruction of minerals (Lisse and Steckloff, 2022). Despite the fact that Phaethon's hemispheres experience similar heating profiles over an orbit , we expect that hypothetical volatilization of elements would be associated with southern latitudes where there are more small particles that are readily liberated from the surface. A higher thermal inertia for Phaethon's northern hemisphere could indicate a high surface coverage of boulders. The thermal fracturing mechanism is most efficient for rocks that are larger than the thermal skin depth (Molaro et al., 2020) and near the equator ; we expect it to be more relevant in the northern equatorial region that has a high thermal inertia and a larger diurnal temperature range. On the other hand, smaller grains that are likely found in the southern latitudes are easier liberated from an asteroid surface. Instead of thermal fracturing, it's possible that small particles are lost via radiation pressure and/or electrostatic forces at small heliocentric distances (e.g., Bach and Ishiguro, 2021;Kimura et al., 2022). Certain elements and/or minerals are known to be volatile at large temperatures (Springmann et al., 2019;Masiero et al., 2021;Lisse and Steckloff, 2022), and thus could drive dust activity at small heliocentric distances. A volatile mechanism is consistent with the relatively low abundance of sodium found among Geminid meteors Abe et al., 2020).
The DESTINY + spacecraft will not, unfortunately, visit Phaethon during its perihelion passage, but instead at 1 au where no dust activity has been detected Ye et al., 2021). Hypothesized activity mechanisms can be investigated with the information collected by DESTINY + Dust Analyzer (DDA) instrument Arai and DESTINY + Team, 2019) and images of the surface . The latter could involve identification of recently disturbed sites by mapping Phaethon's surface colors and boulder coverage, and interpretation of morphological features. The DDA will measure the mass, charge, density, and composition of collected particles throughout the entire mission (e.g., Krüger et al., 2019), with the aim of identification of dust that originated from Phaethon. This may be particularly useful for testing the impact-driven activity hypothesis (Szalay et al., 2019).

Spectral Variation
Spectral slope differences in the UV region were noted by (Licandro et al., 2007) to potentially indicate variations in surface properties. Using published spectra and new observations from the 2017 close-Earth approach, Lee et al. (2019) found no clear evidence for surface heterogeneity from spectral slope values. Specifically, no outlier observations were identified across many sampled longitudes of Phaethon, but the authors remarked that spectral slopes showed a slight decrease within the range of uncertainty when northern latitudes became more visible to Earth. This scenario is supported by variations in B − V colors, which show a distinct change when the sub-observer latitude is larger than +20 • (Tabeshian et al., 2019). Spectral observations gathered within a few days of Lee et al. (2019) by Lazzarin et al. (2019) showed a clear decrease in spectral slope from one night to the next. This change was explained to be due to the changing viewing aspect, in which the surface portion above +70 • disappear from view over the observations, implying that the spectral differences are due to latitudinal variation in surface properties.
Because of the lack of absorption features among B-types, the groups in de León et al. (2012) were mainly distinguished by spectral slope with Phaethon representing the bluest (most negative) slope. Linking asteroid spectra to meteorite types, however, can be confounded by factors such as grain size which has been shown to affect the spectral slope. A more negative spectral slope is correlated with an increase in grain size and increased thermal alteration (Johnson and Fanale, 1973;Cloutis et al., 2012). Laboratory spectra of meteorite chips in these meteorite groups are less red (more blue-sloped) than particulates and thermal alteration decreases the spectral slope (Cloutis et al., 2012). Thus, these variations seen in reflectance spectra are consistent with particle size differences. In addition, both red and blue spectral slopes have been measured for CM chondrite Murchison caused by changes in surface roughness (Binzel et al., 2015). Finally, the spectrum of ultra-blue (153591) 2001 SN263 resembles that of a chip of the thermally altered CY chondrite Y-82162 (Perna et al., 2014). Thus, Phaethon's blue slope may be a result of the large grain size of its regolith, surface roughness, thermal alteration, or any combination of the three.

Polarimetric Variation
Linear polarization (P r ) measurements of Phaethon have been gathered and presented in several works to date. Polarimetric observations were gathered in the autumn of 2016 by Ito et al. (2018) and throughout Phaethon's 2017 apparition by Borisov et al. (2018), Devogèle et al. (2018), Shinnaka et al. (2018), and Okazaki et al. (2020). Differences in the phase polarization curve from these two years suggest some level of heterogeneity Shinnaka et al. (2018); Okazaki et al. (2020). We summarize the dates in Table 6 and, when possible, extrapolate P r trends to α = 100 • . Notably, the 2016 and 2017 polarimetric phase curves are offset at large phase angles by roughly 10-12%.
In order to reconcile these differences in polarimetric properties in the context of our findings we calculate the surface distribution of reflected sunlight that was observed during the 2016 and 2017 sightings. For this exercise we use the same dates as the NEOWISE(2) and NEOWISE(4) observations that approximately match the dates of the P r measurements of these works. Akin to the thermal flux distribution depicted in Fig. 5, we show the latitudinal distribution of reflected sunlight observed from Earth on 2016 Oct 2 and 2017 Dec 17 in Fig. 8. There is a clear difference in the latitude coverage of reflected sunlight between the two dates, with the 2016 observations highly skewed towards the northern hemisphere and the 2017 observations covering both hemispheres.
The polarimetric phase curve based on 2016 observations is steeper (i.e., a larger P r variation with changing α ) than the curve constructed from 2017 observations. Because higher P r at larger phase angles are associated with larger particle size, we show that the differences in Phaethon's polarization are consistent with larger grain sizes in the northern hemisphere which are in agreement with the conclusions in this work. It is also interesting to note that rotational variations in P r are on the order of ∼ 5%, which may indicate longitudinal heterogeneity  that was not found in spectrophotometric observations (Lee et al., 2019). Searching for longitudinal variation in thermal inertia is technically difficult and beyond the scope of this work, but the residuals from the quasi-hemispherical model fit (Fig. 6) are possibly consistent with this scenario. As we noted earlier, the equatorial ridge may consist of material transported from higher latitudes, making it a likely location for longitudinal changes in surface properties.

Relationship to 2005 UD
Evidence for an association to (155140) 2005 UD has been suggested on orbital Ohtsuka et al. (2006), spectroscopic Licandro et al. (2007); Kareta et al. (2018); Devogèle et al. (2020), and polarimetric similarities Devogèle et al., 2020). Dynamical arguments (Ryabova et al., 2019) and observations of dissimilarities in near-infrared spectral slopes (Kareta et al., 2021) have been used to refute claims that this asteroid pair formed from the splitting of a larger progenitor body in the past.
The thermal inertia of 2005 UD was determined by Devogèle et al. (2020) to be 300 +120 −110 JK −1 m −2 s −1/2 using NEOWISE data at 3.4 and 4.6 µm. Observations from 2 sightings were included when 2005 UD was at 1.36 au and 1.03 au from the Sun and the thermal emission dominated the W2-band and was significant in the W1-band for the warmer sighting. Because 2005 UD does not have a constrained shape model or spin axis coordinates, it's not currently possible to determine independent thermal inertia for these epochs, as done here for Phaethon. Thus, it is not clear if this reported thermal inertia is more representative of the data collected at 1.36 au or 1.03 au or a weighted combination of both. We suspect that the warmer surface temperatures and better signal-to-noise at 1.03 au may skew the thermal inertia value towards being more representative at this heliocentric distance.
Comparing 2005 UD's thermal inertia to the values derived in this work, one can infer that its surface properties may be a mixture of Phaethon's two components (Sec. 4.3), which have higher and lower thermal inertias. The grain size estimated from its thermal inertia are in the range of 1-10 mm (Devogèle et al., 2020). This value lies between the grain sizes found to best represent Phaethon's northern and southern areas. Interestingly, the polarization measurements of 2005 UD (P r (α = 100 • ) ≈ 50 − 55%; Ishiguro et al., 2022) most closely match the 2016 polarimetry (Table 6), which are more representative of the Phaethon's  (Table 6).
northerly areas (Fig. 8). It is possible that 2005 UD formed out of material ejected from Phaethon's northern and southern hemisphere. Binzel et al. (2015) showed that Bennu's near-infrared spectral slope varied from positive and negative values, which they attributed to fine-grained material in the equatorial region. Spectral mapping analyses of OSIRIS-REx observations later supported this claim. Furthermore, our discussion about Phaethon's spectral characteristics mentions grain size and thermal heating effects on the spectral slope. The maximum temperature that 2005 UD currently reaches is just under 1000 K and its perihelion has been larger in the recent past, suggesting that its thermal history is different from Phaethon over the past several thousand years . Because the surface roughness of 2005 UD is not currently known (Devogèle et al., 2020), it could be a yet unaccounted for factor when comparing the near-infrared spectra of the two unusual objects.

Conclusions
Based on our TPM modeling and analysis of Phaethon's rich thermal infrared dataset we conclude the following: 1. When incorporating the detailed non-convex shape, Phaethon's size as determined from thermal infrared observations (D eff = 5.4 ± 0.1 km) is consistent with the radar-derived size. From this effective diameter we revise the bulk density of Phaethon to be 1580 ± 450 kg m −3 , which is slightly lower than the value reported in Hanuš et al. (2016) of 1670 ± 470 kg m −3 . 2. Phaethon's surface is heterogeneous, with a coarse-grained, or boulder-dominated, northern hemisphere with a higher thermal inertia, compared to southern latitudes dominated by fine grains. The quasi-hemispherical model estimates the north-south boundary latitude to be in the range of -30 • and +10 • . 3. The different surface regions described by the quasi-hemispherical surface model (Sec. 4.3) are consistent with previous evidence from spectral and polarimetric observations (Sec. 5.2.1, Sec. 5.2.2) that Phaethon is heterogeneous. Smaller particles that are likely found among southern latitudes are more easily ejected than larger particles in the northern hemisphere.
Once the view factor from facet j to i (F j→i ) is calculated using the procedure outlined above we use the simple reciprocity relationship, a j F j→i = a i F i→j , (A.4) to calculate the view factor, F i→j , from facet j to facet i. This further saves computational time, as the number of view factors that need to be calculated is effectively cut in half.