Dynamical Evolution and Thermal History of Asteroids (3200) Phaethon and (155140) 2005 UD

The near-Earth objects (NEOs) (3200) Phaethon and (155140) 2005 UD are thought to share a common origin, with the former exhibiting dust activity at perihelion that is thought to directly supply the Geminid meteor stream. Both of these objects currently have very small perihelion distances (0.140 and 0.163 au for Phaethon and 2005 UD, respectively), which results in them having perihelion temperatures of or exceeding 1000 K. NEO population models compared to discovery statistics suggest that low-perihelion objects are destroyed over time by a temperature-dependent mechanism that becomes relevant at heliocentric distances < 0.3 au. Thus, the current activity from Phaethon is relevant to the destruction of NEOs close to the Sun, which most likely has produced meteor streams linked to asteroids in the past. We model the past thermal characteristics of Phaethon and 2005 UD using a combination of a thermophysical model (TPM) and orbital integrations of each object. Temperature characteristics such as maximum daily temperature, maximum thermal gradient, and temperature at varying depths are extracted from the model, which is run for a predefined set of a and e. Next, dynamical integrations of orbital clones of Phaethon and 2005 UD are used to estimate the past orbital elements of each object. These dynamical results are then combined with the temperature characteristics to model the past evolution of thermal characteristics such as maximum (and minimum) surface temperature and thermal gradient. We find that dwarf planet (2) Pallas is unlikely to be the parent body for Phaethon and 2005 UD, and it is more likely that the source is in the inner part of the asteroid belt in the families of, e.g., (329) Svea or (142) Polana. The orbital histories of Phaethon and 2005 UD are characterized by cyclic changes in e, resulting in perihelia values periodically shifting between present-day values and 0.3 au. Currently, Phaethon is experiencing relatively large degrees of heating when compared to the recent 20, 000 yr. We find that the subsurface temperatures are too large over this timescale for water ice to be stable, unless actively supplied somehow. The near-surface thermal gradients strongly suggest that thermal fracturing may be very effective at breaking down surface regolith. Observations by the DESTINY+ flyby mission will provide important constraints on the mechanics of dust-loss from Phaethon and, potentially, reveal signs of activity on 2005 UD.


