ASTERIA—Asteroid Thermal Inertia Analyzer

Thermal inertia estimates are available for a limited number of a few hundred objects, and the results are practically solely based on thermophysical modeling (TPM). We present a novel thermal inertia estimation method, the Asteroid Thermal Inertia Analyzer (ASTERIA). The core of the ASTERIA model is the Monte Carlo approach, based on the Yarkovsky drift detection. We validate our model on asteroid Bennu plus 10 well-characterized near-Earth asteroids (NEAs) for which a good estimation of the thermal inertia from TPM exists. The tests show that ASTERIA provides reliable results consistent with the literature values. The new method is independent of TPM, allowing an independent verification of the results. As the Yarkovsky effect is more pronounced in small asteroids, the noteworthy advantage of ASTERIA compared to TPM is the ability to work with smaller asteroids, for which TPM typically lacks input data. We used ASTERIA to estimate the thermal inertia of 38 NEAs, with 31 of them being sub-kilometer-sized asteroids. Twenty-nine objects in our sample are characterized as potentially hazardous asteroids. On the limitation side, ASTERIA is somewhat less accurate than TPM. The applicability of our model is limited to NEAs, as the Yarkovsky effect is yet to be detected in main-belt asteroids. However, we can expect a significant increase in high-quality measurements of the input parameters relevant to ASTERIA with upcoming surveys. This will surely increase the reliability of the results generated by ASTERIA and widen the model’s applicability.


INTRODUCTION
Despite their great importance, knowledge of the physical properties of asteroids lags far behind the current rate of their discoveries.To better understand asteroids' properties, one of the most critical parameters to estimate is their surface thermal inertia.For instance, the value of thermal inertia can give essential constraints on the type of material present on the surface.However, as the asteroid's thermal inertia varies in a considerable interval of values spanning more than an order of magnitude, it cannot be reliably assumed a priori.Yet, the determination of the surface thermal inertia is a challenging step.
Data provided by space missions to asteroids and comets are extremely valuable information.Among them, NASA's Dawn space probe visited asteroids Vesta and Ceres (Russell et al. 2015), while ESA's Rosetta mission performed a detailed study of comet 67P/Churyumov-Gerasimenko (Glassmeier et al. 2007).The most recent examples are two robotic sample-return projects: JAXA's Hayabusa2 and NASA's OSIRIS-REx, exploring Ryugu and Bennu, respectively (see, e.g.Watanabe et al. 2019;Lauretta et al. 2019).The level of detail, reliability, and precision of the thermal inertia values provided by these missions is outstanding (see, e.g.Capria et al. 2014;Rognini et al. 2020;Grott et al. 2019;Rozitis et al. 2020).However, such data is available only for a small sample of objects, and it is difficult to extrapolate these results to other objects or to derive global properties valid to a population level from a small set of well-known asteroids.Nevertheless, the spacecraft visits provided data that are key for testing and validating the other models used for thermal inertia estimations.
Typically, the thermal inertia is estimated by thermophysical models, making use of thermal infrared observations (see, e.g.Delbo' et al. 2015, and references therein).Such observations, however, are available for a limited number of asteroids.Until several years ago, the thermal inertia was estimated for a very small number of asteroids (Delbo' et al. 2015).The situation started to change recently, primarily thanks to the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) and its post-cryogenic survey extension known as NEOWISE (Mainzer et al. 2014a).Jiang et al. (2020) investigate the size, thermal inertia, surface roughness, and geometric albedo of 10 Vesta family asteroids using the Advanced Thermophysical Model based on the thermal-infrared data acquired by WISE.Similarly, Jiang & Ji (2021) derived thermal inertia for 20 Themis asteroid family members from thermophysical modeling based on the WISE/NEOWISE observations.MacLennan & Emery (2021) used infrared data for a set of 239 objects observed by the WISE to estimate their size, albedo, thermal inertia, and surface roughness.Marciniak et al. (2021) used the light curves datasets in combination with the thermal infrared data, mainly from the Infrared Astronomical Satellite (IRAS, Neugebauer et al. 1984), AKARI (aka as ASTRO-F or IRIS -InfraRed Imaging Surveyor) (Usui et al. 2011), and WISE in an optimization process using the Convex Inversion Thermophysical Model.The authors analyzed the properties of 16 slowly rotating asteroids, including estimating their thermal inertia.Hung et al. (2022) used thermal flux measurements obtained by the WISE/NEOWISE during its fully cryogenic phase to obtain new thermophysical model fits of 1847 asteroids, deriving their thermal inertia, diameter, and Bond and geometric albedos.Jiang et al. (2023) studied the surface properties of asteroid (704) Interamnia, combining the mid-infrared wave bands observations from the Subaru telescope along with those obtained by IRAS, AKARI, and WISE/NEOWISE.Alí-Lagoa et al. (2020) used Herschel Photodetector Array Camera and Spectrometer data taken during the Asteroid Preparatory Programme for Herschel, ASTRO-F and the Atacama Large Millimeter/submillimeter Array (Müller et al. 2005) for a thermo-physical characterization of 18 large main belt asteroids.
We briefly reviewed only the most recent works as they provide the most currently available thermal inertia estimates.Nevertheless, we still have thermal inertia estimations for only a few hundred asteroids, primarily large main-belt asteroids, and some near-Earth asteroids (NEAs).Therefore, it is still an essential task, and the surface thermal inertia (TI) estimation for any new object is very important.Alternative methods for the TI estimation, using other kinds of observations, would be of great scientific value.
Recently, Fenucci et al. (2021) introduced a statistical method for the thermal estimation of NEAs based mainly on ground-based observations.The method was successfully applied to determine the thermal inertia of two rapidly rotating NEAs, asteroids 2011 PT and 2016 GE1 (Fenucci et al. 2021(Fenucci et al. , 2023)).The method requires, in fact, a determination of the semi-major axis drift produced by the Yarkovsky effect (see, e.g.Bottke et al. 2006;Vokrouhlický et al. 2015, for a review), a thermal effect caused by sunlight that strongly depends on the surface properties.Measurements of the Yarkovsky effect can be performed through orbit determination employing astrometric and/or radar observations (Chesley et al. 2003;Milani & Gronchi 2009), which are typically ground-based.Then, by matching the measured semi-major axis drift with that predicted by a physical model, the method produces a distribution of the surface thermal inertia Γ, possibly giving further constraints on other physical parameters.The Yarkovsky effect has been detected for hundreds of NEAs (Farnocchia et al. 2013;Del Vigna et al. 2018;Greenberg et al. 2020), and more and more detections are foreseen in the future because of two reasons: 1) the length of the observational arc can be only extended when NEAs are re-observed, and 2) new ground-based surveys like the Vera Rubin Observatory (Jones et al. 2018) and the ESA Flyeye telescope will provide more chances of re-observe known NEAs or recover the lost ones.The determination of the Yarkovsky effect is also starting to become a routine operation for the automatized orbit determination systems provided by the ESA NEO Coordination Centre1 (NEOCC), the NEODyS2 service, and the NASA JPL Small-Body Database3 .This situation allows us to apply the method by Fenucci et al. (2021) to many NEAs.An essential advantage of this model is that it can estimate the thermal inertia of comparatively smaller asteroids than the other methods.However, despite being successfully applied for rapidly rotating objects, the method proposed by Fenucci et al. (2021) still requires several improvements as well as additional validations before it can be safely applied to a large number of objects.
Motivated by these facts, one of the aims of the "Demystifying Near-Earth Asteroids" (D-NEAs) project funded by The Planetary Society4 is to develop further methods for the characterization of NEAs based on ground-based observations.In this paper, we propose extensions of the concept introduced in Fenucci et al. (2021Fenucci et al. ( , 2023) ) in order to make the thermal inertia estimation more reliable and accurate.In particular, we implemented different Yarkovsky effect models suitable for different dynamical and physical characteristics scenarios.Moreover, we refined the population-based modeling of poorly constrained or unknown physical parameters needed for the Yarkovsky modeling.We implemented all the methods and models described in the present paper in the Asteroid Thermal Inertia Analyzer (ASTERIA).The ASTERIA software v1.0.0 (Fenucci et al. 2023) is publicly available to the scientific community under the CC BY-NC-SA 4.0 license, and it can also be retrieved from a dedicated online repository5 .
The paper is structured as follows.In Sec. 2, we describe the general idea and Monte Carlo method for thermal inertia estimation.In Sec. 3, we give the details of the Yarkovsky effect models included in the ASTERIA software.Sec. 4 is dedicated to the explanations on modeling of the physical parameters entering the Yarkovsky effect model.In Sec. 5, we presented a basic introduction to the software, while Sec.6 is devoted to the model and code validation on Bennu and 10 other well-characterized NEAs.In Sec. 7, we have presented new thermal inertia estimates for 38 NEAs.Finally, we provide a summary and main conclusions in Sec. 8.

