Thermal properties of large main-belt asteroids observed by Herschel PACS

Non-resolved thermal infrared observations enable studies of thermal and physical properties of asteroid surfaces provided the shape and rotational properties of the target are well determined via thermo-physical models. We used calibration-programme Herschel PACS data (70, 100, 160 $\mu$m) and state-of-the-art shape models derived from adaptive-optics observations and/or optical light curves to constrain for the first time the thermal inertia of twelve large main-belt asteroids. We also modelled previously well-characterised targets such as (1) Ceres or (4) Vesta as they constitute important benchmarks. Using the scale as a free parameter, most targets required a re-scaling $\sim$5\% consistent with what would be expected given the absolute calibration error bars. This constitutes a good cross-validation of the scaled shape models, although some targets required larger re-scaling to reproduce the IR data. We obtained low thermal inertias typical of large main belt asteroids studied before, which continues to give support to the notion that these surfaces are covered by fine-grained insulating regolith. Although the wavelengths at which PACS observed are longwards of the emission peak for main-belt asteroids, they proved to be extremely valuable to constrain size and thermal inertia and not too sensitive to surface roughness. Finally, we also propose a graphical approach to help examine how different values of the exponent used for scaling the thermal inertia as a function of heliocentric distance (i.e. temperature) affect our interpretation of the results.