Introduction
Current knowledge of the observed near-Earth asteroid (NEA) population indicates a deficiency of objects with perihelion distance q < 0.3 au, compared to population prediction models (Granvik et al., 2016). This observation-model mismatch among low-q objects suggests that they are actively destroyed at small heliocentric distances by one or more mechanisms unaccounted-for in the model. Because low-q asteroids are subject to intense solar radiation, their surface temperatures reach nearly 1000 K-the largest for any solid surface in the Solar System. This fact suggests that a temperature-dependent process, in particular, becomes relevant at small heliocentric distances and destroys asteroids on timescales shorter than 1 Myrthe dynamical lifetime of NEAs.
Observational monitoring of the largest low-q asteroid, (3200) Phaethon (q ≈ 0.140 au), reveals consistent and repeated activity during perihelion passage in the form of an extended dust tail Hui and Li, 2016). This consistent mass-loss from Phaethon, along with its orbital similarity to Geminid meteors, is strong evidence that it supplies the Geminid meteor stream (Williams and Wu, 1993;Jenniskens, 2006). It is the first asteroid discovered to be associated with a meteor stream (Gustafson, 1989). The orbital similarity of the smaller (155140) 2005 UD (q ≈ 0.163 au) to both Phaethon and the Geminid stream suggests that it may be genetically related.
Phaethon has an effective diameter of 5 − 6 km with a rotation period of ∼ 3.6 hours (Taylor et al., 2019;Hanuš et al., 2018). Its spectral classification as a B-type (Licandro, J. et al., 2007) combined with numerical simulations (de León et al., 2010) suggest its parent body to be the main-belt dwarf planet (2) Pallas, which is the largest B-type asteroid. The geometric albedo of p V ≈ 0.13 Masiero et al., 2019) is slightly higher than most B-types but consistent with Pallas collisional family members in the main-belt Alí-Lagoa et al. (2013. Phaethon's thermal inertia of 600 ± 200 J m −2 K −1 s −1/2 (Hanuš et al., , 2018 gives a characteristic grain size of 1 − 2 cm. Masiero et al. (2019) estimated the near-infrared to visible reflectance ratio to be p IR /p V = 0.8 ± 0.06, which is consistent with the Pallas collisional family Alí-Lagoa et al. (2013). Using measurements of the semi-major axis drift rate, Hanuš et al. (2018) calculated a bulk density of 1670 ± 470 kg m −3 by modeling the thermal recoil forces off the surface (i.e., the Yarkovsky effect). Anomalous brightening and evidence for a dust tail were observed in NASA STEREO images during Phaethon's 2009Phaethon's , 2012, and 2016 perihelion passages (Jewitt and Li, 2010;Li and Jewitt, 2013;Jewitt et al., 2013;Hui and Li, 2016). Current mass estimates for the dust tail (∼ 3 × 10 5 kg) over the narrow ∼ 2 d window of activity at perihelion suggest mass-loss rates (∼ 0.1 − 3 kg s −1 ) that are insufficient to supply the mass of the Geminid stream Hui and Li, 2016). Non-detection of activity elsewhere in the orbit (Jewitt et al., 2019) thus strongly suggests that a vast majority of the Geminid stream was produced when activity was orders of magnitude greater than currently observed.
The Apollo asteroid 2005 UD (hereafter, UD) has been inferred to be associated with Phaethon on similarity of their orbital parameters Ohtsuka et al. (2006) and colors . R-band photometry were acquired by Jewitt and Hsieh (2006), which was used to construct a lightcurve and estimate the object's rotation rate and elongation. The lightcurve modeling yielded a rotation period of 5.249 hr, with a peak-to-trough range of 0.4 mag, although the authors point out that a rotation period of 6.719 hr cannot be ruled out. From the photometric range an axis ratio of a/b = 1.45 ± 0.06 was determined. Using the ∼5.2 hr rotation period and elongation, the authors also calculated a minimum bulk density of 570 kg m −3 that is required for self-gravitation counteracting the centripetal forces. Masiero et al. (2019) used NEOWISE data to more precisely estimate the size (D ≈ 1.2 km), geometric albedo (p V ≈ 0.14), crater fraction (≈ 0.38), and infrared reflectance (p IR /p V ≈ 1.41). Recently, Devogèle et al. (2020) presented a comprehensive lightcurve dataset confirming a period of 5.235 hr, a polarimetric phase curve that is consistent with Phaethon's, and a thermal inertia of 300 +120 −110 J m −2 K −1 s −1/2 using NEOWISE data. This thermal inertia gives a characteristic regolith grain size of 0.9 − 10 mm (Devogèle et al., 2020). Ryabova et al. (2019) backwards integrated the orbital elements of Phaethon and UD over 5000 yr and found no evidence for common origin, contradicting the conclusion of Ohtsuka et al. (2006).
As one of the most reliably active meteor shower events, the Geminid stream has been extensively observed and characterized using several methodologies (Jenniskens, 2006). Ryabova (1999) estimated the age of the stream to be less than a few thousand years, which was corroborated by Beech (2002) who placed constraints of 1000 − 4000 yr. A series of papers has used radar and video observations to construct a detailed model of the structure of the stream (Ryabova, 2007(Ryabova, , 2012(Ryabova, , 2016(Ryabova, , 2017Ryabova and Rendtel, 2018). The ejection velocity was calculated using the apparent width at 1 au by Ryabova (2016) to be 1 km s −1 . Estimates of the mass of the entire stream range from 1.6 × 10 13 kg (Blaauw, 2017) to ∼ 10 14 kg (Beech, 2002). Consistent with these values, Ryabova (2017) estimated the mass of the stream to be M Gem ∼ 10 13 − 10 15 kg, which is comparable to the mass of Phaethon (M P haethon ≈ 1.46 × 10 14 kg, using ρ bulk = 1670 kg m −3 ). Most recently, Battams et al. (2020) noted a direct detection of trail of dust in Parker Solar Probe images that overlaps with Phaethon's orbit. This dust trail was observed when Phaethon was near aphelion, implying that dust fills the entire orbit (with M Gem ∼ 0.4 − 1.3 × 10 12 kg) and an old age (≥ 150 yr). This mass estimate is close to the upper limit on recent (dating back 40 yr) mass-loss of 10 12 kg by Ye et al. (2018). Monitoring of lunar impact ejecta show that meteoroid impactors of 0.2−2 cm exist in at least two distinct sub-streams (Szalay et al., 2018). This could be taken to indicate multiple mass-shedding events from the parent body separated by a large time period, or ejection events from at least two parent bodies.
The exact cause of Phaethon's activity is unknown, but it is likely that its mass-loss is related to the depletion of low-q asteroids in general. (Jewitt, 2012) details many possible activity-causing mechanisms for active asteroids including: radiation pressure, electrostatic repulsion, sublimation, thermal fracturing, dehydration, rotational instability, and meteoroid impacts. Of these mechanisms, Jewitt (2012) favored one of the temperature-dependent processes (i.e., sublimation and thermal fracturing) for Phaethon because of the nature of activity observed during perihelion passage-which we investigate in this work. The proposed JAXA DESTINY+ 1 mission is scheduled to fly by these two bodies in the upcoming decade (Arai et al., 2018;Arai and Destiny+ Team, 2019). In-situ observations of dust activity from Phaethon, and possible activity from 2005 UD, will serve as crucial indicators as to the nature of activity and ultimate destruction of these bodies.
Rotational instability has been proposed as a mechanism causing activity among some asteroids (Jewitt, 2012). Surface instability occurs when outward centrifugal forces are equal to the self-gravitation at the surface (Scheeres, 2015). Torques caused by anisotropic thermal emission from an aspherical body can cause rotational acceleration, referred to as the YORP effect (Rubincam, 2000). We do not consider mass-loss via YORP to be an intrinsically temperature-dependent process because the centripetal force that is directly responsible for mass-loss is not temperature-dependent. We calculate the critical disruption rotation period for Phaethon and UD using, where a/b is the approximate prolate ellipsoid axis ratio, and ρ bulk is the bulk density. The shape of UD has been estimated to be a/b ≈ 1.4, whereas Phaethon roughly resembles an oblate ellipsoid with a/b ≈ 1. Assuming ρ bulk = 1670 kg m −3 for both objects (Hanuš et al., 2018) we find P crit of 2.56 h and 3.58 h for Phaethon and UD, respecively. Thus, neither asteroid is at high risk for rotational disruption, however we do note that if ρ bulk = 1090 kg m −3 is assumed, would place both objects at their respective critical rotation period. Typical YORP timescales for NEOs are on the order of 10 4 yr to 10 6 yr (Rossi et al., 2009) and the rotational acceleration is proportional to (a 2 √ 1 − e 2 ) −1 . This orbital parameter dependence means that objects with lower perihelia will, on average, spin-up (or down) slower than an NEO with a similar semimajor axis and zero eccentricity 2 . For instance, an asteroid with e = 0.9 will take ∼ 40% longer to change its rotation rate, compared to an object in a circular orbit. It is very improbable that YORP spin up is the cause of disruption for low-q asteroids, which typically have larger e. Granvik et al. (2016) mentioned that solar tidal effects, which could disrupt an asteroid, are relevant only at heliocentric distances r < 0.02 au. Theoretically, the combined tidal and rotational accelerations of a rapidly-rotating body could increase this disruption limit to greater distances. Under solar tidal forces, the disruption distance (given in solar radii) for a rotating spherical asteroid is: where ω rot and ω orb are the rotational and simultaneous orbital rotation rates of the body 3 , G is the universal gravitational constant, and ρ is the bulk density of the Sun (1410 kg m −3 ). Due to the strong cubic dependence of the tidal force on distance to the massive body, tidal disruption is relevant only for very small distances. A few active asteroids in the Main Belt were likely subjected to a disrupting impact (Jewitt, 2012), thus it is reasonable to consider that Phaethon's observed activity and/or contribution to Geminid stream is due to meteoroid impacts. Wiegert et al. (2020) proposed that meteoroid collisions of SOHO comets and near-Sun asteroids like Phaethon are the cause of anomalous brightening and their ultimate destruction. Monitoring of dust impacts on the Parker Solar Probe  have shown that most of the dust production occurs at heliocentric distances around Phaethon's perihelion . Based on the orbital distribution of meteoroids (< 1 m) orbiting within Mercury's orbit, Wiegert et al. (2020) suggest that a impact-driven process is the cause of destruction for asteroids q < 0.3 au. In this scenario, a proposed stable ring of dust at low heliocentric distances supplies the meteoroids. We note that this hypothetical ring of debris may be unstable, as solar radiation pressure may be too great to allow particles to exist for long timescales. Regardless of the existence of this ring, Phaethon's impact likelihood is orders of magnitude greater than a circular orbit of similar semimajor axis (Szalay et al., 2019). Using a model of Phaethon's dust environment, Szalay et al. (2019) predict that impact-driven activity should be highly asymmetric and that most impacts should occur on the daylight hemisphere near the sunrise terminator. The activity of Phaethon has only been observed within a narrow time range of two days from perihelion and has been consistent in magnitude. If a model of meteoroid impacts is to be the preferred explanation, it must explain this observational constraint.
Lastly, Kimura et al. (2014) propose that hydroxylated dust aggregates, because of their large surface area to volume ratio, carry enough charge to be readily lofted and swept away by radiation pressure on sub-km bodies. They suggest that this novel mechanism is responsible for inner debris disk of large particles orbiting Fomalhaut near 1 au (Mennesson et al., 2013). We note here there that, since Fomalhaut is ∼ 16.6 times more luminous than the Sun, the energy received at 1 au is equivalent to the solar insolation at 0.24 au. That distance is very close to the upper limit range of meteors with perihelion distances q < 1 au (Wiegert et al., 2020) and close to the range for which asteroid thermal destruction is expected Granvik et al. (2016).
The discovery and monitoring of main-belt comets (MBCs; Hsieh and Jewitt, 2006)-a subset of active asteroids for which the repeated activity is likely to be sublimation driven (Hsieh et al., 2010(Hsieh et al., , 2011(Hsieh et al., , 2018bblurs the boundary between the concepts of icy comets and rocky asteroids. MBCs are distinguished by their periodic and consistent activity that occurs around perihelion, suggesting that sublimation of volatiles (most likely water ice) is the cause. A similar argument can be made for Phaethon, as it's activity is present during perihelion passage. The clear factor that separates Phaethon from MBCs is the drastically higher maximum surface temperatures and shorter duration of dust ejection. If volatile sublimation is responsible for Phaethon's activity, and the destruction of low-q NEAs, modeling must show that it is possible for water ice to exist somewhere in the body and become activated by rising temperatures at perihelion.
Cyclic insolation changes will cause heterogeneous thermal expansion and a corresponding a stress field within a rock. If the local stress exceeds the material strength then material failure will occur, ultimately resulting in the degradation of the rock. Experimental heating of meteorite samples show that material breakdown can occur on rates consistent with hypothesized regolith generation timescales (Delbo et al., 2014). Evidence from modeling (Molaro et al., 2017;Hazeli et al., 2018) and experimental studies (Delbo et al., 2014) suggests that this process could be more efficient than impacts at breaking down asteroid regolith (Basilevsky et al., 2015), particularly for asteroids that come close to the Sun and that have fast rotation periods (Molaro and Byrne, 2012;Molaro et al., 2015). Thus, thermal fracturing may contribute to the destruction of low-q asteroids such as Phaethon and 2005 UD. A thermophysical model (TPM) can be used to calculate thermal gradients and the efficiency of thermal fracturing during perihelion passage.
The OSIRIS-REx mission to study (101955) Bennu is the first mission to an active asteroid (Lauretta et al., 2019b), although it was not planned as such. Measured surface temperatures across the surface of Bennu were used by Rozitis et al. (2019) to estimate the probability that volatiles were present. The authors found through temperature observations and modeling that small patches in Bennu's polar regions are cold enough for water ice to possibly exist. However the occurrence of particle ejections near the equator are inconsistent with sublimation-driven activity . Rozitis et al. (2019) also considered thermal cycling as the mechanism for Bennu's particle ejection. Particles originated from areas when the local time was in the late afternoon which Rozitis et al. (2019) noted was when a negative thermal gradient was present between a warmer subsurface and cooling surface. Exfoliation features on several boulders on Bennu's surface (Molaro et al., 2020) are clear indicators that thermal fracturing is a prominent mechanism acting on NEAs.
Our aims in this paper are thus: 1) investigate potential temperature-dependent mechanisms causing Phaethon's activity, 2) hypothesize as to the temperature-dependent destruction of low-q asteroids, and 3) speculate on the surface properties of Phaethon and UD ahead of the DESTINY+ mission. To do so, we quantitatively characterize the current and past temperature characteristics of Phaethon and 2005 UD by combining thermophysical and dynamical techniques. The temperature-dependent mechanisms considered in this work-sublimation, thermal fracturing, and dehydration-act differently according to the temperature characteristics. Sublimation of volatiles and material dehydration will become relevant after a threshold temperature is exceeded, and the processes of thermal fracturing is expected to be more efficient when the thermal gradient is increasingly large. We thus model the maximum temperatures and thermal gradients for these two objects. Our results are then used to postulate on the broader context of probable temperaturedependent processes affecting this subset of objects.

Dynamical Modeling 2.1.1. Orbit Computation
Samples of orbital elements for the present epoch for Phaethon and 2005 UD are obtained in a two-step process by utilizing the OpenOrb software (Granvik et al., 2009). We begin by computing the best-fit orbital elements in the least-squares sense, as well as the covariance matrix that describes the uncertainty hyperellipsoid for the orbital elements. We make a conservative assumption for the astrometric uncertainty by using σ RA = σ Dec = 0.5 for all astrometry available. Then we sample the resulting uncertainty hyperellipsoid to generate orbits that reproduce the available astrometry. For the sampling we model a Gaussian distribution centered on the least-squares solution with σ j = 2σ nominal,j where j refers to one of the six orbital elements (a, e, i, Ω, ω, M 0 ). This approach is used in order to focus the orbit sample on the best-fitting orbits while ensuring that a significant fraction of the uncertainty region is covered. Note that we do not automatically accept each sample orbit, but require that each one also agrees with the available astrometry given our conservative estimate for the astrometric uncertainty. In addition, we also require that ∆χ 2 = χ 2 − χ 2 LS < 20.1, where χ 2 corresponds to the orbit in question and χ 2 LS corresponds to the best-fit least-squares orbit. The chosen upper limit corresponds to the 3σ limit for a six-dimensional problem. The dynamical model includes gravitational perturbations by the 8 planets, Pluto, The Moon (DE430; Folkner et al., 2014), and the 25 largest asteroids (BC430; Baer and Chesley, 2017). We account for the relativistic perihelion shift, but do not account for non-gravitational perturbations. We discard observations as outliers if the RA and/or Dec residuals are > 2 .

N-body Orbital Simulation
Our integrator of choice was the SWIFT RMVS integrator, which can efficiently handle a close encounter between massless objects and the massive objects in our simulation (Levison and Duncan, 1994). Since there is no straightforward way to perform backward integrations with the SWIFT package, we adopted negative values for all time parameters and modified the relevant checks in the main loop of the integrator accordingly.
We chose a small time step equal to 10 −3 yr, which is an appropriate value for integrating orbits that reach small heliocentric distances. We stopped the integration after 1 Myr, at which point the divergence of the orbits has made it impossible to make a reasonably accurate assessment of the evolution of the perihelion distance q. In addition, we integrated 100 random orbits for 5 Myr to check whether their dynamical evolution would change qualitatively. We record the orbital elements every 100 yr.
We also investigated the importance of general relativity in the orbital evolution of these NEOs, due to the low q values. We included post-newtonian effects in a Bulirsch-Stoer integrator, incorporating the commonly used relation for relativistic corrections derived by Newhall et al. (1983). Next, we integrated the nominal orbits of Phaethon and UD for 0.5 Myr both with and without relativistic corrections. We found that the effect of close encounters in the secular evolution of both objects is much more significant than the small change in the precession of the argument of perihelion (ω) induced by the relativistic effects. Hence, the results that follow are computed without relativistic corrections.

Thermophysical Modeling
We describe here the methods used to model the temperatures for an asteroid during an entire orbital revolution, which we then apply to both Phaethon and UD. The TPM developed by MacLennan and Emery (2019), which is used to calculate diurnal surface temperatures, is modified in order to account for the unique circumstances associated with the longer heating timescale, time-varying insolation, and changing size of the solar disk throughout the orbital revolution of each object. We refer to the model after these modifications as the orbital TPM (orbTPM) and to the unmodified version of MacLennan and Emery (2019) as the diurnal TPM (diTPM). The orbTPM is run for an entire orbital period instead of for individual positions along the orbit. This feature of the orbTPM accounts for the rapid changes in insolation around each object's perihelion passage and for the deeper thermal wave propagation due to insolation changes throughout the orbit. In Sec. 2.2.2 we investigate the significance and importance of accounting for these phenomena, as well as the possible implications for the interpretation thermal infrared emission of NEAs in high-e orbits.

Model Description
The orbTPM numerically calculates surface and sub-surface temperatures, T , on a rotating spherical object revolving around the Sun. The solution to the one-dimensional heat diffusion (Fourier's Law) equation, is computed using the surface insolation as the only energy input. Eq. 3 is evaluated at discrete depth (x) and time (t) interval by using the Crank-Nicolson finite-difference method (Press et al., 2007). The thermal conductivity (k), density (ρ), and heat capacity (c s ) are assumed to be constant with x and t and are treated as bulk effective properties of the regolith. The product of these three thermophysical properties are commonly expressed as the thermal inertia: Γ = √ kρc s . The energy balance for a surface facet with a solar incidence angle of z i and a bolometric Bond albedo of A, is used as the upper boundary condition for solving Eq. 3. The left hand side of Eq. 4 is the absorbed insolation where the solar flux at 1 au, S , is taken to be 1367 W m −2 and is scaled using the inverse square of the heliocentric distance (r). The right hand side is the energy radiated assuming a blackbody with emissivity ε, with Stefan-Boltzmann constant (σ), and thermal energy conducted into the subsurface. The lower boundary condition is taken to be dT dx = 0, implying zero heat transfer, at a sufficiently large depth-as discussed below.
Additionally, in order to partly account for the Sun being a disk, the surface insolation is altered near local sunrise and sunset of a surface facet. This is done, crudely, by using a simple function to modulate the input energy when the solar incidence angle is less than half of angular diameter of the Sun (δ ): This scheme does not explicitly approximate the Sun as a disk, but does approximate a sunset or sunrise well enough for our application to objects with rotation periods much smaller than their orbital periods. More importantly, the approximation produces reasonably accurate facet temperatures in the polar regions, which is relevant to estimating temperatures over an entire orbit. Surface temperatures are calculated for a set of discrete latitude values as the object rotates about its spin axis and orbits the Sun. This means that we do not explicitly calculate temperatures for a complete set of latitude and longitude values at each point in time. Because the rotation timescales for asteroids are significantly shorter than the orbital timescale, the time variable can be substituted for the longitude when a complete map of surface temperatures is needed. The depth variable, x, is discretized into increments small enough such that the heat transport from one to the next ensures model convergence. We set the total number of depth steps to be 300, in order to model heat transport at large depths -which is many factors greater than the diurnal insolation variation. The thermal skin depth, l s , gives the length over which the temperature variation changes by a factor of e = 2.718 . . . over a specified timescale, τ : Typically, TPMs that are constructed to interpret observations collected over a diurnal timescale (τ ∼ 10 hours) must model subsurface heat transport to a sufficient depth for which there is no thermal gradient (and thus no heat transport). Because we wish to model temperatures for an entire orbit (τ ∼ 10 5 hours), we must account for subsurface heat transport associated with insolation changes throughout an orbital period. Additionally, for the sake of increasing the computation efficiency, we increase the spacing between subsequent depth steps as the depth increases such that ∆x i = ∆x 0 × 1.02 i . Thus, the (i + 1)th depth step is 2% larger than the ith depth step. This can be done without sacrificing accuracy because the thermal variations become exponentially smaller with depth and decreases the necessary number of depth steps. On the other hand, the size of discrete time steps is held constant throughout the orbit of the object and set according to the number of time steps per rotation. The number of time steps per rotation is chosen to be large enough such that for large diurnal temperature changes (i.e., at perihelion) the temperature difference between each time step is small enough to ensure numerical stability. For the most extreme cases the diurnal temperature range can exceed 400 K, for which we find that 1800 time steps per rotation is sufficient. After each orbital period we adjust the temperature profile by a particular value-based off of conservation of energy-in order to achieve faster convergence. The model was run for 20 orbital periods with constant orbital elements, which allowed for temperature convergence (the temperature at a given orbital position is < 1 K different from the previous orbital revolution). A final check for model convergence was performed by ensuring that the orbit-averaged surface temperature for a given latitude was < 1 K different from the temperature at the largest depth step (a consequence of energy conservation).

Implications and Caveats
As a sanity check, we compare the surface temperatures from the orbTPM implementation with surface temperatures computed from the diTPM at a given heliocentric distance. For this demonstration, we model the flux from a set of asteroids with varying e and with a = 1.27 au, A b = 0.05, Γ = 450 J m −2 K -1 s −1/2 , P rot = 4 hours, and zero spin obliquity using both models. We find discrepancies in the maximum and minimum surface temperatures throughout the orbit, as well as the diurnal temperature range at small heliocentric distances. In all cases, the discrepancy is higher for larger e.
The discrepancy in diurnal temperature range, which is mainly influenced by thermal inertia and rotation period of the object, increases from nearly zero at > 1.5 au to nearly 15 K at the perihelion point for e = 0.9 (top right plot in Fig. 1). The behavior is likely due to the storage of thermal energy as the surface experiences higher levels of insolation. The discrepancy is greater after perihelion than when an object is approaching from larger heliocentric distances. Maximum and minimum temperature discrepancies show a trend opposite to that of the diurnal temperature range-temperatures calculated using orbTPM are cooler for small heliocentric distances and vice versa further away from the Sun. As with the diurnal temperature range, the discrepancy in maximum temperatures is larger post-perihelion, than pre-perihelion. Generally speaking, the discrepancies seen here are due to the time lag difference between the propagation of the diurnal and orbital thermal wave and how much insolation is absorbed throughout the orbit. These discrepancies, particularly in the maximum temperature, imply that asteroid diameters and albedos calculated over diurnal timescales-particularly for high-e NEAs-may be wrong by a large fraction.
The orbTPM has a few built-in assumptions that sacrifice real-world accuracy for computational convenience. First, the model assumes that the thermophysical properties of the regolith do not vary with depth, or across the surface of an asteroid. For example, Lunar regolith core samples show that deeply buried material is more compacted, which increases the thermal conductivity, bulk density, and effective heat capacity (Keihm et al., 1973). The result of compacted sub-surface material on daytime (i.e. maximum) surface temperatures and thermal gradients is negligible, because the net heat flow is downward from the surface. Spacecraft missions to asteroids (Robinson et al., 2002;Mazrouei et al., 2014;Sugita et al., 2019;Lauretta et al., 2019a) show that their surfaces are non-homogeneous and that regolith grain sizes (and thus, thermal inertia) can vary by a factor of several. Thus, the orbTPM temperature calculations represent an average of the temperature characteristics on the surface. Hanuš et al. (2016) showed that thermal torques acting on Phaethon will cause a drift in its spin axis. This drift is mostly in the longitude at pericenter (i.e., Figure 8. in Hanuš et al., 2016), or in our model the position of the solstice point, with minimal change in the spin obliquity (around ±10 • ). The uncertainty in both the orbit, spin, and thermophysical parameters make it difficult to accurately predict the spin axis orientation at a given moment in the past, but Hanuš et al. (2016) demonstrated that the solstice point oscillates with a range of ∼ 80 • centered around a value close to the present-day spin orientation. Thus, while the current spin axis direction may be known to high accuracy, it is difficult to estimate, with similar accuracy, the past evolution of spin axis.
Finally, it has been theoretically and experimentally shown that thermophysical parameters such as thermal conductivity and heat capacity are dependent on the temperature of a material (k ∝ T 3 ; c s ∝ T ). In addition, latent heat capacity may be relevant if volatiles exist (albeit deeply buried) beneath the surface of asteroids. At certain temperatures (i.e. when sublimation occurs) the latent heat becomes relevant (e.g., Schörghofer and Hsieh, 2018) and thermal energy will be sequestered in the subsurface, subsequently lowering the surface temperature. The effect of latent heat is the opposite from a situation of a radioactive heat source that contributes thermal energy and would raise the surface temperature at the surface. While these temperature-dependent factors are certainly relevant for the calculation of temperatures, we wish to study here the evolution of the temperature characteristics as the orbits evolve through time. The perihelion distance, for example, will have a much larger influence on the maximum surface temperature (T max ∝ q −2 ) than will thermal inertia or heat capacity.

Orbit Solutions
The entire publicly available astrometric data set for Phaethon comprises a total of 5545 optical observations spanning from 1983-10-11 to 2020-02-20. We correct the astrometry for zonal errors in reference star catalogues with the model developed by Farnocchia et al. (2015) and use only these observations for which corrections can be performed. The resulting orbit solutions are based on 4675 optical observations obtained between 1994-11-29 and 2020-02-20 out of which we discarded 48 observations as outliers. The least-squares orbit solution and the derived orbit sample obtained for Phaethon are barely in a few-sigma statistical agreement with solution #712 computed and reported by the Jet Propulsion Laboratory (JPL; top panels in Fig. 2). The remaining, fairly minor, differences in absolute terms between the orbit solutions are likely a consequence of different astrometric weighting schemes and the facts that we did not account for non-gravitational effects nor did we include radar observations in our solution (only one of the 8 radar observations included in the JPL solution are publicly available). We decided to omit non-gravitational forces in orbit computation because it is not clear how comet-like non-gravitational forces should be accurately accounted for in long-term orbital (backward) integrations. That is, we do not know when the observed comet-like activity on Phaethon began or how the activity evolved before the object was discovered.
For UD the entire publicly available astrometric data set comprises a total of 697 optical observations spanning from 1982-11-06 to 2019-10-30. After correcting for zonal errors in reference star catalogues we are left with 668 optical observations over the aforementioned time span. The resulting orbit solutions are based on 666 optical observations which means that we discarded two observations as outliers. The least-squares orbit solution and the derived orbit sample obtained for UD are in statistical agreement with solution #91 computed and reported by JPL using the same input data (bottom panels in Fig. 2).
We will next feed the computed orbits (later referred to as clones) to an orbit integrator and propagate them backwards in time.

Late Dynamical Evolution
As is typical for NEOs in general, also the dynamical history of Phaethon is characterized by numerous close encounters with planets in the inner Solar System: the semimajor axes a of Phaethon clones suffer small jumps and quickly diverge as we integrate their orbits further into the past (Fig. 3).
It is also apparent that all the clones undergo similar variations in eccentricity e during the last 3 − 4 × 10 5 yr, as the mean value of e rises to its current value, a result which is in accordance with Hanuš et al. (2016). When we go further back in time, the eccentricities of individual clones can either increase or decrease, though the mean eccentricity of all clones remains relatively constant until the end of our integration. The apparent stability is due to the complete mixing of the secular phases of the clones which hides the complex behavior of each individual clone. In practice this means that we cannot exactly determine the eccentricity of Phaethon further back in time than 3 − 4 × 10 5 yr.
As far as the evolution of inclination i is concerned, due to close encounters with the planets, the secular phases of the clones' orbits become fully mixed after a few hundreds of thousands of years. In addition, if we track the clones' orbits further back in time, until t = −5 Myr, the mean value of i becomes ∼ 33 • (Fig. 4), which is also the proper inclination of (2) Pallas. Calculating proper elements for NEOs is not possible using the method by (Milani, Andrea and Knežević, Zoran, 1990), because the orbits intersect those of the terrestrial planets. Using the mean inclination as a proxy for the proper inclination is therefore a reasonable alternative.
Integrating the orbits of 2005 UD and its clones we find a qualitatively similar orbital history as Phaethon (Fig. 5) suggesting a common origin. Currently, 2005 UD is at a different secular phase compared to Phaethon (Fig. 6). Considering the difference in secular phase, Hanuš et al. (2016) propose that these objects must have separated earlier than 100 kyr ago. In addition, Ryabova et al. (2019) suggest this fragmentation cannot have happened during the last minimum of q since none lies within the stream of fragments formed by a breakup of the other body. The mean value of i for all 2005 UD clones for t < −2 Myr is also found to be ∼ 33 • (Fig. 4).

Early Dynamical Evolution and Source Regions
Direct backward integrations only allow us to obtain statistically meaningful information about Phaethon and UD clones during the last 10s or 100s of kyr. Hence to study the early dynamical evolution including their likely source regions in the asteroid belt we need an alternative approach which relies on direct integrations of orbits from the asteroid belt then into and within the NEO region.
Let us first assume that either UD is a fragment of Phaethon which, in turn, is a fragment of Pallas or both Phaethon and UD are fragments of Pallas. By analyzing the orbital integrations described by Granvik et al. (2017), that started with an unbiased sample of nearly 80,000 known asteroid orbits in the asteroid belt. Focusing on Pallas-like initial orbits-|a − 2.77 au| < 0.05 au, |e − 0.23| < 0.1, and |i − 34.8 • | < 4 • -we find that it takes a test asteroid 34 ± 28 Myr or 45 ± 29 Myr to reach the NEO region (q = 1.3 au) if its diameter D = 0.1 km or D = 3.0 km, respectively, assuming a typical value for the Yarkovsky drift. Hence the size of the object does not significantly affect the timescale, and the 100-Myr integration is long enough to produce meaningful results.
In total 60 test asteroids with D = 0.1 km reach an entrance route to the NEO region, out of which 30 reach the 5:2 mean-motion resonance (MMR) with Jupiter, 16 reach the so-called outer part of the ν 6 secular resonance, and 13 reach the 8:3 MMR with Jupiter. For one of the test asteroid the delivery mechanism remains unspecified. Among the test asteroids with D = 3.0 km we find 8 that reach the NEO region. Six of these reach the outer part of the ν 6 secular resonance and two reach the 5:2 MMR with Jupiter.
To follow the evolution of the test asteroids to Phaethon-like orbits-|a−1.271 au| < 0.05 au, |e−0.890| < 0.05, and |i − 22.260 • | < 1.0 • -or to UD-like orbits-|a − 1.274 au| < 0.05 au, |e − 0.872| < 0.05, and |i − 28.675 • | < 1.0 • -we use the continuation of the aforementioned orbital integrations in NEO space (Granvik et al., 2016(Granvik et al., , 2018. To account for the destruction at small perihelion distance we consider two different critical perihelion distances, 0.005 au (essentially an impact with the Sun) and 0.076 au (average super-catastrophic destruction distance). As a result of cloning test asteroids entering the NEO region through the outer part of the ν 6 secular resonance or the 8:3 MMR with Jupiter, a total of 204 test asteroids can be linked back to Pallas-like initial orbits. Not a single one of these 204 test asteroids reach Phaethonlike or UD-like orbits when integrated in the NEO region. In fact, there are only 3-4 test asteroids in each scenario that enter the NEO region through the 5:2 MMR complex which covers the region where the Pallas family resides (Tables 1 and 2).
To better understand the source regions for Phaethon and UD we compute the relative likelihoods of the source regions or escape routes from the asteroid belt as defined in (Granvik et al., 2018). The results suggest that it is unlikely that Phaethon or UD would have resided on a Pallas-like orbit in the asteroid belt, because the likelihood for the 5:2 MMR complex is at most around 5% (Tables 1 and 2). An inner asteroid belt source appears much more likely based on the dynamical evidence alone.
The typical lifetimes for test asteroids in the NEO region (q < 1.3 au) until they reach either a Phaethonlike or UD-like orbit range from about 10 Myr to about 80 Myr depending on the entrance route or source region (Tables 1 and 2). For the ν 6 resonance complex and the Hungaria group, that is, the most likely candidates to explain the origin of both Phaethon and UD, there typical lifetimes range from 20 Myr to 80 Myr. Similar timescales are also found for outer-belt sources.
Let us now go see assume that an inner asteroid belt source for Phaethon/UD is the most likely. We now ask: does the relatively high inclination of Phaethon place a limitation on the inclination of the origin? An examination of the initial inclination distribution of test asteroids that eventually evolve into Phaethon-like orbits shows a range from about 0 • to 50 • (left plot in Fig. 7). Excluding the Hungaria and Phocaea groups as sources for Phaethon (right plot in Fig. 7), the inclination range reduces to about 0 • to 30 • with more than 90% of the test asteroids having i < 10 • . The small semimajor axis of both Phaethon and UD implies that these objects have evolved for a relatively long time in the region of the terrestrial planets and they have thus had plenty of time for close planetary encounters that have changed their inclinations. Hence we posit that focusing on the recent average inclination for Phaethon and UD is not very useful for determining their source region or parent body. We cannot exclusively rule out a high-inclination source, however.

Application of orbTPM to Phaethon and 2005 UD
We ran the orbTPM for two sets input parameters which represent the nominal surface properties of Phaethon and UD: the Bond albedo A b , thermal inertia Γ, rotation period P rot , and spin axis orientation, which is expressed in terms of the obliquity ( ) and solstice point (λ ν ), that is, the true anomaly at which northern hemisphere experiences summer solstice (Table 3). For each object, the orbTPM was run for a predefined set of semimajor axes (a) and eccentricities (e): values of a range from 1.1 au to 1.6 au using increments of 0.1 au, and values of e range from 0.50 to 0.95 with increments of 0.05. Additionally, calculations are also executed using the current nominal orbits of both objects. The effect of the spin obliquity on surface temperatures is apparent when comparing the results, particularly near perihelion and aphelion in the polar regions (Fig. 8).
The calculated temperatures resulting from one run of the orbTPM contains a large amount of information about the temperature characteristics of an object such as surface and sub-surface temperatures across several latitudes and for several thousand points along the orbit. Thus, in order to distill this information and best characterize the temperature characteristics we calculate a few diagnostic parameters for each latitude value. In particular, we select the extrema (maximum and minimum) of: 1) T the surface temperature, 2) dT /dx, the thermal gradient (the surface to one diurnal skin depth), and 3) dt/dt the time rate of change in surface temperature. Each of these three parameters is calculated for 1000 points evenly spaced in time throughout the orbit of the object.

Analysis of Current Temperature Characteristics
The surface properties and current orbital configuration of Phaethon and UD result in both similar and contrasting temperature characteristics. We study various parameters as a function of true anomaly (ν), which is not linear with time or heliocentric distance. Since UD's spin axis is not constrained well enough, we analyze the case for which solstice occurs at perihelion in order to compare to Phaethon and investigate the temperature dependencies on spin axis orientation. The best estimate for Phaethon's spin axis orientation implies a solstice point very close to ν = 270 • , that is, 90 • before perihelion, so we chose this scenario to analyze here.
Generally speaking, the maximum surface temperature reached by Phaethon (∼ 1050 K) is larger than that of UD (∼ 975 K), due to the lower perihelion distance of the former (Fig. 8). While the perihelion distance largely determines the value of surface temperatures, the spin-axis orientation (the obliquity, ) primarily influences the distribution of surface temperatures. For both objects, the maximum temperatures are reached at perihelion. Maximum global temperatures can be found at equatorial latitudes for Phaethon, due to the location of the sub-solar point at perihelion, and our spin-axis choice for UD puts the maximum temperatures near 75 • N. When approaching perihelion a nuanced pattern is observed in which the hottest temperatures happen on either side of perihelion, and a local temperature minimum occurs between these maxima, just before perihelion. Because the northern hemisphere is heated on Phaethon's approach toward the Sun, the nighttime temperatures are higher than the nighttime temperatures on the southern hemisphere, which experiences local sunset at aphelion. This thermal behavior is caused by the interplay between the changes in insolation and illumination angle for each hemisphere.
We express the thermal gradient in terms of the temperature change over the topmost thermal skin depth (l s ). Using k = 0.55 W m −1 K −1 , c s = 560 J kg −1 K −1 , and ρ = 3110 kg m −3 (Gundlach and Blum, 2013;Hanuš et al., 2018) and it's rotation period of 3.6 h we obtain l s = 2.56 cm for Phaethon. The thermal inertia of UD is half that of Phaethon, and with P rot = 5.2 h we obtain l s = 1.53 cm The thermal gradient throughout the diurnal skin depth shows a similar pattern as the maximum temperature history. Extreme thermal gradients occur either at perihelion (for northern latitudes) or briefly before and after perihelion (for southern latitudes) for both objects. Positive and negative values of the thermal gradient represent gradients occurring during the daytime and nighttime, respectively (Fig. 8). The equator exhibits the maximum global thermal gradient regardless of orbital position or axial tilt, as the diurnal temperature range are maximized there. Less extreme thermal gradients occur at higher latitudes and are, in general, larger during the daytime hours. When regions are cast into permanent shadow, such as the northern pole of Phaethon after perihelion, the thermal gradient quickly approaches zero. Conversely, the southern pole region of Phaethon experiences a significant thermal gradient increase when it is exposed to sunlight just before perihelion.
On each object the latitude primarily determines the rates of temperature change. The maximum rates of change are always found at the equatorial latitudes, with the polar regions exhibiting the lowest values. This latitude dependence is independent of spin-axis orientation because the equator will always experience a more pronounced day/night cycle when compared to polar regions. Furthermore, surface temperature rates of change are strongly dependent on the rotation period of the object, with higher rates of change correlating with shorter rotation periods. Both objects experience a similar daily temperature range, but Phaethon's faster rotation rate leads to maximum temperature rates of change which are over 1.5 times that of UD.
We have shown that the solstice point will influence the timing of peak perihelion temperatures whereas the obliquity changes where on the surface the peak temperatures are found, particularly in the polar regions. A solstice point at perihelion results in one of the polar regions existing in permanent shadow, which will drastically reduce the total insolation that it receives over an orbit, compared to the rest of the surface and to other spin orientations. This effect can be seen by computing the orbit-averaged temperatures across the surface (Fig. 9). Permanent sunlight at perihelion for the northern polar regions results in higher temperatures at perihelion, and lower orbit-averaged temperatures. Comparing the orbit-averaged temperatures between Phaethon and UD, we can see that the solstice point is more important that the spin obliquity (middle panel of Fig. 9). The southern polar region shows a drastic drop in the peak temperature, due to permanent shadowing at perihelion, but little change in the minimum temperature. Comparing Phaethon to UD we can see that a larger obliquity has a greater influence on the magnitude of the effects described here.
Lastly, we extracted the sub-surface temperature for different latitude bins at various multiples of the annual skin depth. The temperatures at these depths are controlled by the changes in insolation over an entire orbit. Thus, most of the time the temperatures remain constant with an increase during perihelion passage. In Fig. 10 we track this temperature dataset over the 20 days proceeding and leading up to perihelion passage. As expected, the minimum temperatures at each latitude are very close the orbit-averaged temperature value (Fig. 9). Northern latitudes, because they are sunward-facing during perihelion approach, are warmer than southern latitudes until a few days after perihelion. Shallower depths exhibit a larger temperature range, with the timing of peak temperatures occurring closer to perihelion. Temperatures in the southern hemisphere begin to increase right around perihelion because of the drastic increase in insolation when these latitudes come out of shadow. For northern latitudes, which are cast into shadow at perihelion, the situation is reversed and temperatures reach a maximum at (for shallower depths) or sometime after (deeper in the sub-surface) perihelion.

Thermal History
We now move on to characterize the change in temperature characteristics for Phaethon and UD across dynamical evolution timescales. To do so, we must distill the temperature information calculated from the orbTPM runs into a single value for each orbital configuration (i.e., for a given a and e pair). We do this by taking the parameters presented in section 3.2.1 (surface temperature, thermal gradient, and temperature rate of change) and integrate over the time for which any part of the surface meets a particular criterion (i.e., T max > 800 K), per orbital revolution. We perform this analysis on each of the orbTPM runs of predefined a and e, and interpolate the values for desired orbital elements of each clone.
We compute the history of maximum temperatures of the Phaethon and UD clones for the past 100 kyr, because the orbit evolution is predictable over that time interval for both objects (Fig. 11). The effect that evolving eccentricity has on the maximum temperature is very apparent in the periodic evolution of the largest temperature values. The uncertainty, shown by different clone tracks, can also be seen to gradually increase when there is a planetary encounter-in which case the semimajor axis is changed slightly (4 kyr ago for Phaethon and 13 kyr ago for UD). The results show that both Phaethon and UD undergo periods of cooler temperatures every ∼ 20 kyr during which surface temperatures remain under 1100 and 1000 K (most recently 9 to 13 kyr for Phaethon and 3 to 8 kyr for UD).
In the context of this work it is important that we also consider the relevance of the YORP torque on the spin axis for an object, and its effect on the temperature characteristics and thermal history. YORP spin axis drift can alter both the obliquity and solstice point, at potentially different rates. Such change in the spin axis will alter the temperature characteristics of the object. Additionally, since Phaethon (and possibly UD) is losing mass then the YORP coefficient will be affected. Asymmetric mass loss (i.e. from a localized region) would result in a change in the moment of inertia. Inevitably, this would cause a reorientation of the spin axis as the body will seek out the lowest-possible energy spin state.

Discussion
In order to place our analysis of the current temperature characteristics and thermal history of Phaethon and 2005 UD into context, we discuss possible temperature-dependent mass-loss mechanisms: volatile sublimation and thermal stress cycling. We mention, when possible, the likelihood of observing activity from Phaethon or UD during the DESTINY+ flyby inSec. 4.4. The location of activity and the position of the object in its orbit will provide important observational constraints on the mass-loss mechanism.

Dynamical Origins
The similarity of the orbital elements of Phaethon and UD is unlikely to be coincidental, and the resemblance of their reflectance spectra in visible wavelengths implies that the two bodies were once part of the same object . Furthermore, the similarity of the inclinations of Pallas and Phaethon can be interpreted as supporting the arguments that have been made in favor of the genetic relationship between Phaethon and Pallas (de León et al., 2010). However, the results presented here suggest that a relationship between Phaethon and Pallas is highly unlikely given that the likelihood of an origin for the former (and possibly UD) in the inner asteroid belt is much greater than an origin in the outer belt. We have also explicitly shown that using the inclination to argue for or against certain parent bodies is misleading, because the long time spent in the NEO region have allowed Phaethon and UD to significantly change the inclination they had when entering the NEO region through numerous planetary encounters. Hence the inclination-based argument for Pallas being the parent body of Phaethon is weak (de León et al., 2010). Note also that this long evolution also makes it highly unlikely that Phaethon and UD would have separated before entering the NEO region. Therefore we consider it much more likely than not that Phaethon and UD separated relatively recently, but earlier than 100 kyr ago Ryabova et al., 2019).
Although we cannot rule out the possibility that Phaethon and Pallas would be genetically related, the connection is highly unlikely and we therefore have to look for other parent bodies for Phaethon. Given that the evolutionary NEO model by Granvik et al. (2018) suggests that the most likely escape route for Phaethon is the ν 6 resonance complex (a combination of the ν 6 secular resonance and the 4:1 and 7:2 meanmotion resonances with Jupiter) followed by the 3:1 complex, we focus our search in the inner asteroid belt, that is, the region bounded by the ν 6 secular resonance and the 3:1 mean-motion resonance with Jupiter.
Most of the primitive (C-type and B-type) asteroid families in the inner asteroid belt have nearby secular or mean-motion resonances and could be the source family for an H = 17.5 NEO (Fig. 12). Two of these families, Svea and Polana, are particularly interesting to us. The high-inclination (i ≈ 16 • ) Svea family is adjacent to both the ν 6 and 3:1 resonances and is capable of delivering H=17.5 NEOs through both the ν 6 and the 3:1 resonances. If we assume that the Yarkovsky effect has caused the spreading of Svea family members, then its more likely that the ν 6 has delivered Phaethon and UD, as they are both retrograde rotators. While (329) Svea is itself classified as a C-type a notable fraction of the family members are Btype (Morate, David et al., 2019). Morate, David et al. (2019) posit a possible slope-size correlation among primitive asteroids, with smaller bodies (H > 13 mag) having a wider distribution of spectral slopes than larger asteroids (H ≈ 10 mag). This correlation broadens the criteria for spectral similarity among primitive asteroids and thus likelihood of a given object to have originated in a given family-thus Phaethon's origin is consistent with the Svea family despite Svea having a C-type spectrum. We note that the inclination of the Svea family is relatively close to the inclination of Phaethon and UD. The Polana family, with the B-type asteroid Polana the largest member, is a low-inclination (i < 7 • ) family which is thought to be the source family for asteroid Bennu (Bottke et al., 2015). Bennu is actively ejecting mm-to-cm-scale particles at a heliocentric distance ranging from 0.90 au to 1.36 au (see Hergenrother et al., 2020, and references therein). If a mechanism similar to that driving the activity on Bennu is the reason for activity on Phaethon, it would explain the lack of detected activity: the ejected particles are too large, and/or the scattering cross section of ejected particles is too low to be detected by optical telescopes in the visual and near-infrared wavelengths. Furthermore, if both Bennu and Phaethon indeed belong to the Polana family, then other NEOs from the same origin may also exhibit low-level activity at 1 au.

Volatile Sublimation
Detection of hydrogen-bearing volatiles-in the form of adsorbed molecular water (H 2 O) and/or hydroxyl (OH)-has been made on the Moon's polar regions using neutron spectroscopy (Feldman et al., 1998(Feldman et al., , 2000Lawrence et al., 2006), and across the entire lunar surface using near-infrared spectroscopy (Clark, 2009;Bandfield et al., 2018). The southern polar region may harbor water ice deposits in the near-subsurface as evidenced by the LCROSS impact experiment (Gladstone et al., 2010;Schultz et al., 2010), which created a plume in which ≈ 155 kg of water ice was spectroscopically detected . On the other hand, molecular hydroxyl has been detected across the entire lunar surface (Hibbitts et al., 2011;Poston et al., 2015;Lawrence et al., 2015;Grumpe et al., 2019). Additionally, evidence for possible water ice deposits have been detected on Mercury's permanently shadowed polar (northern and southern) crater floors using radar albedo and polarization (Harmon et al., 2001), laser altimeter Deutsch et al., 2017), neutron spectroscopy (Lawrence et al., 2013), and visible imagery (Chabot et al., 2014). These observed lines of evidence, when combined with modeling efforts , strongly support the existence of water ice at Mercury's poles. Features referred to as hollows-typically associated with impact craters-can be found across the surface of Mercury (Thomas et al., 2014;Rubanenko et al., 2019), and have been hypothesized to be formed via sublimation. The locations of hollows are correlated with sun-facing crater slopes and higher latitudes. This correlation mimics the detection of hydrogen-bearing voltiles on pole-facing slopes at latitudes lower than 50 • for the Moon (McClanahan et al., 2015). Lastly, observations of water vapor release from the Lunar surface during meteor shower events are consistent with a hydrated layer at least 3 m-deep overlaid by a ∼ 8 cm layer of desiccated regolith (Benna et al., 2019). The evidence that volatiles existing at or near the surface of the Moon and Mercury despite their large surface temperatures supports the possibility of hydrogen-bearing volatiles in the subsurface of low-q asteroids.
Here, we map the temperature characteristics at perihelion from Fig. 8 onto the radar shape model of Phaethon (S. Marshall, personal communication; Taylor et al., 2019). To do so, we compute the surfacenormal coordinates (as described in MacLennan and Emery, 2019) of each facet of the radar shape model. The temperature characteristic from a given latitude bin from orbTPM are then mapped to the appropriate surface-normal latitude on the radar shape model, performing an interpolation between consecutive latitude bins when necessary. We caution that this method produces an approximation of the temperature characteristics for Phaethon, as shadowing and self-heating effects (e.g. Rozitis and Green, 2011) are not included in the orbTPM, which we leave for future investigations.
The investigations of hydroxyl/water detections on the Moon and Mercury suggest that the conditions may exist for ice to survive buried deep beneath the polar regions of airless bodies orbiting at r ≤ 1 au (Lawrence, 2017). Unlike the Moon and Mercury, Phaethon does not have any permanently shadowed regions due to its high spin obliquity. In general, asteroids with a rubble-pile structure imply that craters with steep walls cannot exist, in which case removes the possibility of having permanently-shadowed patches. For Phaethon, even the high-relief topography at its poles will reach temperatures large enough to vaporize water ice (top row of Fig. 13). Thus if the activity is explained by sublimation, the ice must be buried deep enough where the lifetime is greater than the orbital period, but shallow enough to be activated during perihelion passage.
Evidence of mineral hydration has been detected on Bennu  and it is also possible that Phaethon also contains hydrated primitive material. Although the non-detection of a 3 − µm hydration feature (Takir et al., 2020) does not support the idea that Phaethon's surface contains hydrated minerals or volatiles, we consider the possibility that hydrogen-bearing volatiles could exist at some depth below the surface. A concept called the "buried snow line" was introduced by Schörghofer (2008) to describe the depth at which water ice could exist. To estimate whether ice could survive deep beneath the surface of Phaethon and UD we calculate the orbit-average surface temperature from our model run for the nominal orbit of each object. We found that the orbit-average temperatures were correlated with eccentricity, as expected, but did not vary drastically (±3 K) across kyr timescales-most likely due to the small variations in semimajor axis. We utilize the idealized thermochemical model presented in a series of works by Schörghofer (2008), Schörghofer (2016), and Schörghofer and Hsieh (2018). This series of models were used to calculate the theoretical lifetime of ice on spherical main-belt asteroids, but can be synthesized with our estimates for the temperatures on Phaethon and UD. Using the theory and results laid out in the aforementioned papers, we estimate the buried snow line depth, lifetime of ice survival, and mass-loss rates on Phaethon and UD.
The time spent at these small heliocentric distances is just around two weeks (less than 3% of the time is spent at r < 1 au) with the majority of orbit spent having colder temperatures. The implication of the short heating period is that the temperature deep beneath the surface of low-q objects is low enough to allow for ice to survive across the lifetime of an object, with perihelion temperatures triggering sublimation activity. The thermal history of the polar regions of Phaethon and UD, in particular temperature maxima and minima as well as orbit-averaged temperature. By estimating these values we can assess whether or not volatiles could potentially be stable for certain orbital configurations, and at what perihelion distance sublimation may trigger activity. Using Eq. 6 with τ = 4.522 × 10 7 s, the orbital period of Phaethon, and the regolith thermophysical properties from Sec. 3.2.1 we obtain l s = 1.51 m. Similarly, we use τ = 4.542 × 10 7 s for UD and obtain l s = 0.76 m.
The middle panels of Fig. 13 show that orbit-averaged surface temperatures in Phaethon's northern polar region are ≈ 195 K. Using this temperature, a grain size of 1 cm, and thermal conductivity of 10 −3 W m −1 K −1 , we calculate the desiccation time-the time for all ice to sublimate-to be t des = 50 Myr, assuming a 5 km diameter for Phaethon. This is on the same order of magnitude as the estimated delivery time of Phaethon to the NEO region, as presented in this work. Schörghofer and Hsieh (2018) estimates that half the volume of ice loss happens after 11% of the desiccation time has passed, at which point the buried snow line has will have retreated to ∼ 80% of Phaethon's radius. This means that any ice at depths < 200 m in the north polar region cannot cannot exist after 5.5 Myr. If a 100 µm regolith grains are assumed, then the mean-free path of the water vapour decreases and this depth cutoff changes to < 20 m, which is at least an order of magnitude larger than the annual skin depth. Schörghofer and Hsieh (2018) calculated shallow (1 − 3 m) depths for the survival of ice over the lifetime of the Solar System within main-belt asteroids. Sublimation-driven MBC activity is hypothesized to be initiated by small impacts that expose this shallowly buried ice. Because typical annual skin depths for these objects are on the same order as the depth at which ice can exist, increased solar heating at perihelion is thought to cause sublimation-driven activity of water ice at perihelion. We have demonstrated above that water ice could only exist at depths much larger than the annual skin depth. Thus, this is inconsistent with the main-belt comet-like activity. However, deeply buried ice in Phaethon's core could be activated from an impact that removes tens of meters of regolith.
Using the orbTPM and the model of Schörghofer and Hsieh (2018) we completely rule out the possibility of long-lived water-ice in the topmost 20 meters in the northern pole of Phaethon. A remaining possibility for the existence of a buried hydrogen-bearing volatile is if it were actively generated on short timescales. In this scenario molecular hydroxyl forms via the creation of . OH radicals in crystal defects from solar wind bombardment (Farrell et al., 2015). Evidence for such a process exists for the Moon (Ichimura et al., 2012;Jones et al., 2018) and interplanetary dust particles (Bradley et al., 2014). Zhuravlev (2000) reviewed many experiments which have shown that hydrogen-oxygen bonds can exist for temperatures up to 463 K. This temperature is not reached at an annual skin depth beneath the surface, even for the equator of Phaethon.
Could Phaethon and other low-q asteroids possess a buried layer of polar water-ice (and/or other volatile) that is actively supplied by the implantation of hydrogen supplied via the solar wind? As mentioned before, outgassing during meteoroid impacts on the Moon's surface indicate a hydrated layer a few meters deep (Benna et al., 2019). The temperature at these depths are in the range of 220 − 255 K, which is around the range we calculated for three annual skin depths (195 − 270 K) in the polar regions of Phaethon (Fig. 10). Could a layer of hydrated material that is actively supplied by the solar wind exist at these depths? This hydrogen implantation process could theoretically act on igneous and metamorphic asteroids, which are thought to have lost all their volatiles via radiogenic heating in the early Solar System (Starukhina, 2001). These asteroids are rich in silicates, which have experimentally been shown to form OH bonds when bombarded with protons having energies similar to the solar-wind (Schaible and Baragiola, 2014). Reflectance observations of NEAs have shown that spectral features in the 3-µm region due to OH/H2O are not uncommon, regardless of spectral type and despite their high surface temperatures relative to MBAs. This is consistent with a solar wind supply of hydrogen, but does not indicate that volatiles exist in the subsurface of these bodies.
For primitive asteroids, OH may already be trapped within the crystal lattice of (hydrated) minerals such as phyllosiciates. Laboratory heating of phyllosilicates have shown that the release of OH − occurs over a time period of a few dozen hours. However, experiments that heat a sample rapidly may or may not demonstrate that energy released from within the mineralogical structure is large enough to directly eject mass from an asteroid surface.
Observations and characterization of more low-q objects (spin rate, colors/spectra, thermal inertia, and shape), and monitoring of activity or anomalous brightening will provide important insight into the possible volatile content of these objects. For example, the sungrazing object 322P/SOHO1 has not displayed a coma or tail during its past five perihelion passages, although anomalous brightening has been consistently detected after perihelion. The size and rotation period imply a bulk density of over 1000 kg m −3 which is outside the range of all known comets, and its colors are most similar to a V or Q-type asteroid rather than a typical comet (Knight et al., 2016). The orbital elements of 322P place it close to the Jupiter-family comets, but its semimajor axis is very close to the 2:1 mean-motion resonance with Jupiter which could imply that 322P originated in the main asteroid belt. Could other objects from the asteroid belt be orbiting among known comet populations? The work of Hsieh et al. (2020) show a plausible contribution from the Themis family into this comet population via the 2:1 resonance with Jupiter. With the many challenges and uncertainties associated with dynamical modeling of NEOs, follow-up physical characterization may break the degeneracy between a main-belt or cometary origin for individual low-q objects.

Thermal Stress Cycling
El Mir et al. (2019) constructed a thermomechanical model to estimate the lifetime of a 10-cm-diameter rock on asteroids of various rotation periods and heliocentric distances. This model can be scaled in order to estimate lifetimes of a particular size rock on an asteroid with any orbital configuration. El Mir et al. (2019) showed that thermal fracturing is most efficient for rocks equal to the thermal skin depth, because this is the length scale at which a significant thermal gradient exists. We note that large thermal gradients lead to large thermal stresses when they occur across particles of a size near the thermal skin depth. Our estimate for Phaethon's skin depth of 2.55 cm is quite close to its characteristic grain size of 1-2 cm, however UD's skin depth (1.53 cm) is just above the upper limit of possible grain size range of 0.9-10 mm.
As with the case for Phaethon, extreme thermal gradients for 2005 UD are found among the lower latitudes. It is interesting to note that while positive thermal gradients are less than that of Phaethon's, negative thermal gradients found on UD are similar in magnitude between both objects. However, the large assumed spin obliquity for UD's orbTPM runs caused the timing of the extreme values to occur at ν ± 60 • from perihelion (Fig. 8). According to the results of El Mir et al. (2019), thermal fracturing is less efficient for rocks larger or smaller than the thermal skin depth. From this fact, its possible that thermal fracturing may indeed be an efficient process on Phaethon but less likely so for 2005 UD.
Boulders larger than the thermal skin depth are more susceptible to the formation of large thermal gradients. It is probable that Phaethon and UD have several boulders on their surface that are not captured by thermal inertia estimates. UD's thermal inertia is very close to that of Bennu's: 350±50 J m −2 K −1 s −1/2 , as measured by OSIRIS-REx (Dellagiustina et al., 2019). Analysis of OSIRIS-REx images and thermal data showed an unexpectedly high fraction of boulders with high-porosity (and thus, low thermal inertia) on Bennu's surface (Lauretta et al., 2019a)-not unlike that found for the C-type asteroid Ryugu , which has not been observed to be active. Phaethon's thermal inertia is even higher, increasing the probability that it has many large boulders on the surface. It is entirely possible that Phaethon and UD's surfaces harbor boulders that are currently contributing to active mass-loss.
If the mechanical stress from thermal expansion is converted completely into kinetic energy, particle ejection can occur. The areas most likely to eject regolith on Phaethon are nearest the equator, since this region experiences the greatest temperature as a result of extreme diurnal insolation changes (bottom row of Fig. 13). At perihelion, Phaethon's equatorial region has a diurnal temperature range just over 500 K, maximum temperature rate of change of nearly 60 K h −1 , and temperature difference of 275 K over a skin depth. Jewitt and Li (2010) estimate for the ejection velocity as a function of temperature difference to be v ej = α∆T ηY ρ , where α is the thermal expansion coefficient, η is the efficiency, and Y is Young's modulus. We can use ∆T = 275 K from the results of our modeling (Fig. 8) to estimate an upper limit to the ejection velocity via thermal fracturing of v ej = 16 m s −1 . This value is within the range of estimated ejection speeds , and well above the gravitational escape velocity (1.5 − 2 m s −1 ) of Phaethon (Hanuš et al., 2018).
Thermal breakdown is predicted to be most efficient at lower latitudes where there are large diurnal temperature changes. In a study by Hamm et al. (2019) the authors note that even for highly oblique spin orientations, thermal fracturing should be most efficient at latitudes ≤ 30 • from the equator. This can be seen in the thermal gradient results from our model for both Phaethon and UD (Fig. 8). Breakdown via thermal cycling was proposed by Graves et al. (2019) to be the spectral freshening mechanism for asteroids with q < 1 au. These asteroids exhibit an inverse relationship between perihelion distance and spectral slope (a proxy for space weathering) and Graves et al. (2019) were able to model the relationship using an approximation for the breakdown efficiency. If thermal fracturing is indeed more efficient for low-q objects, then grain sizes from thermal inertia estimates should be smaller for asteroids with lower q. Unfortunately, we only have thermal inertia estimates of Phaethon and UD, which shows a positive correlation with perihelion. However, this could indicate Phaethon's loss of smaller dust particles in the past, increasing the fraction of larger rocks on the surface that raises the thermal inertia. More objects should be used to make any conclusions about a possible correlation between thermal inertia and perihelion distance.
Laboratory experiments that replicate the solar heating environment at 0.15 au have thus far not been performed. We encourage a replication of the heating rates (up to 25 K min) estimated in this work on meteorite samples or analog materials to investigate whether thermal fracturing can eject large quantities of material in low gravity conditions.

Mission Flyby Predictions
Given our investigation into the thermal history and using a detailed shape model of Phaethon from the orbTPM we predict the maximum surface temperatures, and the diurnal temperature range during perihelion passage across the surface. Doing so allows us to glean knowledge about the extreme temperature characteristics and how it may cause surface heterogeneity that can be observed by the DESTINY+ mission.
A small number of carbonaceous chondrite (CI and CM) meteorites have been shown to have experienced secondary episodes of secondary heating (beyond initial aqueous alteration) that resulted in dehydration of their mineral constituents (Nakamura, 2005, and references therein). Additionally, it has been suggested that heating of CK chondrites-which display a range of petrologic types representing different metamorphic textures-is a result of solar radiation (Chaumard et al., 2012). Many heating experiments show that mineral transformation(s) occur at different temperature dependent stages. Laboratory experiments have shown that dehydration of phyllosilicates can occur after heating at 400 K for a week. This is very similar to the exposure time and temperatures experienced by Phaethon and UD when they pass perihelion. Because most or all of the surface exceeds this threshold desiccation temperature, hydrated minerals should not exist anywhere near the surface.

Phaethon's Surface
Here we speculate on the surface features (spectral and morphological) of Phaethon and inference into the nature of activity and mass-loss that can be observed with the instruments aboard the DESTINY+ spacecraft. Spectral and morphological observations can be used to infer on past activity from an active region on the surface.
Using spectra collected to-date, and new observations from the 2017 close-Earth approach, Lee et al. (2019) found no significant evidence for surface heterogeneity. No outliers 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 the northern latitudes were visible to Earth. Spectral observations in which the northern hemisphere were gathered by Lazzarin et al. (2019) within a few days of Lee et al. (2019) and were able to detect a clear change in spectral slope. Lazzarin et al. (2019) notes that northern latitudes > 70 • disappear from view during their observations, implying that the spectral difference is due to the change is viewing aspect. The observations showing a lack of hydration features (Takir et al., 2020) mostly sampled the northern latitudes, suggesting that the slope changes are independent of any potential surface hydration.
Generally speaking, Phaethon's B-type spectrum bears resemblance to both CK meteorite spectra (Clark et al., 2010) and heated CI/CM (Licandro, J. et al., 2007). 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. Specifically, Johnson and Fanale (1973) show that blue slopes can be exhibited by coarse-grained carbonaceous chondrite meteorite samples that show positive slopes when ground into fine-grained powders. Thus, Phaethon's blue slope may be a result of the large grain size of its regolith, and establishing a relationship to a parent body or meteorite type should take the grain size/spectral slope relationship into account. Establishing a clear link between a specific meteorite group and Phaethon would provide insight into the formation and heating of B-type asteroids (de León et al., 2012) and constraints on the efficiency of the solar heating mechanism which is likely responsible for post-metamorphic dehydration of chondritic material.

Phaethon's Activity
Monitoring of the locations and timing of dust activity observed from Phaethon's surface will be indicative of the cause. Direct dust ejection via thermal fracturing will most likely cause dust activity from the equatorial regions, whereas ejection from ice sublimation would originate from the polar regions. Knowledge of the timing of Phaethon's activity, and any low-level dust ejection, is possibly restricted by the brightness limits of remote observations using telescopes. Thus, if dust activity extends beyond the current ∼ 2 days at perihelion then it has implications for the dust ejection mechanism and mass-loss rate.
It is also possible that dust ejected from Phaethon's polar regions could be due to electrostatic and radiation pressure forces. During perihelion passage, the solar illumination conditions could be ideal for the electrostatic charging of dust grains at the day-night terminator. Radiation pressure is relevant for carrying away electrostatically-levitating small dust grains and may be strong enough to push large grains out of Phaethon's gravity if directed near-parallel to the surface (Jewitt, 2012).
Observations of Phaethon's dust tail place the window of activity to at least a few days around perihelion, with a tail mass of ∼ 3 × 10 5 kg and loss-rate of ∼ 0.1 − 3 kg s −1 . Given this tail mass, Jewitt et al. (2019) points out that this is orders of magnitude too low to resupply the Geminid stream which has a mass on the order of 10 13 − 10 15 kg. The photometric observations of Phaethon's dust tail are sensitive to ∼ 1 µm grains, which was used to estimate the mass of the ejected plume and loss rate . Jewitt (2012) points out that the µm-sized grains that are capable of ejection through thermal stress cycling are quickly swept away from the surface via radiation pressure. Modeling performed by Moorhead (2020) suggests that Geminid particles smaller than 3 µm are ejected from the solar system.
On the other hand, Szalay et al. (2018) calculated sizes of observed Geminid meteoroid impacts on the Moon to be 0.2 cm < d grain < 2 cm, suggesting that the Geminid debris stream consists of dust that spans many orders of magnitude in size. Additionally, thermophysical modeling suggests that Phaethon's characteristic regolith grain size is 1 cm < d grain < 2 cm (Hanuš et al., 2018). Combining lunar impact events, visible camera observations, and radar monitoring of meteors Blaauw (2017) calculated the mass index of the Geminid stream to be s = 1.68 ± 0.04, indicating that most of the mass is contained in larger particles. The cumulative mass function can be expressed as log 10 N c = (1 − s) log 10 m + c 1 , where N c is the number of particles larger than mass, m, and c 1 is a constant. Estimates of the mass of Phaethon's dust tail, when assuming that it consists exclusively of ∼ 1 µm particles, range from 2.5 − 4 × 10 5 kg Jewitt et al., 2013). Jewitt et al. (2013) mentions that larger (mass dominant) particles could escape detection in the optical observations. In the following, we explore the implications of such a scenario on the estimated mass ejected during active events.
Here we adopt a power distribution for the dust tail of Phaethon in order to re-evaluate mass estimate of the ejected material. We assume a maximum grain radius of 2 cm, which is slightly less than the size largest-known Geminid (3.6 cm Blaauw, 2017). The characteristic grain radius of Phaethon's regolith from thermal analyses is 2 cm, implying that smaller grains are less common (Gundlach and Blum, 2013). The mass estimate from Jewitt et al. (2013) serves as an anchor point for the size distribution-the integrated mass of 1 µm grains in the tail. Using the mass index function from Blaauw (2017) we extrapolate the mass contribution from particles over the range 1 µm − 2 cm. Doing so, we estimate that the mass-loss of the tail to be ∼ 4 × 10 7 kg. It is important to note here that the dust observed from Phaethon is distinct from the Geminids that are currently interacting with the Earth-Moon system and we are assuming that the size-distribution of ejected particles have not significantly changed over the past 1000 yr. Because larger particles in this size distribution are less abundant, but contribute a larger fraction of the total mass, our estimate is a factor of 2 orders of magnitude larger than Jewitt et al. (2013). Subsequently, the mass loss rate (∆M/∆t) presented here is larger than Jewitt et al. (2013) by the same factor, implying a reduction of the time needed to supply the entire Geminid stream with meteoroids. This time is approximately 1×10 3 yr, which is consistent with the dynamical lifetime of the stream (Ryabova, 2007). From this analysis we claim that the mass-loss rate is at least two orders of magnitude larger than previous estimates.
In-situ observations performed by the DESTINY+ spacecraft will further constrain the mass-loss rate more accurately and confirm whether or not mm-to-cm-scale particles are ejected from Phaethon's surface. For now, we can attempt to extrapolate this mass-loss rate to approximate the total amount of material that Phaethon has shed in the past 1000 yr or so. If we assume that the dust lost rate is held constant but activity only is present for a 2 day period in each revolution we can estimate the mass lost over a given time period. For this, we assume that the semimajor axis is constant but account for changes in eccentricity calculated using our numerical integrations. Using this approach, we compute that 2.7 × 10 10 kg of material has been shed from Phaethon in the past 1000 yr. This value is larger than previous estimates, but still 2-3 orders of magnitude less than the mass of the Geminid stream. Thus, in order to explain the remaining mass of the Geminid stream Phaethon must not be not the only parent body of the Geminids, and/or Phaethon's mass-loss rate was much higher at some point in the past.
Our mass loss rate estimate, if accurate, also has implications regarding the orbital drift modeling of Phaethon. Hanuš et al. (2018) assumed that all of the observed drift in Phaethon's orbit was due to the Yarkovsky effect (thermal photon emission recoil) to estimate its bulk density. Using the mass loss estimates from Jewitt et al. (2013) the authors concluded that the recoil force from Phaethon's activity contributed, at most, to 10% of the orbital drift rate of −6.9 × 10 −4 au Myr −1 . Following (Hanuš et al., 2018), if our mass-loss rate is adopted along with an ejection velocity of 16 m s −1 , the drift rate from perihelion activity is 4.1 × 10 −4 au Myr −1 (nearly 60% of the observed drift rate). As outlined in (Appendix B Hanuš et al., 2018), this contribution assumes that the mass is being ejected in a narrow jet out of the north pole. If the activity originated from the equator and was uniformly ejected in time, for example, then the average recoil force would be negligible. Nevertheless, we strongly encourage a future investigation into the possibility that mass-loss from Phaethon may contribute in a non-negligible way to the orbital drift. If this is the case, then the mass estimate of Phaethon would need to be reassessed. If a significant fraction of the orbital drift is due to the dust activity, then the bulk density (or mass) estimate will in actuality be larger than the Hanuš et al. (2018) value.

2005 UD: An Active Asteroid?
No observations of activity currently exist for 2005 UD-which, if the object does undergo some degree of mass-loss, may be due to its small size and/or low degree of activity. Compared to Phaethon, particle ejection from UD's lower surface gravity field would not require as much force. Detailed physical characterization has been performed using many methods (Devogèle et al., 2020), but UD's shape is not as well characterized as Phaethon's. Furthermore, it is difficult to make accurate predictions on the temperature characteristics due to the uncertainty in the spin axis orientation of UD. If our assumption of a 70 • obliquity is correct, we can use our modeling results to speculate on possible dust activity.
If UD displays dust activity during flyby observations, the locations will be provide important constraints on the mechanism. We have demonstrated, for example, that thermal gradients are maximized for the equatorial region. Thus, greater dust activity observed at these orbital positions would strongly support thermal fracturing as the source of activity. This logic holds true for the assumed larger spin obliquity for UD, but is independent of the solstice position-as this would not alter the thermal characteristics of the latitudes nearest the equator where thermal cycling is thought to be most efficient . Similarly to Phaethon, ice sublimation would theoretically be most likely to occur near the northern pole of UD-if volatiles are somehow present beneath the surface. The orbit-averaged temperature across the surface is mostly dependent on the solstice point (Sec. 3.2.1), which is unknown for UD.
Our dynamical modeling results show that UD has reached perihelion distances of approximately 0.14 au in the past. This value matches Phaethon's current perihelion and the heliocentric distance at which activity occurs. If UD is subject to the same temperature-dependent mass-loss process as Phaethon, then it's possible that activity from UD was prevalent when the perihelion was smaller in the distant past. If true, then UD's original size could have been much larger than today. This has interesting implications regarding the hypothetical separation of UD from Phaethon.
Asteroid pairs-which are two asteroids that exist in orbits closer than can't be produced other than from a single body-are thought to have formed via rotational fission. Theoretically, a smaller section of a rapidly rotating body is ejected due to large centrifugal forces overcoming the strength of self-gravity. The orbits of Phaethon and UD are not currently close enough in orbital space to be classified as a pair, however this could be due to orbital drift due to, e.g., mass loss or planetary encounters. This also raises the possibility that many more asteroid pairs remain unidentified through current methodologies. It is also possible that the tidal forces during a close planetary encounter could provide the means for the formation of an asteroid pair. This is possible, for instance, because the orbital histories of both objects exhibit multiple close encounters with the terrestrial planets.

Conclusions
Following the results of our investigation into the dynamical and thermal history of Phaethon and 2005 UD we conclude the following.
We confirm the results of previous dynamical studies of Phaethon and UD, but our interpretation is markedly different. We do not find compelling evidence that Pallas would be the parent body of these objects. Rather our results points to a parent body in the inner asteroid belt which contains families such as Svea and Polana.
We also confirm the finding of previous dynamical studies of Phaethon (e.g., Hanuš et al., 2016) that the most recent minimum q value occurred ∼ 2000 yr ago and may coincide with the birth of the Geminid meteor stream (Fig. 3).
Volatiles, even those buried up to 200 m below the surface, are highly unlikely to have survived 5 Myr and no ice can survive anywhere in Phaethon longer 50 Myr. Thus if activity-driving, hydrogen-bearing volatiles exist beneath the surface of Phaethon, they must be actively generated somehow-perhaps by similar process(es) hypothesized to supply ice on the Moon and Mercury.
The regoliths of Phaethon and 2005 UD are particularly susceptible to thermal fracturing, which we theoretically predict can eject particles up to 2 cm from the equatorial region. Observations from the DESTINY+ mission will aid in resolving the cause of Phaethon's perihelion activity.
Assuming a particle size distribution for Phaethon's dust tail, we estimate that the mass-loss rate is at least two orders of magnitude larger than previously calculated.     : Inclination distribution for the initial orbits in the asteroid belt that eventually evolve to Phaethon-like orbits. The left plot includes test asteroids from all source regions and escape routes and for the right plot we removed test asteroids belonging to the Hungaria and Phocaea groups.       Granvik et al. (2018), and the time (∆t) from q = 1.3 au until (a, e, i) are similar to Phaethon's current orbit. For ∆t we provide the 10%, 50% (that is, median, in bold), and 90% quantiles, and N is the number of test asteroids used for the analysis. Two different assumptions are used for the critical distance at which asteroids are destroyed close to the Sun: q = 0.005 au and q = 0.076 au (see text for further details).  Granvik et al. (2018), and the time (∆t) from q = 1.3 au until (a, e, i) are similar to 2005 UD's current orbit. For ∆t we provide the 10%, 50% (that is, median, in bold), and 90% quantiles, and N is the number of test asteroids used for the analysis. Two different assumptions are used for the critical distance at which asteroids are destroyed close to the Sun: q = 0.005 au and q = 0.076 au (see text for further details).   Huang et al. (2020)