MONTE CARLO METHOD OF THERMAL INERTIA ESTIMATION
The method of thermal inertia estimation introduced by Fenucci et al. (2021) relies on the assumption that a measure (da/dt) m of the Yarkovsky effect obtained from astrometric data (see, e.g.Farnocchia et al. 2013;Del Vigna et al. 2018;Greenberg et al. 2020) is available for the object we are interested in.Models for Yarkovsky estimation by orbit determination (see, e.g.Milani & Gronchi 2009) can evaluate the value of the semi-major axis drift accurately, provided that the observational arc is long enough and the observations are of good quality.These models are typically independent of the physical characteristics of the object.In turn, theoretical models of the Yarkovsky effect rely on the knowledge of the physical properties of the asteroid (see, e.g.Bottke et al. 2006).If the Yarkovsky effect is independently determined from astrometry, the theoretical model of the Yarkovsky effect is constrained, i.e., its solution is known.Therefore, the theoretical model could be inverted, and at least one of its parameters could be estimated.
Yarkovsky effect models depend on several parameters, which are usually the semi-major axis a and the eccentricity e of the asteroid orbit, the diameter D, the density ρ, the thermal conductivity K, the heat capacity C, the obliquity γ, the rotation period P , the absorption coefficient α, and the emissivity ε.If all the mentioned parameters but the thermal conductivity K are fixed, then the model vs. observed Yarkovsky drift equation can be solved for K, thus determining the thermal conductivity value that allows the model da/dt( • ) on the lefthand side of Eq. (1) to reproduce the observed semi-major axis drift.The choice to determine the parameter K is dictated by the fact that it is the most uncertain one, contrary to the other parameters that can be either estimated or modeled with a fair degree of accuracy.In fact, the thermal conductivity K strongly depends on the characteristics of the material present on the surface of the asteroid, and it can range from ∼ 10 −5 W m −1 K −1 for very porous regolith material (Krause et al. 2011;Sakatani et al. 2017), up to ∼ 50 W m −1 K −1 for very conductive materials such as irons (Noyes et al. 2022, and references therein).The Monte Carlo (MC) method presented here assumes that all the parameters but K have some characteristic distribution, and then Eq. ( 1) is solved for K over a large set of combinations of input parameters.These parameters are randomly generated according to the assumed distributions.Solutions of Eq. ( 1) are searched for in a user-defined interval [K min , K max ], and they are computed with a bisection method with a tolerance of 10 −11 W m −1 K −1 .Moreover, usually, there is more than one K solution of Eq. (1) for a fixed combination of the remaining parameters.By collecting all the solutions obtained for given input combinations, the MC method produces a probability distribution function (PDF) of K in the output.A PDF for thermal inertia Γ is then obtained by converting each K solution using the relation and the corresponding values of ρ and C.
The level of accuracy of the estimation of Γ that the MC method described above can achieve depends, therefore, on two factors: 1) the accuracy of the Yarkovsky effect modeling, and 2) the reliability of the distributions assumed for the input parameters.In our software, we implemented different Yarkovsky effect models that are suitable for different dynamical and physical scenarios.The details of these models are described in Sec. 3.
The distributions considered for the orbital and physical parameters also depend on what properties of the studied asteroid are known.Because our method assumes the availability of a measurement of the Yarkovsky effect by astrometry, it implies that the orbital parameters are all determined with very good accuracy, with uncertainties in the semi-major axis typically of the order of ∼ 10 −9 au and of ∼ 10 −8 in eccentricity.Therefore, the orbital parameters are always assumed to be fixed at their nominal values, since their uncertainties produce negligible variations on the predicted Yarkovsky drift.On the other hand, physical properties may be estimated or constrained utilizing different kinds of observations.However, the number of NEAs for which physical properties are determined is still small, and models relying on the general properties of the whole NEA population are needed.In Sec. 4, we describe the models adopted for the input physical parameters.

YARKOVSKY EFFECT MODELS
The code provides the user with three different models for the computation of the Yarkovsky effect: 1) an analytical model for circular orbits, 2) a semi-analytical model for eccentric orbits, and 3) a semi-analytical model for 2-layered bodies.In addition, variable thermal inertia along the orbit is also implemented in models 2) and 3).We describe hereafter the details of these models.

Analytical model for circular orbits
The analytical model for the Yarkovsky effect was described in Vokrouhlický (1998Vokrouhlický ( , 1999)).The model assumed the asteroid to be a spherical body moving on a circular orbit around the Sun, that rotates around a fixed axis.The boundary conditions for the heat diffusion equation are linearized, and the drift in semi-major axis (da/dt) due to the Yarkovsky effect is given by two distinct components: the seasonal effect (da/dt) s , and the diurnal effect (da/dt) d .These effects have the following form: In the above equations, α is the surface absorption coefficient, Φ is the radiation pressure coefficient (e.g., Vokrouhlický et al. 2015), γ is the spin axis obliquity, and ω rev is the orbital frequency.R ′ s and R ′ d are the scaled (non-dimensional) values of the radius R, defined as where l s , l d are the penetration depths of the seasonal and diurnal thermal waves given by The penetration depths l s and l d depend on the thermal conductivity K, the heat capacity C, and the density ρ of the asteroid.Additionally, the two length scales depend on the respective frequencies: (i) the spin frequency ω rot = 2π/P in the case of the diurnal effect where P is the rotation period, and (ii) the orbital frequency ω rev in the case of the seasonal effect.The thermal parameters Θ s , Θ d also depend on the physical and thermal characteristics of the object, and they are defined as In the above equations, σ is the Stefan-Boltzmann constant, ε is the emissivity and T ⋆ is the subsolar temperature, defined by εσT 4 ⋆ = αE ⋆ , with E ⋆ being the solar radiation flux at the distance of the body.The function F in Eqs. ( 3) and ( 4) is given by It depends on both the corresponding scaled radius and the thermal parameter.The coefficients k 1 , k 2 , k 3 are positive analytical functions of the scaled radius, and their complete definition can be found in Vokrouhlický (1998Vokrouhlický ( , 1999)).
Figure 1 shows these coefficients as a function of R ′ .For R ′ large enough, all the coefficients approach the value 1/2.It is worth noting that the seasonal component always produces an inward migration, while the direction of migration for the diurnal component depends on the obliquity γ, which is maximized for γ = 0 • , 180 • .It is also important to note here that a double-peaked distribution is often seen in the thermal inertia estimates (see, e.g., Fig. 3), and it is a consequence of Eq. ( 8).In most cases, one of the peaks could be discarded either based on the number of possible solutions or rejected due to unphysical values of thermal inertia (e.g., Γ > 2000 J m −2 K −1 s −1/2 ).

Semi-analytical model
When an asteroid's orbital eccentricity is large, the model described in Sec.3.1 may not be accurate for estimating the Yarkovsky effect.In this case, we can compute the instantaneous semi-major axis drift caused by the Yarkovsky effect as where a is the semi-major axis of the orbit, n is the mean motion, v is the heliocentric orbital velocity, and f Y is the instantaneous value of the Yarkovsky acceleration.The term f Y is computed with analytical formulas by Vokrouhlický et al. (2017), again assuming a homogeneous spherical shape of the asteroid and a linearization of the surface boundary condition.This acceleration is given by where f Y, d , f Y, s are the diurnal and the seasonal component, respectively.The diurnal component is expressed as In Eq. ( 11), n = r/r is the heliocentric unit position vector, and s is the unit vector of the asteroid spin axis.In addition, where S = πR 2 is the cross-section of the asteroid, E ⋆ is the solar radiation flux at a heliocentric distance r, m is the asteroid mass, and c is the speed of light.The coefficients γ 1 , γ 2 are expressed as where Θ d is defined by Eq. ( 7) and k 1 , k 2 , k 3 are the same coefficients as in Sec.3.1.The seasonal component is given by where N is the unit vector normal to the orbital plane, and γ1 , γ2 have the same expressions as Eq. ( 13), but evaluated with the scaled radius R ′ s of Eq. ( 5) and thermal parameter Θ s of Eq. ( 7).The average Yarkovsky drift da/dt is then obtained by averaging the instantaneous Yarkovsky drift of Eq. ( 9) over an orbital period as where ℓ is the mean anomaly.The integral at the right-hand side of Eq. ( 15) is computed with the trapezoid rule (see, e.g.Stoer & Bulirsch 2002).To prevent numerical errors at high eccentricity, we computed the integral of Eq. ( 15) by using a fixed step in the eccentric anomaly u, which provides a better sampling of the orbit near the perihelion.