Introduction
Non-resolved observations of asteroid thermal infrared (IR) emission provide information about the asteroid sizes and the thermal properties of their surfaces. By taking into account asteroid shape, rotation, and the geometry of the observations, thermo-physical models (TPMs) can be used to compute surface temperatures and fit thermal properties, such as thermal inertia, to the IR data. The number of asteroids with a thermo-physical characterisation has increased greatly over the last two decades thanks to the ever growing number of available shape and rotational models -which are necessary input for the TPM -and great observational efforts, both in the visible and the thermal IR (for recent reviews see e.g. Delbo et al. 2015;Ďurech et al. 2015;Mainzer et al. 2015, and references therein).
The number of available shape models, by which we refer to both shape and rotational properties, is dominated by the several hundred models derived from the inversion of non-resolved optical light curves (e.g. Ďurech et al. 2010;Ďurech et al. 2018;Hanuš et al. 2011Hanuš et al. , 2013 following the method by  (see also Kaasalainen et al. , 2002. Increasingly more sophisticated algorithms have begun to combine various data types in the inversion (Carry et al. 2010;Ďurech et al. 2011;Viikinkoski et al. 2015; Bartczak & Dudziński 2018), even thermal IR data Müller et al. 2017), and there are now dozens of models based also on stellar occultations, radar, and adaptive optics observations (e.g. Ďurech et al. 2015;Benner et al. 2015, and references therein).
As the availability of asteroid shape models continues to increase, we can readily make use of a large collection of thermal IR data from several catalogues like those provided by the Infrared Astronomical Satellite (IRAS), AKARI, or the Wide-field Infrared Survey Explorer (WISE). Up to 2018, targeted observations from ground-and space-based facilities (e.g. Müller 2002;Müller & Blommaert 2004;Emery et al. 2006;Müller et al. 2017;Landsman et al. 2018) and/or from all-sky thermal IR surveys (e.g. Delbo' & Tanga 2009;Alí-Lagoa et al. 2014;Rozitis et al. 2014;Hanuš et al. 2015Hanuš et al. , 2016Bach et al. 2017;Marciniak et al. 2018) had been used to model ∼100 asteroids with TPMs. With the recent addition of another 100 (Hanuš et al. 2018b), WISE data are now the single largest source of asteroid thermal inertias.
Our aim in this work is to provide a thermo-physical characterisation of main-belt asteroids (MBAs) with Herschel Photodetector Array Camera and Spectrometer (PACS) data taken during the Asteroid Preparatory Programme for Herschel, ASTRO-F and the Atacama Large Millimiter/submillimiter ar-1 Large MBAs were too bright for WISE and their partially saturated fluxes might not be optimal for thermo-physical modelling. On the other hand, partial saturation can be corrected for and reasonably suited for thermal modelling purposes (e.g. Masiero et al. 2011;Mainzer et al. 2015). 2 http://archives.esac.esa.int/hsa/legacy/UPDP/SBNAF_MBA/SBNAF_MBA_Release_Note.pdf 3 For (21) Lutetia (see Sect. A.4), we also included most of the data analysed in O' Rourke et al. (2012). provide all relevant details about the production of the IRDB 4 , including the colour correction approach, sources and references of each catalogue, and useful auxiliary quantities such as observation geometry, and light-travel time. For completeness, Table B.1 provides all previously unpublished PACS observations used in this work.

Thermo-physical model implementation
Our TPM implementation is the one used by  based on that of Delbo' et al. (2007) and Delbo' & Tanga (2009). Said version was upgraded to account for the effects of shadowing, but global self-heating was neglected for this work (see e.g. the discussion about average view factors by Rozitis & Green 2013). Below we provide only a brief summary of the technique that we use and the approximations that we make. Our methodology was described in more length in Marciniak et al. (2018) and Marciniak et al. (2019). There, details about the thermo-physical modelling of each target were provided in separate sections. Here, we only present some relevant plots in the main text and include all TPM-related plots and additional comments in Appendix A.
We take a shape model as input for our TPM (see Table 1) with the main goal of modelling the surface temperature distribution at epochs at which we have thermal IR observations and constrain the target's diameter, thermal inertia and, whenever possible, surface roughness. From the known geometry of observation we identify which surface elements (usually triangular facets) were visible to the observer at each epoch and compute the model fluxes from the surface temperature distribution.
To account for heat conduction towards the subsurface, we solve the 1D heat diffusion equation for each facet and we use the Lagerros approximation when modelling the surface roughness via hemispherical craters of different depth coveraging 0.6 of the area of each facet (Lagerros 1996, Lagerros 1998, Müller & Lagerros 1998, Müller 2002. We also consider the spectral emissivity to be 0.9 regardless of the wavelength (see e.g. Delbo et al. 2015).
For each target, we estimated the Bond albedo (necessary input) as the average value obtained from the different radiometric diameters available from AKARI and/or WISE (Usui et al. 2011;Alí-Lagoa et al. 2018;Mainzer et al. 2016), and all available H-G, H-G 12 , and H-G 1 -G 2 values from the Minor Planet Center, Oszkiewicz et al. (2011), or Vereš et al. (2015. This approach leaves us with two free parameters, the scale of the shape (interchangeably called the diameter, D) and the thermal inertia (Γ). The diameters and other relevant information related to the TPM analyses of our targets are provided in Table  2. Whenever the data are too few to provide realistic error bar estimates, we report the best fitting diameter so that the models can be scaled and compared to the scaling given by the occultations. On the other hand, if we have multiple good-quality thermal data (with absolute calibration errors below 10%) then this typically translates to a size accuracy of around 5% as long as the shape is not too extreme and the spin vector is reasonably well established. This rule of thumb certainly works for large MBAs like the Gaia mass targets. We do not consider the errors introduced by the pole orientation uncertainties or the shapes (see Hanuš et al. 2015 andDudziński 2019), so our TPM error bars estimates are lower values.  Table 2 summarises the thermophysical properties obtained for each target in our sample. For comparison, we also include the results for the spheres with the same rotational properties as the shape models. All plots produced during the modelling and further details about some selected targets are provided in Appendix A. In this section we provide a summary of the results and focus on several aspects relevant to the sample and/or a larger catalogue of thermo-physical properties retrieved from the literature. For example, the thermal inertia versus size plot ( Fig. 1) shows that our new thermal inertias fall within the same range of those of other large MBAs found in previous works.

Results
Most ADAM shape models only required a small re-scaling of the size to fit the data with lowχ 2 min , which we take as an additional confirmation of their high quality. However, in the cases of (65) Cybele, (18) Melpomene, and (54) Alexandra, we required re-scalings of the order of 10% to fit the data, which are larger than expected for the available high-quality thermal data and shape models (e.g. Delbo et al. 2015). Thus, regardless of their quality, it is still worthwhile keeping the scale as a free parameter when performing thermo-physical modelling using already-scaled shape models in order to check potential issues. For example, the ADAM shape of Pallas with a fixed scale could not reproduce well PACS data taken closer to pole-on, as illustrated in Figs. 2 and 3. The first one shows the observationto-model ratios (OMR) versus aspect angle for our best fit, that is Γ = 30 J m −2 s −1/2 K −1 and a re-scaling of +3%. The second one shows the best fit obtained with a fixed scale, which led to Γ = 8 J m −2 s −1/2 K −1 . Although this fit is also formally acceptable, the OMRs are only close to one at equatorial views (aspect angles around 90 degrees). No other combination of thermal inertia and roughness could bring the other ratios closer to unity, which suggests that the shape model could be less accurate at pole-on views. In this particular case, this could also explain the lower thermal inertia needed to fit the data when we kept the scale fixed, because the data taken at higher phase angles were coincidentally the ones taken at non-equatorial sub-observer latitudes, which makes it more difficult for the shape model to reproduce the thermal phase curve (see Fig. 4).
We found formally good fits and thermal properties within the expected range for the two targets with SAGE models available for our study as well, albeit with slightly higher error bars. Additional SAGE models scaled using stellar occultation chords and further discussion are provided by Podlewska-Gaca et al. (2020). It is worth mentioning that the TPM did not help us to favour any of the two mirror shape solutions of (20) Massalia. For each shape model, we also fitted the data using spheres with the corresponding rotational parameters with the aim of comparing the resulting thermo-physical parameters (see Table  2). On the one hand, this approximation seems to provide reasonable estimates for the diameters that would be expected for large objects with relatively low-amplitude light curves. On the other hand, spheres lead to scales ∼5% larger on average than the ADAM shapes, and several of the thermal inertia values were up to an order of magnitude higher than those of the corresponding shape models. With the current sample we could not identify any single cause that should lead to such systematic effects, but perhaps a future larger sample could help explore this further. Also, this could be used to create a reliable benchmark to estimate thermal inertias of objects with more limited shape models.

Thermal inertia, pole-on views, and rotational variability
Only those asteroids with highly oblique rotational axes (or alternatively, unusually high orbital inclinations) can be seen with aspect angles close to pole-on, that is, 0 or 180 degrees for the north and the south pole. In these cases, seasonal effects on the surface temperature distribution are expected, since the depth at which the heat wave penetrates is higher at close to pole-on illuminations. In principle, data taken in such circumstances could require a model accounting for thermal inertia variation with depth (such as those used for the Moon; Hayne et al. 2017), but our OMR versus aspect angle plots do not suggest problems in fitting data taken at close to pole-on views (e.g. see Fig. 2). We did not find any correlation between Γ and the pole ecliptic latitude either. The vast majority of OMR versus rotational phase plots suggest that the shape models and the assumption of homogeneous thermal properties throughout the surface can still reproduce the IR data well (see , third panel from the top).

Constant thermal inertia with temperature
Given the dependence of conductivity on temperature, thermal inertia obtained from IR data taken at an average heliocentric distance r needs to be normalised to some reference heliocentric distance (i.e. temperature), usually 1 au, following where α = −0.75 if we consider a radiative conduction term in the thermal conductivity κ. Marsset et al. (2017) analysed thermal IR data of asteroid (6) Hebe taken at a wide range of heliocentric distances and found indications that higher thermal inertias fitted low-r data better. However, even though our sample includes some highly eccentric asteroids, our OMR versus heliocentric distance plots show no systematics or biases due to the Article number, page 3 of 28 A&A proofs: manuscript no. pacs_accepted_arxiv Our new values (1σ error bars). Compilation  Fig. 1. Thermal inertia normalised at 1 au versus size. The PACS targets (black symbols) follow the trend in the compilation by Delbo' et al. and Hanuš et al. (2018b). Observation to model ratios vs. aspect angle for the best fitting TPM model to the PACS data of (2) Pallas. The best fit was obtained by Γ = 30 J m −2 s −1/2 K −1 and a re-scaling of 1.03 (see Table 2). assumption that Γ = constant (see e.g. Fig. 5), so it seems that the temperature range covered by most MBAs might be such that this effect is not critical. Rozitis et al. (2018) studied three highly eccentric near-Earth objects observed by WISE at widely different r and fitted the exponent α separately for each case. This work clearly showed that the physical interpretation of thermal inertia and the comparison between different objects and populations without good constraints for α is quite complicated. Here we propose a graphic approach to analyse a large catalogue of thermal inertias with the necessary assumption that all asteroids can be modelled using a single value of α, which is arguable. The aim, nonetheless is to further explore the impact of the value of this parameter on our interpretation of the thermal inertias with a larger catalogue of values (see also Szakáts et al. 2020).
Instead of normalising the Γ-values to 1 au, we can plot the thermal inertias compiled in Delbo et al. (2015) plus those in Hanuš et al. (2015), Hanuš et al. (2018a), Hanuš et al. (2018b), Marciniak et al. (2018), and Marciniak et al. (2019) as a function of heliocentric distance and compare them with different curves with Γ 1au = 10, 50, . . . , 2000 and a given value of α using Eq. 1. With such a plot, we can visually estimate Γ 1au for each asteroid as the curve going through its corresponding point. On the one hand, α = −0.75 ( Fig. 6) leads to the conclusion that most D > 10-km asteroids (i.e. the green, brown and yellow symbols) have thermal inertias at 1 au between 10 and 250 J m −2 s −1/2 K −1 (dotted line), whereas those of most sub-km objects (pink symbols) lie between ∼200 and 1000 J m −2 s −1/2 K −1 . On the other hand, with α = −2.2 (Fig. 7), we would conclude that there is a much higher degree of thermal inertia variability amongst D > 10 km asteroids, since Γ 1au spans the range < 50-2000 J m −2 s −1/2 K −1 .
Finally, we also examined plots like Figs. 6 and 7 but using the rotational period instead of size in the colour axis: they do not show any appreciable excess of slow rotators in the regions of higher Γ 1au or, conversely, fast rotators in the regions with lower Γ 1au , regardless of the value of α. This is consistent with the conclusions of Marciniak et al. (2019), whose sample shows that slow rotators do not always present higher thermal inertias (cf. Marciniak et al. 2018, who had a smaller sample, and Harris & Drube 2016 based on an empirical relation between thermal inertia and the infrared beaming parameter η used in the near-Earth asteroid thermal model of Harris 1998

Discussion
One of the original motivations for this work within the framework of the SBNAF project was to test how the TPM could be used to evaluate the quality of shape models. While the comparison of the results with those obtained for a sphere is always a useful benchmark for discussion of bad or borderline acceptable fits, we found a relatively small systematic offset of the equivalent diameters obtained for the spherical models. However, after a thorough examination of each individual case, we failed to identify any cause that could explain why the spheres would occasionally provide good approximations of the thermal inertia, while being up to an order of magnitude higher in other cases. Neither the irregularity of the shape nor the particular orientation of the spin axis seem to produce worse results, at least within our sample.
We fully neglect the uncertainties in the shape models even though they would certainly contribute to the final error bars of our TPM parameters (see Hanuš et al. 2015). In this sense, our error bars are lower limits. One of the reasons for this limitation is that the methods to estimate errors in the shape are either labour intensive, computationally costly, or both, but there are already some works that explore this direction. Hanus et al. (2015) proposed a way to assess the effects of shape uncertainty in TPM modelling by bootstrapping the visible data used to obtain a shape model from light-curve inversion and producing a set of shape models, each one of which is in turn modelled with a TPM to provide statistics for the inferred TPM parameters. More recently Bartczak & Dudziński (2019) used an approach featuring millions of slightly perturbed shape clones to examine how the uncertainties and inaccuracies of the inversion models map over the whole surface. As the number of shape models determined from adaptive optics observations and stellar occultations continues to rise, we could reach a point where we have sufficient ground-truth information about the shapes to have a more empirical estimate of the errors of thermal inertia and diameter. Podlewska-Gaca et al. (2020) provide a larger sample of SAGE models scaled with stellar occultation chords and, whenever possible, thermo-physical modelling.
Finally, we note that our thermal inertias for Ceres and Pallas are higher than previous values (Müller & Lagerros 1998; see also Secs. A.1 and A.2). Thermo-physical modelling of Uranian satellites by Detre et al. (submitted to Astronomy & Astrophysics) suggests this could be a small systematic trend: the radiometric diameters are 3-5% larger than the ones derived from direct methods. We have not identified a clear unique cause for such an offset, but it could be related to the absolute flux calibration of the PACS data. On the other hand, we obtained only slightly higher values than O' Rourke et al. (2012) for Lutetia and Capria et al. (2014) for Vesta. Another source of error might be the assumption of a constant emissivity of 0.9, but the need for much lower thermal emissivities has so far been pointed out at significantly longer wavelengths (Müller & Lagerros 1998, Müller & Barnes 2007).

Conclusions
In this work we derive thermo-physical properties of 18 large main-belt asteroids, 12 of which did not have any previous constraints on thermal inertia. This was enabled by previously unpublished Herschel PACS data (Table B) and recently available state-of-the art non-convex shape models (the term also includes the rotational properties). Most of our targets' shape models (see Table 1) were derived from a combination of adaptive optics, occultation and optical light-curve data using the ADAM algorithm ; three of them, which we used as benchmarks, from direct images from spacecraft, and two of them from only optical light curves via SAGE (Bartczak & Dudziński 2018).
We find that the ADAM shapes can reproduce the thermal IR data in spite of the usual simplifying assumptions of the thermo-physical modelling -constant properties over the surface, and thermal inertia also independent of depth and temperature 5 -even in the case of (10) Hygiea, which has been recently shown to have a variable albedo over the surface (Vernazza et al. 2019). We optimised the scale of the shape models (i.e. the scale was a free parameter in our thermo-physical model, as usual) and found that most cases required an average re-scaling of 5%, which serves as a cross-validation. Nonetheless, there were notable exceptions, especially (65) Cybele (see Table 2), that require ∼10% scaling, which means these targets are interesting for follow-up observations and modelling.
From the example of (2) Pallas, a target extensively observed at a wide range of aspect angles due to the high obliquity of its rotation axis, we also examined how potential inaccuracies in the shape can bias the results of thermo-physical modelling when the scale of the shape model is not allowed to vary, so we suggest that results from both approaches (fixed and fitted scale) should always be compared and a careful examination of the modelling results should be performed on a case-by-case basis. Also, more work is needed to quantify and parameterise shape model errors (Hanuš et al. 2015;Bartczak & Dudziński 2019) so that they can be propagated into the thermo-physical properties derived from them in a practical way.
In addition, we found relatively low thermal inertias in the same range as previous results , and references therein) for similarly-sized asteroids, which continues to support the notion that these bodies are covered by fine regolith. This, rather than composition, is the dominant effect governing the thermal emission of the large asteroids. Although the peak of the emission of main-belt asteroids is closer to the 10-micron region of the spectrum, high-quality data at longer wavelengths like PACS data are also extremely useful to determine goodquality sizes and thermal inertias; objects with large data sets especially benefit from the fact that surface roughness is not as dominant as thermal inertia at this wavelength range.
Acknowledgements. The research leading to these results has received funding from the European Union's Horizon 2020 Research and Innovation Programme, under Grant Agreement number 687378 (SBNAF). C.K., R.S., A.F-T., and G.M. Table 2. Summary of TPM results. D 0 is the sphere-equivalent diameter of the ADAM or in-situ shape models. Spherical model refers to a sphere (∼3000 facets) with the same spin axis. The symbols D, Γ (SI units = J m −2 s −1/2 K −1 ) and ρ denote the best-fitting diameter, thermal inertia and surface roughness (rms) of the corresponding model. To normalise Γ at 1 au (Γ 1au ), we took the mid-point (r) between the shortest and longest heliocentric distance at which the data were taken. The surface roughness (rms) were not constrained at the 1σ level unless otherwise stated. For more information we refer to Appendix A and 0.6 Low χ 2 m but unrealistically high Γ, perhaps because shape is elongated Article number, page 7 of 28 Appendix A: Thermo-physical model analysis In this section we provide all observation-to-model ratio (OMR) plots to help further examine the fits. Table A.1 links each target to its corresponding plots. Some targets warrant a short subsection with relevant comments in addition to those of Table 2.  The 180 PACS data points are fitted with a very lowχ 2 min of 0.2 with Γ = 25 J m −2 s −1/2 K −1 , an extremely high roughness of rms∼1, and a rescaling of about +2% of the original ADAM shape. This thermal inertia is higher than but compatible within the error bars with classical TPM results (10 ± 10 J m −2 s −1/2 K −1 ; Müller & Lagerros 1998) and Dawn-based analyses (< 20 J m −2 s −1/2 K −1 ; Rognini 2018). If we fit the data while keeping the Dawn stereo-photoclinometry (SPC) shape fixed, we obtain a lower thermal inertia of 6 +9 −6 J m −2 s −1/2 K −1 and roughness values lower than rms∼0.35 are rejected at the 3σ level.
We find systematically higher thermal inertias when we model only PACS data and optimise the diameter, although the results are not statistically significantly different (1-σ regions in the χ 2 vs. Γ plots overlap). Actually, for targets with such a rich PACS data set as Ceres, we can fit data from the three PACS filters separately. We found Γ = 35, 40, and 50 J m −2 s −1/2 K −1 with the 70, 100, and 160 µm data. Although the trend is not statistically significant -our χ 2 minima still overlap within the ranges of formally acceptable fits -it could be caused by the fact that the longer wavelengths probe deeper and perhaps more compacted layers of the subsurface. However, it could also be an artefact related to our assumption of a constant emissivity of 0.9 if the (spectral) emissivity at the PACS wavelengths was lower, but previous work by Müller (2002) suggests that the effect should only be significant at wavelengths much longer than that of PACS. Indeed, Müller et al. (2014) were able to reproduce even Herschel Spectral and Photometric Imaging Receiver (SPIRE) data at 250, 350 and 500 µm with the same modelling approach we use. Nonetheless, we leave it to future work to revisit the data with additional observations and better constraints on spectral emissivity.

Appendix A.2: (2) Pallas
The ADAM shape fits the 80 PACS data with a very lowχ 2 min (lower than 0.1) with extremely high roughness and a thermal inertia of 30 J m −2 s −1/2 K −1 . Nevertheless, we find indications of small shape model inaccuracies given it cannot reproduce all data equally well when we fix the scale (see the main text).

Appendix A.3: (3) Juno
The OMR plots present systematics within the 10% error margins. From the aspect angle plot, we infer that the northern hemisphere part of the shape model could have an artefact, and the wavelength plot shows a slight "convex up" curvature. The SAGE model provides a borderline formally acceptable fit and shows more scatter in the OMR plots and larger parameter error bars, which means the shape is not optimal (the rotational parameters are virtually the same). The trends in the OMR plots are slightly blurred to the eye due to the higher scatter. Nonetheless, the best-fitting values of size and gamma are fully compatible within the error bars, but the ADAM ones are still better constrained. The same can be said about the scatter in the plots corresponding to the sphere (which does slightly better than SAGE in terms of χ 2 , but not statistically significantly).

Appendix A.4: (21) Lutetia
For this target we incorporated most of the data featured in O' Rourke et al. (2012) in our analysis, for example Spitzer Infrared Array Camera (IRAC) fluxes, with the notable exception of Herschel SPIRE fluxes, which we did not model in this work. As expected, our results for Lutetia are fully consistent with the previous work. O'Rourke et al. suggest that a localised inaccuracy in the shape model, which is after all based on data that did not cover 100% of Lutetia's surface, is responsible for the inability of the TPM model to fully reproduce the October 17 IRAC thermal light curve (see their Fig. 4). Here we also considered an offset in the rotational zero-phase as a possible explanation, so we repeated the analysis with a delay of 20 degrees. As Fig. A.11 shows, the maxima and minima of the model were better aligned with both IRAC light curves but the overall flux levels of the October 17 one were still not matched. Nevertheless, this mismatch is small, on average 5% (pink OMRs at aspect angle close to 40 degrees in the fourth panel of Fig. A.12), and a possible localised inaccuracy in the shape-model still cannot be ruled out.

Appendix B: Herschel PACS observations
We refer the reader to the caption of    Müller et al. (2005b), Nielbock et al. (2013), Balog et al. (2014), andMüller et al. (2014) for details on the data set and observing modes (Obs. mode), and Kiss et al. (2014) on the reduction approach, which required non-standard strategies for a subset of the observations. These fluxes will also be retrievable from the SBNAF database (Szakáts et al. 2020) along with the other data and all relevant references. The start and end times of the observations are denoted by t s and t e , λ is the reference wavelength of the band, f λ and σ f λ the calibrated in-band flux and its error, f λ the colour-corrected flux density, σ f λ the absolute calibration error, r the heliocentric distance, ∆ the distance to the observer, and α phase angle. The sign of the phase angle is given by that of the vector product between versors from the asteroid to the Sun and the observer.

Target
Obs. ID Obs. mode t s (JD) t e (JD) Band