Semi-analytical 2-layer model
In-situ observations made by spacecrafts in the last two decades showed that asteroids have many different surface structures.Some of them, like Bennu or Ryugu, show very rough surfaces composed of boulders spanning orders of magnitude in diameter, basically without any area covered by fine regolith (Lauretta et al. 2019;Watanabe et al. 2019;Sugita et al. 2020).Other objects like Itokawa or Eros present large smooth areas covered with fine dust and regolith (Cheng 2002;Miyamoto et al. 2007), similar to the surface of the Moon.At the same time, their interior may be constituted by much larger boulders, thus breaking the homogeneity assumption made in the Yarkovsky models described in Sec.3.1 and 3.2.This last case is particularly important for correctly estimating the Yarkovsky effect because regolith layers are typically composed of loose agglomerates of small rocky particles with large porosity, decreasing the thermal conductivity K with respect to typical values of bare rocks.Since the thermal parameter K can largely change the magnitude of the Yarkovsky effect, the 2-layer model is essential to model regolith-covered asteroids correctly.
In Vokrouhlický & Brož (1999), the authors considered a spherical body composed of a core with parameters ρ 2 , C 2 , K 2 , covered by a shell of thickness h with parameters ρ 1 , C 1 , K 1 .In order to find a solution to the heat conductivity equation, linearized boundary conditions on the outer shell were considered.In addition, a continuity condition on the temperature distribution T and the heat flux when passing between the layers was also assumed.Under these assumptions, the authors found an analytical solution for the seasonal Yarkovsky effect and discussed the difficulties in finding a solution for the diurnal effect for meter-sized objects.The problem of finding analytical solutions to the heat diffusion equation on non-homogeneous bodies was also discussed in Čapek & Vokrouhlický (2012).A complete analytical model for the Yarkovsky effect on 2-layered asteroids is unavailable.
Due to the above difficulties, we used a simplified approach to define a 2-layer Yarkovsky model.We assumed the surface material to have a density of ρ surf , while the interior to have a density ρ.To define the thickness of the top layer, we assumed that the penetration depths of the seasonal and diurnal thermal waves are smaller than the real regolith thickness.At least roughly speaking, in this case, the magnitude of the Yarkovsky effect depends only on the heat transport in the top layer.
The second layer represents the body interior, excluding only the small top layer.It is assumed, however, that the second (interior) layer does not affect the temperature distribution at the surface.It is taken into account only when asteroid mass is computed.
Therefore, the thickness of the top layer could be approximated with the length of the more relevant penetration depth, which is the one related to the more important component of the Yarkovsky effect.In most cases, this is the diurnal penetration depth l d , except when an asteroid's obliquity is close to 90 degrees when it is assumed that the top layer corresponds to the seasonal penetration depth l s .
The two penetration depths are obtained as Once the top layer has been constrained, we used its adopted parameters to compute the Yarkovsky effect with Eq. ( 9).
We recall that this Yarkovsky model is valid only as long as the penetration depths of Eq. ( 16) are smaller than the actual regolith thickness, which is typically satisfied when R ′ s and R ′ d are >> 1.Therefore, it is the user's responsibility to make use of this model with more caution by checking that the values of l s and l d obtained in output are reasonably smaller than the presumed thickness of the surface material.Generally, the approach used for the 2-layer model is an approximation.In particular, we found that it generally tends to underestimate the dust layer's role.Therefore, in its current form, the 2-layer model may only indicate how significantly a possible presence of the surface dust layer may affect the thermal inertia estimation by the ASTERIA.

Variable thermal inertia
Thermal inertia depends on the thermal conductivity and heat capacity through Eq. ( 2).Laboratory analyses on lunar regolith samples (Keihm 1984) and meteorites (Opeil et al. 2020) showed that both thermal conductivity K, and heat capacity C depend on the temperature.Therefore, the thermal inertia also depends on the temperature, which is, in turn, related to the heliocentric distance r.Accordingly, accounting for its variation along the asteroid orbit is important to estimate the thermal inertia correctly.
For predominantly radiative heat transfer, as expected in a fine regolith on an airless body, the thermal inertia scales with Γ = Γ 0 r −3/4 , (17 where Γ 0 is the value of the thermal inertia at 1 au.In a more recent study, Rozitis et al. (2018) found that the trend of Eq. ( 17) found by theoretical models does not always agree with the observations, but more extreme values of the exponent may be possible.Similar conclusions were also drawn by MacLennan & Emery (2021).
To account for the variable thermal inertia in our model, yet keeping it flexible concerning the exact dependence on the heliocentric distance, we assume that the thermal inertia Γ varies with the heliocentric distance as The exponent β is left as a free parameter to be chosen by the user.The variable thermal inertia is implemented for the Yarkovsky models described in Sec.3.2 and 3.3.On the purely technical side, we underline that when the variable thermal inertia model is used, output from our model is thermal inertia Γ 0 , which is a value normalized to r = 1 au.

Non-sphericity effects
The above-described Yarkovsky models assume a spherical shape of the body.However, it is known that asteroids are not perfectly spherical but rather elongated, in some cases even very elongated bodies (see, e.g.McNeill et al. 2019, and references therein).Vokrouhlický (1998) found that in the general case, the diurnal Yarkovsky component for spheroids and a sphere does not differ more than 30%.However, in some cases, the difference could be up to a factor of 2 or 3.Such an example is fixed spin axis orientation at high obliquity.
As the Yarkovsky detection is performed from relatively short time intervals (up to a few tens of years), it is reasonable to consider the orientation of the asteroid spin axis fixed in space.
Assuming a biaxial ellipsoid shape (spheroid), rotation around the shortest axis, and high obliquity, an elongated object has a smaller average cross-section than a spherical body of the same equivalent diameter.This implies that elongated bodies are subject to a smaller Yarkovsky effect than spherical ones of the same equivalent diameter.Therefore, to account for the impact of the body's non-sphericity, we introduced a correction factor ξ into the model.We found that the non-sphericity correction factor could be reasonably well approximated as ξ = f −0.3 .The factor is defined to mimic the results presented in Vokrouhlický (1998).The value f is the ratio of the a/b axes and is computed from the maximum light curve amplitude ∆m using the relation ∆m = 2.5 log f .As we assume a > b > 0, the axis ratio f and the correction factor ξ are always in the (0, 1) interval.
As the non-sphericity under the given assumptions decreases the Yarkovsky effect, and the correction factor expressed the magnitude of the decline, to compensate for the non-sphericity, we divided the measured drift in the semi-major axis da/dt with the factor ξ. In other words, we artificially increased the measured Yarkovsky drift to match the one expected for a spherical object having the same equivalent diameter.We point out to use this scaling with caution as it might break down for extreme light curve amplitudes ∆m 2.0.

Albedo distribution
The albedo p V is typically determined from the infrared observations.It is nowadays available for a large number of asteroids thanks to space surveys dedicated to this purpose, such as the IRAS (Tedesco et al. 2002), the WISE (Masiero et al. 2011;Nugent et al. 2015), and the AKARI mission (Alí-Lagoa et al. 2018).
When albedo and its associated uncertainty are known, the user can select an option that uses the nominal value and assumes a Gaussian distribution for the measurement error.
If p V is not known, we model this parameter by combining the NEA absolute magnitude and orbital distribution by Granvik et al. (2018), and the NEA albedo distribution by Morbidelli et al. (2020).A first description of this model for the NEA 2011 PT was given in Fenucci et al. (2021).Given the orbital elements (a, e, i) of the NEA and its absolute magnitude H, the model by Granvik et al. (2018) provides 7 probabilities P s (a, e, i, H), s = 1, . . ., 7. Each of these numbers corresponds to the probability that the NEA arrived in the NEA region through one of the 7 seven transport routes from the main belt: 1) the ν 6 secular resonance; 2) the 3:1 Jupiter mean-motion resonance (JMM); 3) the 5:2 JMM; 4) the Hungaria region; 5) the Phocaea region; 6) the 2:1 JMM; 7) the Jupiter Family Comets (JFC) region.In the NEA albedo distribution by Morbidelli et al. (2020), asteroids are divided into three albedo categories: c 1 corresponds to dark objects with p V ≤ 0.1, c 2 corresponds to moderately bright objects with albedo 0.1 < p V ≤ 0.3, and c 3 corresponds to bright objects with p V > 0.3.For each source region, indexed with s, the model provides the fraction of objects p s (c i ), i = 1, 2, 3 belonging to the category c i that arrived from the source region s.Numerical values are given in Table 1.To define a PDF for the albedo of a specific object, we first define an albedo PDF p s (p V ) for each escape route s.This PDF is uniform in the categories c 1 and c 2 , and exponentially decaying in c 3 (Morbidelli et al. 2020) as The albedo PDF p(p V ) is then defined as

Diameter and density distributions
The input diameter and density distributions could be either population-based or user-provided.The diameter of an asteroid can be determined with different methods, including radar observations (Ostro 1985), thermal infrared modeling (Harris 1998), and stellar occultation (Millis & Elliot 1979).If the asteroid's diameter is measured, the user can specify the mean value and the 1-σ uncertainty to model D with a Gaussian distribution.
On the other hand, if the diameter is not available, it could be estimated by the conversion formula (see, e.g.Bowell et al. 1989;Pravec & Harris 2007) The density is even more challenging to measure.Typical methods to determine the mass of an asteroid are through close encounters with smaller bodies, spacecraft tracking, and the orbit of a satellite (Carry 2012).Then, the density ρ is obtained by combining the estimate for m and that of the diameter D. Therefore, direct density determination is available for a small number of asteroids.Nevertheless, to keep the code flexible, we allow the user to choose a normal distribution for the density ρ, with user-specified mean and standard deviation values.
Another indirect way to have an estimate of the density is through the spectroscopic classification (Tholen 1984;DeMeo et al. 2009;Carvano et al. 2010).Indeed, asteroids with different spectra are supposed to have different compositions.The density of each taxonomic class can be assumed as the same of other asteroids with the same spectral class for which the density was determined with one of the methods mentioned above (see Carry 2012, for a review about asteroid densities), or by the association of the spectral type with known meteorites (see, e.g.Binzel et al. 2019;DeMeo et al. 2022).
When estimates of D and ρ are unavailable, we model these parameters with population-based properties.The model also relies on the fact that diameter D and the density ρ are not uncorrelated because objects with large albedo correspond to smaller asteroids with relatively large density, and objects with low albedo correspond to larger asteroids with relatively low density.Therefore, we generate a joined (D, ρ) distribution using geometric albedo p V obtained as explained in Section 4.1.A sample of p V is obtained first, and then each sample value is converted into D and ρ.The diameter D is obtained through Eq. ( 21), where H is assumed to be Gaussian distributed.The density ρ is instead generated according to the category c j into which p V falls.We associate category c 1 to the C-complex of carbonaceous asteroids, category c 2 to stony asteroids in the S-complex, and category c 3 to bright asteroids of the X-complex.Then, we randomly generate a value of ρ depending on the complex associated with the value of p V .The densities of the complexes are assumed to be distributed according to a log-normal distribution, with parameters listed in Table 2.
Table 2. Parameters of the log-normal distributions of the density for the three different complexes of asteroids.

Surface density distribution
The two-layer model of Sec.3.3 also needs the distribution of the density of the surface material ρ surf in input.As discussed above, this model is supposed to be used only in the case of a body covered by a regolith layer.Even though regolith comprises a large variety of dust and gravel sizes, it generally refers to an unconsolidated, loose, porous, and heterogeneous set of small rocky materials.
Thanks to in-situ observations by spacecrafts and to the samples returned from the Moon surface during the Apollo missions, the regolith density has been constrained by direct measurements (Mitchell et al. 1972).The average density of the lunar surface regolith is about 1300 kg m −3 .Recent works based on the lunar far side obtained values for the upper top level as low as about 500 kg m −3 (Xiao et al. 2022).We allow the user to use only a Gaussian distribution for ρ surf , by specifying the mean and standard deviation values.The critical point is that the two-layer model makes sense only if the surface density is significantly lower than the bulk density.Therefore, values above 1300 kg m −3 should be avoided.For example, we consider the values of 1100±100 kg m −3 , as a good compromise.Nevertheless, the users are welcome to test other values.

Obliquity distribution
The obliquity of an asteroid is typically determined by light curve inversion method that makes use of photometric observations (Kaasalainen & Torppa 2001;Kaasalainen et al. 2001Kaasalainen et al. , 2002;;Durech et al. 2015), or by radar observations (Hudson 1994;Hudson et al. 2000;Magri et al. 2007).When determinations are available, the user can decide to use a Gaussian distribution with a mean value equal to the nominal estimated obliquity, and standard deviation equal to the 1-σ uncertainty.
When the obliquity is unknown, we rely again on the global properties of the NEAs population, and we model this parameter with the NEA obliquity distribution found by Tardioli et al. (2017).The distribution of γ is best fitted by a quadratic PDF in cos γ of the form where a = 1.12, b = −0.32,and c = 0.13.The resulting probability density function of the NEAs' obliquities derived from Eq. ( 22) is shown in Fig. 2. The adopted distribution has about a 2:1 ratio between retrograde and prograde rotators (La Spina et al. 2004) and represents the recent obliquity model of near-Earth objects population (Tardioli et al. 2017).In practice, this should not affect the final result produced by the ASTERIA model.The values of γ not compatible with the measured Yarkovsky drift are always rejected.

Rotation period distribution
The rotation period P is an essential parameter for correctly estimating the Yarkovsky effect because it may change the magnitude by a factor of several.It is typically evaluated by analyzing the light curve, and most of the measurements can be found in the online Asteroid Light curve Database6 (LCDB; Warner et al. 2009).
While the rotation period of NEAs generally shows a trend with the diameter D, modeling P in different size categories based on NEA population properties would produce poor results in estimating the thermal inertia Γ with the MC method.This is because the distribution of P has a large standard deviation, meaning that a population-based distribution is not very representative of the true rotational state the asteroid is in.Additonally, Pál et al. (2020) found that the number of bodies with long rotation periods is underestimated by the ground-based surveys, implying that the LCDB catalogue could be biased.For these reasons, we decided not to model P with population-based information, and we limited ourselves to providing the possibility of using a Gaussian distribution model with userspecified parameters.Therefore, the period needs to be known, or the model should be run for at least several discrete values.

Heat capacity
The heat capacity value C depends on the physical characteristics of the materials composing the asteroid.It can vary within a factor of several (Delbo' et al. 2015).Typical values assumed for rocky and regolith-covered mainbelt asteroids range from 600 to 750 J kg −1 K −1 , while they go as low as 500 J kg −1 K −1 for iron-rich asteroids (Farinella et al. 1998).In addition, the heat capacity shows to increase with increasing temperatures, and therefore, it could be higher for NEAs, though not exceeding ∼1200 J kg −1 K −1 (e.g.Opeil et al. 2020).However, there is insufficient data to extrapolate a general distribution valid for the whole NEA population.Therefore, the model keeps the heat capacity C fixed, and the user can specify its value.We anyway suggest keeping C within the reasonable ranges discussed above.We also note that the thermal inertia estimations by the ASTERIA model are not sensitive to changes in C.However, the code also provides an opportunity to estimate the thermal conductivity.In the latter case, caution is needed when selecting the input heat capacity, and we suggest running the code with several different values.

Absorption coefficient and emissivity
The absorption coefficient is defined as α = 1 − A, where A is the Bond albedo.The Bond albedo A is related to the geometric albedo p In the above formula, q is the phase integral defined by where G is the slope parameter.If the user provides a measurement of the geometric albedo p V , Eq. ( 23) is used to produce a distribution of α, which therefore follows Gaussian distribution.On the other hand, the population-based albedo distribution of Eq. ( 20) is used to produce a distribution of α.In all the computations, the phase angle is assumed to be G = 0.15.The emissivity ε can be constrained by measurements performed on meteorite samples.Ostrowski & Bryson (2019) collected the measurements of ε of 61 meteorites, and all the objects but one have a value ranging between 0.9 and 1, with an average value of 0.984.Since the interval in which the emissivity may change is small, we expect changes in the estimated Yarkovsky effect to be small, and therefore we decided to keep this parameter fixed in the MC method.The value can be specified by the user, but if there is no other reason to believe it to be smaller than 0.9, we suggest using the average value ε = 0.984 found for meteorite samples.

Measured Yarkovsky effect distribution
As mentioned in Sec. 2, the availability of a measurement of the Yarkovsky effect obtained by astrometry is a fundamental requirement for our MC method.Therefore, we always assume the semi-major axis drift (da/dt) m to be Gaussian distributed with a mean value corresponding to the nominal value, and standard deviation equal to the 1-σ uncertainty.These values can be specified by the user in the input of the software.

Open-source software and repository
All the methods and models described in Sec. 2, 3, and 4 are implemented in the ASTERIA software v1.0.0 (Fenucci et al. 2023), which is written in fortan90 programming language.A publicly available repository is also maintained on GitHub7 .The ASTERIA software is released under the Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) license.Automatic scripts, working under UNIX systems, for the compilation of the source code with the GNU gfortran or Intel ifort compilers, are included in the distribution.
Once compiled, the ASTERIA software provides the user with two executable drivers: 1) the gen distrib.xdriver, implementing the methods described in Sec. 4, and 2) the gamma est mc.x driver, implementing the MC method described in Sec. 2. The driver gen distrib.xneeds to be executed first in order to produce the input distributions of the physical parameters entering the Yarkovsky modeling.The files generated are then used by the gamma est mc.x driver for the MC method of thermal inertia estimation.This separation was done to keep the code more flexible.The user is strongly encouraged to use the gen distrib.xdriver for generating the input distributions, but one can also produce their own files with parameters' distributions and use them as input for the gamma est mc.x driver.The driver gamma est mc.x for the MC estimation of thermal inertia has been parallelized by using the OpenMP API, and it is, therefore, able to run on multi-CPUs machines.Detailed instructions for the compilation and for the usage of these drivers can be found in the Software User Manual included in the ASTERIA software repository.

MODEL VALIDATION
Model validation is a key procedure for quantifying the model's accuracy, i.e., determining how well the model describes the actual conditions.Validation aims to quantify confidence in the model's predictive capability by comparison with available data.
To validate our model, we selected a set of eleven asteroids with thermal inertia estimates available in the literature as determined by means of other models and techniques.The group essentially consists of the asteroid Bennu plus ten other near-Earth asteroids.The asteroid (101955) Bennu, a target of NASA's OSIRIS-REx mission, is the most suitable object for this purpose.It has practically all parameters relevant to the ASTERIA model very well constrained (Table 3), including the thermal inertia, which is the output of the model.Additional objects selected for our validation purposes are those NEAs well-characterized by ground-and space-based telescopes.Therefore, we first validated our model on Bennu as the most reliable reference8 , then extended it to ten other NEAs.The validation on Bennu aims to show that our model can reproduce the thermal inertia estimate when the input parameters are well-known.In contrast, the test on the other objects primarily seeks to show that the model gives reliable results in more realistic situations when only some input parameters are known and that the results are consistent with the available estimates in the literature.

Model validation on asteroid Bennu
A preliminary model and ASTERIA code testing on asteroid Bennu have already been performed by Fenucci et al. (2023).This work investigated how various unknown input parameters affect the final thermal inertia estimates, and showed that the model can reproduce reasonably well Bennu's thermal inertia.Building on these findings, here we extended this analysis, investigating the other relevant factors that may affect the match between our estimate and the reference value(s).In particular, the factors that were not considered by Fenucci et al. (2023) are the role of the non-spherical shape and of the spatial variations of Bennu's surface thermal inertia.Referent thermal inertia estimates for asteroid Bennu, used to benchmark our model, are listed in Table 4.
For the input parameters, we did not use population-based results but generated corresponding distributions based on the values listed in Table 3.In addition to the parameters listed in the table, we adopted a heat capacity value9 of C = 750 J kg −1 K −1 , following Rozitis et al. (2020).
Also, as there are different variations and options in our model, as described in Section 3, we analyzed which one should be the most appropriate for asteroid Bennu, and set up our nominal model settings.In particular, we implemented into the ASTERIA two Yarkovsky models based on circular and eccentric orbit assumptions.In this respect, we underline that the main advantage of the circular model is that it is fast and provides results almost instantaneously.In contrast, the eccentric-orbit-based model is a much more time-consuming procedure.With an eccentricity of about 0.2, the orbit of Bennu is not very eccentric, and even a circular model might be appropriate.However, aiming to obtain as accurate results as possible, we used solely the eccentric Yarkovsky model for Bennu.
Similarly, we decided not to use a two-layer density model, as no fine regolith (dust) has been found at the surface of Bennu (Dellagiustina et al. 2019).Rozitis et al. (2020) found that thermal inertia variations with temperature (hence with heliocentric distance) at Bennu are small.For that reason, we used constant thermal inertia in our nominal model.
In addition, we considered other aspects that could affect the results in the case of asteroid Bennu.One of these is the correction due to the non-sphericity of Bennu.As the shape of the Bennu is well-defined, we used the measured axis ratio here to compute the non-sphericity factor, instead of obtaining it from the light curve amplitude. 10The a/c ratio of 1.1 yields the correction factor of ξ = 0.95.
Furthermore, as discussed by Fenucci et al. (2021), nonlinearity effects decrease the theoretically predicted semimajor axis drift rate da/dt, with a reduction factor ranging between 0.7 and 0.9.Conversely, thermal beaming effects always increase the semi-major axis drift rate by a factor ranging between 1.1 and 1.5.These two effects tend to compensate for each other, although some residual factors may remain in particular cases.Although we cannot properly resolve these two issues, we attempt to reduce their importance by artificially increasing the uncertainty of the measured semi-major axis drift rate.In this case, we assumed that the uncertainty was 5% of the nominal value, reducing the corresponding detection level to 20-σ, significantly lower than the 1400-σ obtained by Farnocchia et al. (2021).
Finally, for the emissivity, we used a value of 0.95, as in Rozitis et al. (2020).With this in mind, we defined the above-described model settings as nominal.
Running the model with the nominal settings, the resulting median value of the thermal inertia and its corresponding 1-σ lower and upper uncertainties for the higher peak are Γ=381 +55 −56 J m −2 K −1 s −1/2 .The obtained value is generally compatible with all the literature estimates listed in Table 4.Still, we note that our value is shifted towards higher values with respect to other estimates.What could be the reasons for that?The ASTERIA model uses the orbit-averaged Yarkovsky effect to, in principle, estimate the global average thermal inertia.However, the model might, to some degree, depend on the spatial variations of the surface thermal inertia.As such variations have been observed on Bennu (Rozitis et al. 2020), we investigate how they could affect the results.
Rozitis et al. found larger thermal inertia at Bennu's equator than at polar regions.Assuming spatial variations in thermal inertia, the quantity relevant to Yarkovsky is proportional to Γ(ϕ)cos 2 (ϕ)dS, where Γ(ϕ) is the thermal inertia at a given latitude ϕ, and dS is the surface element.The cosine is squared because both the size of the surface element and the projection of its normal onto the orbital plane are proportional to the cosine.
Accordingly, a "global" averaged thermal inertia relevant to the Yarkovsky effect (T I r ) could be obtained from the following equation: If we assume linear variation of TI with latitude, then the ratio of T I r and average thermal inertia T I a results in where T I e and T I p are thermal inertia at the equator and pole, respectively.In the case of Bennu, based on the results by Rozitis et al., this translates into a small decrease of our thermal inertia estimate of about 1.2%, or maximum 2%, given that spatial variation in thermal inertia at Bennu is not linear from the equator to the pole.Considering this effect, our best estimate of the Bennu surface thermal inertia is 374 +54 −53 J m −2 K −1 s −1/2 .The corresponding thermal inertia distribution is shown in Fig. 3. Despite being shifted towards somewhat larger values, the result is fully compatible with the literature estimates listed in Table 4, with all terms within 1-σ of each other.Finally, though Rozitis et al. (2020) found only small thermal inertia variations with heliocentric distance at Bennu, and we did not consider it in our nominal model, we have investigated how applying our variable thermal inertia model would change the result.Therefore, assuming in the ASTERIA model that thermal inertia scales with heliocentric distance as r −0.75 , we found a value of Γ = 366 +53 −52 J m −2 K −1 s −1/2 .This result is even closer to the reference value of Bennu's thermal inertia, which might suggest that variable thermal inertia plays some role for Bennu.Nevertheless, as we do not know how representative the thermal inertia variation exponent of −0.75 is in the case of Bennu, we did not consider this result our best match.

Model testing -further verification on the Near-Earth asteroids
To estimate thermal inertia with the ASTERIA model, several input parameters should be provided.Along with absolute magnitude and orbital parameters, obviously, all the objects must have a Yarkovsky detection.11The distributions of the other parameters could be, in principle, obtained from the population models.However, as noted by Fenucci et al. (2023), if no other parameter is known, the model is not able to provide reliable results, except in special cases, such are super-fast rotators.For this reason, we require here that additionally, the size, albedo, and rotation period of the objects are known.In this case, the density distribution is determined as described in Section 4, but based on the determined albedo values and Gaussian distribution instead of the albedo obtained using the population model.Therefore, the level of the model reliability analyzed here applies to such cases.For objects where the other parameters are available as well, such as bulk density or obliquity, these are neglected for our model validation purposes.
Tables 5 and 6 list the selected NEAs for further verification of the ASTERIA model.The objects are selected so that, on one side, they have reasonably well-determined thermal inertia from TPM, and on the other side, to cover as much as possible variety of orbital and physical properties of NEAs.The Tables 5 and 6 also provide all the relevant orbital and physical parameters for the model validation set of asteroids.
The ASTERIA asteroid thermal inertia estimator gives values that are generally in good agreement with those derived from thermophysical modeling (see Fig. 4 and Table 6).As discussed below for individual cases, all the results obtained with our model are statistically compatible with the corresponding reference literature values.
Below, we briefly discuss the matching on a case-by-case basis.
(1620) Geographos: This object has been the subject of thermophysical modelling by Rozitis & Green (2014), and has trustworthy determined thermal inertia.Therefore, we include this object in the validation list.On the other hand, we note that its Yarkovsky detection has S/N of only ≈2.9, which is typically considered as marginal detection (e.g.Del Vigna et al. 2018).According to Rozitis & Green (2014), the thermal inertia of Geographos is 340 +140 −100 J m −2 K −1 s −1/2 .Our model with nominal Yarkovsky drift from the JPL suggests a higher thermal inertia of ∼500 J m −2 K −1 s −1/2 , with large uncertainties, primarily due to large uncertainty in Yarkovsky drift rate estimation.Nevertheless, the result obtained by Rozitis & Green (2014) is still within the 1-σ error limits of our estimate, and we consider two results to be roughly consistent.A bit larger S/N of ∼ 3.6 of the Yarkovsky detection for Geographos has been obtained by Dziadura et al. (2022), suggesting slightly faster Yarkovsky drift of da/dt = −1.34± 0.37 [10 −4 au Myr −1 ].Plugging this value into the ASTERIA model, we found Γ = 445 +205 −121 , which is fully in agreement with estimate from Rozitis & Green.(1685) Toro: This is the largest asteroid in our sample, with a diameter of about 3.8 km, but with solid detection of the Yarkovsky effect (S/N≈6.6).An estimated thermal inertia of 210 +85 −55 J m −2 K −1 s −1/2 has been recently given by Hung et al. (2022), which is in perfect agreement with our result of 193 +57 −53 J m −2 K −1 s −1/2 .(1862) Apollo: Its thermal inertia has been estimated to be 140 +140 −100 J m −2 K −1 s −1/2 by Rozitis et al. (2013).Our nominal result of 81 +50 −29 J m −2 K −1 s −1/2 is somewhat lower than the one derived by Rozitis et al.. Adopting the physical parameters (geometric albedo p v = 0.20 ± 0.02, and diameter D = 1.55 ± 0.07 km) obtained along with the thermal inertia by Rozitis et al. (2013), our model gives slightly higher thermal inertia of 86 +51 −30 J m −2 K −1 s −1/2 , indicating that the difference in the input parameters is not the main source of the observed difference in the results.Nevertheless, given the significant uncertainties of both results, we consider the two results consistent.We should underline, however, that for this object, the ASTERIA model also provides the alternative solution of higher thermal inertia in the range of 400−650 J m −2 K −1 s −1/2 , which we cannot reject based on the available data.Based on the Yarkovsky effect detection algorithm and interdependence between the bulk density and thermal inertia, Farnocchia et al. (2013) suggested that the TI of (1862) Apollo should be in the range 400−1000 J m −2 K −1 s −1/2 .Therefore, asteroid (1862) Apollo might have higher TI, in line with our alternative solution.
(1865) Cerberus: The Yarkovsky detection at JPL has S/N of only about 2.6, translating into large uncertainties of our thermal inertia estimate.Nevertheless, our value of 998 +625 −321 J m −2 K −1 s −1/2 agrees with the value of 809 +219 −134 J m −2 K −1 s −1/2 , found by Hung et al. (2022).We underline that in the case of the low Yarkovsky detection S/N, the results should be taken with caution.For instance, Yarkovsky detection for (1865) Cerberus has also been provided by the Greenberg et al. (2020).These authors found the semi-major axis drift of −3.75 ± 1.8 × 10 −4 au Myr −1 , therefore with even lower S/N od 2.1.With the Yarkovsky detection from Greenberg et al., our model gives the thermal inertia of 1149 +961 −426 J m −2 K −1 s −1/2 , formally in agreement with the estimate from Hung et al.. (25143) Itokawa: This asteroid has one of the largest known thermal inertia values of 700 ± 200 J m −2 K −1 s −1/2 (Müller et al. 2014a).Though with large uncertainties caused primarily by a small S/N of 3.3 in the Yarkovsky detection, our nominal result of about 686 J m −2 K −1 s −1/2 perfectly matches the literature value.
Table 6.The physical properties of 10 NEAs used to validate the ASTERIA model.In the last two columns, reference literature values of thermal inertia (ΓTPM) and our estimates (ΓAST) are given.In each case, our results refer to the median values and their corresponding 1-σ lower and upper spreads.Notes: The values of absolute magnitudes H are from the JPL, diameters D and albedos pV are from Mainzer et al. (2019), except for (1620) Geographos (Rozitis & Green 2014), ( 25143) Itokawa (Fujiwara et al. 2006;Tatsumi et al. 2018), (29075) 1950 DA (Masiero et al. 2021), ( 99942) Apophis (Brozović et al. 2018), and ( 161989) Cacus (Masiero et al. 2020).The values of maximum light curve amplitudes ∆m and rotation periods P are taken from LCDB (Warner et al. 2009).The reference literature values of thermal inertia ΓT P M are, in order of appearance in the table from top to bottom, from ) 1998 WT24: Based on the thermal-infrared observations and thermophysical modeling Harris et al. (2007) found thermal inertia of 200 ± 100 J m −2 K −1 s −1/2 , that our nominal result of 142 +83 −40 J m −2 K −1 s −1/2 matches very well.
(99942) Apophis: For this object, we found four thermal inertia estimates (Müller et al. 2014b;Licandro et al. 2016;Yu et al. 2017;Satpathy et al. 2022), which are not always consistent or are subject to large uncertainties.Nevertheless, we include it in our validation list.Generally speaking, our estimate of 491 +263 −169 J m −2 K −1 s −1/2 is consistent with three of the values from the literature, with the only exception being a low value of 100 +100 −52 J m −2 K −1 s −1/2 obtained by Yu et al. (2017).
(161989) Cacus: Thermal inertia of 548 +491 −169 J m −2 K −1 s −1/2 estimated by the ASTERIA model match very well the literature value of 650 ± 150 J m −2 K −1 s −1/2 (Hung et al. 2022).Despite a good S/N value of 5.7 in the Yarkovsky detection, our result comes with significant uncertainties caused by uncertainties in the size and albedo of the asteroid Cacus.
(175706) 1996 FG3: There are two available thermal inertia estimates for this object.Wolters et al. ( 2011) obtained 120 ± 50 J m −2 K −1 s −1/2 , while Yu et al. (2014) found a bit lower value of 80 ± 40 J m −2 K −1 s −1/2 .Our nominal result agrees with those estimates, with the mean value being somewhat smaller and in better agreement with the one from Yu et al. (2014).Using the size and albedo from Wolters et al. (2011), our model gives the thermal inertia of 98 +59 −37 J m −2 K −1 s −1/2 , which is in full agreement with the value found by those authors.The above-presented validation results showed that the ASTERIA model is reliable and provides consistent estimates of the thermal inertia with the values available in the literature.The prerequisite is that the diameter and albedo estimates are available for a given object along with the Yarkovsky effect detection.The accuracy of the results obviously depends on the quality of the input parameters, with the model being more sensitive to the Yarkovsky rate and diameter estimations than to the albedo.Nevertheless, as long as the same (or similar values) of the input parameters have been used, the ASTERIA model matches very well the other thermal inertia estimates, and it can be considered a good alternative to thermophysical modeling in proper circumstances.Notes: In the "No. of exposure" column, the numbers in brackets show how many images from each night were used for the period determination.
6.3.ASTERIA model application and criticality of input parameters -an example of asteroid 152671 (1998 HL3) After verifying our model in the previous sub-sections using the well-studied asteroid Bennu, and ten additional well-studied NEAs, here we applied the model to the near-Earth asteroid (152671) 1998 HL3 (hereafter HL3) as a case study of a more general situation.
Contrary to the case of Bennu, we know relatively little regarding the physical parameters of the asteroid HL3.Its known orbital and physical properties relevant to the model, including those determined in this work, are summarised in Table 7.
We also used this subsection to highlight how important is the determination of rotation periods for our model.In this respect, we present our photometric observation of asteroid HL3, aiming to obtain its light curve and derive its amplitude and rotation period.

Observations from Astronomical Station Vidojevica
Observations were collected on the 5th of May 2022, as well as on the 19th and 20th of May 2022 from the Astronomical station Vidojevica (MPC code C89), using the 1.4 m Milanković telescope.We used an Andor iKon-L 2024×2024 pixel CCD camera, which has a field of view of 13.3×13.3arcmin and a pixel size of 13.5×13.5µm.All observations were made in the Standard Johnson-Cousin R-filter.Light curve construction and rotational period determinations were done in MPO Canopus software version 10.8.6.3 (Warner 2021).Period determination was performed using the implemented Fourier Analysis of Light Curves (FALC) algorithm.
The obtained light curve of HL3 is shown in Figure 5.As the object was relatively faint (visual magnitude of 18.4-18.5),the data is noisy, and a unique period solution cannot be reliably derived.The nominal period solution is 5.3 h, but the alternative solutions of 2.6 and 1 h can not be ruled out.Therefore, additional observations of HL3 are needed to establish its rotation period firmly.
The amplitude of our nominal light curve solution is ∆m = 0.1 mag.Given that ∆m depends on the phase angle (Zappala et al. 1990;Gutiérrez et al. 2006), a somewhat larger amplitude is possible as well, but likely not much larger as we observed at phase angle α of ∼32 degrees.To account for the object's non-sphericity in our model, we adopted an amplitude of ∆m=0.15.Following the methodology described in Section 3.5, for the asteroid HL3, we found an axes ratio of a/c = 1.15.These yield a non-sphericity correction factor of ξ = 0.91.

Thermal inertia estimation of asteroid (152671) 1998 HL3
Though the obtained period solution is still unreliable and needs further verification, we obtained a preliminary estimate of the thermal inertia of HL3.We also take the opportunity to highlight once again the importance of an accurate period solution for reliable thermal inertia estimations with the ASTERIA model.Additionally, we demonstrated how the results may change when the variable thermal inertia model is used.
In the case of the HL3, rejecting one of the peaks is not trivial, and none of them could be rejected based on the currently available data.However, given that our aim here is primarily to demonstrate different aspects of the model and the relevance of the input parameters, we adopted the higher peak in our analysis.Adopting the lower peak does not change our overall conclusions.
In the upper panel of Fig. 6, we show the full thermal inertia distribution generated by our model using the parameters listed in Table 7.This distribution corresponds to thermal inertia of 506 +228 −121 J m −2 K −1 s −1/2 .If, however, we also consider the variable thermal inertia in our model, the result drops to 419 +177 −98 J m −2 K −1 s −1/2 .This change is relatively small, though the effect is larger than in the case of Bennu.It implies that the variable thermal inertia plays a limited role for asteroids with orbital eccentricities e ≤ 0.35, with HL3 being a marginal case. 12he lower panel of Fig. 6 demonstrates how our estimates of the asteroid HL3's thermal inertia vary with the rotation period.The results could change up to a factor of a few with this parameter.Keeping all the other parameters fixed, the TI increases with P, as expected.Therefore, the reliability of thermal inertia results obtained with ASTERIA also depends on the accuracy of the rotation period.This fact was also noted by Fenucci et al. (2021Fenucci et al. ( , 2023)), who exploited the fast rotation of two small NEAs 2011 PT and 2016 GE1 to constrain their thermal inertia to low values.

List of priority asteroids for rotation period determination
Having demonstrated how the (un)known rotation period affects the thermal inertia estimation by the ASTERIA model in Table 9, we provide a list of priority targets suitable for observations in the coming years.These objects have all the critical parameters used by the ASTERIA model determined, but the light curve.Therefore, estimating the amplitude of the light curves and determining the rotation periods would allow evaluation of the thermal inertia of those asteroids.As the final step in our analysis, we have selected asteroids from the JPL service for which all the relevant input parameters for the ASTERIA model, including the size and albedo, are available.We found 38 such objects, with all of them being near-Earth asteroids, and estimated their thermal inertia by means of our model.Among those 38 objects, 29 are characterized as potentially hazardous asteroids (PHA), making our results relevant also from the planetary defense point of view.
The complete list of objects, relevant input parameters, and thermal inertia estimations are given in Tables 10 and 11.In Fig. 7, we have shown thermal inertia vs. diameter values for our results and the other available literature estimates.It highlights the ASTERIA model's contribution to thermal inertia estimates of sub-km asteroids.Comments related to some individual cases are given below.
(441987) 2010 NY65 is a very interesting case.Along with the A 2 acceleration component, JPL also found an out-of-orbital plane, the A 3 acceleration component.Seligman et al. (2023) hypothesized that the A 3 component in small asteroids could be due to non-detected cometary-like activity.Therefore, the same reasoning places the object in the category of "dark comets".Its albedo is 0.071±0.014(Mainzer et al. 2011), in agreement with the albedo of cometary nuclei.The very low TI that we found may further support the link to comets, as values of TI below 100 J m −2 K −1 s −1/2 are common for comets (Groussin et al. 2019).However, based on the collected spectrophotometry data, the object was found to be S (or S V ) spectral type (Ieva et al. 2018;Perna et al. 2018), making its real nature quite puzzling.(3200) Phaeton is an active asteroid (Jewitt 2012), and the measured drift in the semimajor axis might not be entirely given by the Yarkovsky effect, because the activity could also contribute to the drift.Therefore, our results are based on the assumption that the possible contribution of the Phaeton's activity to its A 2 acceleration is negligible.This should be a reasonable assumption, given that the activity is expected to cause radial (A 1 ) or out-of-plane (A 3 ) acceleration components, which are yet to be detected according to JPL data.This object also has available thermal inertia estimates in the literature.The thermal inertia of Phaethon has been estimated to be Γ = 600 ± 200 by Hanuš et al. (2016) and Γ = 880 +580 −330 J m −2 K −1 s −1/2 by (Masiero et al. 2019).Both results are consistent with our result within 1 − σ uncertainties, suggesting that our model works well even in the case of active asteroids. 13nother quite intriguing case is (1566) Icarus.With an orbital eccentricity of e = 0.827, it is the second most eccentric orbit in our sample.Our estimation suggests a very low TI of asteroid Icarus.Indeed, an alternative solution from the model suggests the TI > 3000, which, however, seems unlikely.It is classified as a Q-type object (Binzel et al. 2019), suggesting the presence of a fresh surface material.The low thermal inertia, therefore, could be due to loosely packed, highly porous, material at the surface.The unweathered material could be exposed at the surface by different processes, such as planetary tidal forces, micrometeorite impacts, or thermal cracking.Though, in principle, any of these mechanisms could be involved, we believe the high temperature during the perihelion passage plays a key role, as perihelion (1566) Icarus is only 0.19 au from the Sun.Primitive C-type near-Earth objects are typically destroyed close to the Sun (Granvik et al. 2016), while rocky asteroids, as Icarus is supposed to be, are more likely to survive.Still, thermal cracking of the surface material could be a frequent process.
(163899) 2003 SD220 is a complex case.We found a very high TI, but several factors may affect this estimate.The asteroid has a very long rotation period of about 285 hours, and therefore, a high thermal inertia may be required for the Yarkovsky effect to work.However, the slow rotation may also suggest that it could be in a tumbling rotation state (Warner 2016).The maximum amplitude of its light curve is also very high, about 2.2 mag, indicating a very elongated shape.In such cases, our non-sphericity scaling approach described in Sec.3.5 may not be appropriate.However, we have also obtained the thermal inertia neglecting non-sphericity correction and obtained 974 +713 −424 J m −2 K −1 s −1/2 .Despite being somewhat lower, this supports the high thermal inertia solution for asteroid (163899) 2003 SD220.Contrary to the previous example, asteroid 2010 VK139 is a rapid rotator, with a period of only about 108 seconds.Our estimate of its TI of 213 +225 −107 J m −2 K −1 s −1/2 , potentially makes it another member of a group of small super fast rotators characterized by the relatively low thermal inertia, below ∼ 250 J m −2 K −1 s −1/2 .Still, the TI of VK139 is significantly higher than in the cases of 2011 PT, and especially 2016 GE1 (Fenucci et al. 2021(Fenucci et al. , 2023)).
In a more general context, an interesting aspect of thermal inertia investigation is whether typical thermal inertia values could be linked to specific asteroid spectral types (composition).Detailed analysis along these lines is beyond the purpose of this work.We only briefly discuss the case of V-type objects.
Four objects in our sample are classified as V-type asteroids.These are (3908) Nyx, (5604) 1992 FE, (297418) 2000 SP43, and (363599) 2004 FG11 (Thomas et al. 2014;Binzel et al. 2019).Jiang et al. (2020) recently studied the thermal inertia of 10 multi-km Vesta family members, presumed to be V-type asteroids, and found the average thermal inertia of 42 J m −2 K −1 s −1/2 , with all analyzed 10 objects having Γ 100 J m −2 K −1 s −1/2 .Among four V-types studied in our work, three of them have nominal terminal inertia Γ 200 J m −2 K −1 s −1/2 which are consistent with low values found by Jiang et al., taking into account that our estimates refer to near-Earth objects orbiting at smaller heliocentric distance.The asteroid (5604) 1992 FE has also been classified as V-type, but its thermal inertia has been found by MacLennan & Emery (2021) to be quite high of 1000 +2200 −600 , in agreement with our nominal result of 922 +769 −332 J m −2 K −1 s −1/2 .That suggests a large range of possible thermal inertia values also associated with the V-type asteroids, as found in the case for C-and S-type classified objects (Harris & Drube 2020).Still, we note that, in our case, an alternative solution for the (5604) 1992 FE thermal inertia of 153 +85 −71 J m −2 K −1 s −1/2 cannot be ruled out.If the latter estimate is correct, that would indicate that low thermal inertia values preferentially characterize V-type objects.Hung et al. 2022, and references therein).Note that numbers from the literature include multiple instances of the same asteroid.Additionally, our two previous estimates for asteroids 2011 PT (Fenucci et al. 2021) and 2016 GE1 (Fenucci et al. 2023) are shown as red squares and indicated in the plot.

SUMMARY AND CONCLUSIONS
Here, we presented the new Asteroid Thermal Inertia Analyzer (ASTERIA) model and publicly available corresponding software for determining the surface thermal inertia of asteroids.The model allows thermal inertia estimation based mostly on population-modeled input parameters in its basic variant.However, as in general cases, the model may not A general advantage of the ASTERIA model is that it may be applied to smaller asteroids than the thermophysical modeling, because the Yarkovsky effect is more substantial in smaller objects and, therefore, easier to detect.Additionally, thermophysical modeling requires thermal infrared observations and good shape models, which are currently challenging for asteroids below some 100 meters in size.In this respect, we note that among all the asteroids' thermal inertia estimates available in the literature, the results for the three smallest ones are obtained with our model (see Fig. 7).
Based on the astrometric measurements and the detection of Yarkovsky-induced acceleration of orbital motion, it is primarily independent of the most widely used approach for the asteroid thermal inertia estimations based on thermophysical modeling.As such, the ASTERIA may also serve as a benchmark test to independently verify the results derived from the thermophysical modeling, one of the longstanding challenges of TPM models (Hung et al. 2022).
As we found that knowledge of a rotation period is highly relevant for the ASTERIA model to work correctly, we provide a list of high-priority targets (Table 9) and encourage the observers to attempt to obtain light curves of those objects.Determining the rotation periods of those asteroids would allow the thermal inertia estimations by the ASTERIA model.
Finally, we have identified a set of 38 NEAs for which all the input parameters critical for the ASTERIA model to work reliably are available and presented the new thermal inertia for those objects.Among these 38 NEAs, 29 are classified as potentially hazardous asteroids (PHA, Table 10).It makes our results highly relevant from the planetary defense point of view.Our sample of new thermal inertia estimates also includes 31 sub-km-sized asteroids, while there are only 17 other literature values in this size range, highlighting the importance of the ASTERIA model for determining the surface thermal inertia of small asteroids.
Among the limitations, we recall that if the input parameters are not well-defined and have considerable uncertainty, this propagates into significant uncertainties of the thermal inertia estimates by our model.An example of this is the asteroid (1865) Cerberus.Though we obtained results formally in agreement with the reference literature estimate, the situation highlights the model's sensitivity to the Yarkovsky drift determination.Also, for likely elongated objects with large light curve amplitudes, such as Cerberus, our non-sphericity scaling could be rough, contributing to overall uncertainty.Nevertheless, we showed that the obtained results are even in these cases compatible with other estimates.Therefore, results for such objects could be helpful but should be interpreted cautiously.Future better characterization of the input parameters will further improve the precision of our results.Also, our model could be somewhat affected by spatial thermal inertia variation across the surface of an asteroid.In case of significant spatial variations, an additional error of up to about 5% could be introduced.For instance, in the case of asteroid Bennu, we found up to a 2% shift in the result due to spatial variation in thermal inertia from the equator to the polar regions.
Applicability of our model is still limited to the near-Earth asteroids, as the Yarkovsky detections are yet unavailable for main-belt asteroids (Tanga et al. 2023).

Figure 1 .
Figure 1.Coefficients k1, k2, and k3 as a function of the scaled radius R ′ .

Figure 2 .
Figure 2. PDF of the obliquity distribution of the NEAs population.

Figure 3 .
Figure 3. Distribution of asteroid Bennu's thermal inertia Γ generated with the ASTERIA model.The figure shows the results the nominal settings after considering the spatial variations in Bennu's surface thermal inertia.See text for additional details.The median values of Γ and their corresponding lower and upper 1 − σ uncertainties are shown for each of the peaks.

Figure 4 .
Figure 4. Thermal inertia values derived from thermophysical modeling listed in the literature compared to values derived from the ASTERIA model.The blue points are the results for ten NEAs used for model validation purposes, while the black triangle shows the result for Bennu.The red line is the line of equality where the same results should appear.

Figure 5 .
Figure 5. Data for two nights photometric observations of asteroid (152671) 1998 HL3 and the corresponding light curve for the nominal period solution of P = 5.3 h, plotted with MPO Canopus.

Figure 6 .
Figure 6.The thermal inertia estimates for asteroid (152671) 1998 HL3 generated with the ASTERIA model.The upper panel shows the full distribution obtained for our nominal rotation period solution.The lower panel shows how the results change with the rotation period, highlighting the importance of reliable determination of this parameter.The median values of Γ and their corresponding lower and upper 1 − σ uncertainties are shown.The vertical grey area in the lower panel indicates the result for the nominal period solution.See text for additional details.

Figure 7 .
Figure 7. Thermal inertia vs. diameter values for objects with available thermal inertia estimates.Our sample of thermal inertia estimates for 38 asteroids (black triangles) is shown along with the thermophysically derived literature values (orange circles) (see MacLennan & Emery 2021;Hung et al. 2022, and references therein).Note that numbers from the literature include multiple instances of the same asteroid.Additionally, our two previous estimates for asteroids 2011 PT(Fenucci et al. 2021) and 2016 GE1(Fenucci et al. 2023) are shown as red squares and indicated in the plot.

Table 9 .
Priority targets for light curve determination that would allow thermal inertia estimation by the ASTERIA.Objects suitable for observation (visual magnitude < 22, solar elongation > 90 • ) until 2027 are listed.For each object, the month of its brightest visual magnitude has been given.

Table 10 .
The orbital properties of 38 NEAs used in the ASTERIA model for thermal inertia estimations.Source: JPL-Small-Body Database Lookup.