© ESO, 2017

1. Introduction

The study of the visible photometric variation of solar system minor bodies enables us to determine optical light curves (flux or magnitude versus time), for which essential parameters are the peak-to-peak amplitude and the rotational period of the object. Short-term photometric variability of trans-Neptunian objects (TNOs) and Centaurs can be shape-driven (i.e. Jacobi-shaped rotating body, such as Varuna, Jewitt & Sheppard 2002) or be causes by albedo contrasts on the surface of a Maclaurin-shaped rotating spheroid (e.g. the Pluto case). Combinations of shape and albedo effects are also possible and very likely occur within TNOs (e.g. Makemake) and Centaurs. Contact-binaries can also produce short-term photometric variability within the TNO and Centaur populations. The largest amplitudes are usually associated with Jacobi shapes and the smallest ones with Maclaurin shapes with highly variegated surfaces. Large amplitudes can also be associated with contact-binaries and small amplitudes with objects with rotational axes close to pole-on. If we know the rotational properties (i.e. rotation period and amplitude) of a Jacobi shaped object, it is possible to derive the axes ratio of the ellipsoid (i.e. a shape model) and also a lower limit for the density (Jewitt & Sheppard 2002; Sheppard et al. 2008,and references therein), assuming the object is in hydrostatic equilibrium (Chandrasekhar 1987) with a certain aspect angle. On the other hand, if we suspect that the object has a Maclaurin shape we can derive a shape model from the rotational period, but it is needed to assume a realistic density in this case. The majority of the TNOs/Centaurs (~70%) present shallow light curves (amplitudes less than 0.15 mag), which indicates that most of them are Maclaurin-shaped bodies (Duffard et al. 2009; Thirouin et al. 2010). For a couple of special cases (only for Centaurs until now) where more information is available (i.e. long-term changes in light curve amplitudes), the position of the rotational axis can be derived or at least constrained (Tegler et al. 2005; Duffard et al. 2014a; Ortiz et al. 2015; Fernandez-Valenzuela et al. 2017).

Complementary to optical light curves, thermal light curves are a powerful tool to obtain additional information about physical and, in particular, thermal properties of these bodies. At first order, immediate comparison of the thermal and optical light curves enables us to differentiate between shape-driven light curves (the thermal light curve is then correlated with the optical one) and light curves that are the result of albedo markings (the two light curves are anti-correlated). Furthermore, quantitative modeling of the thermal light curve enables us to constrain the surface energetic and thermal properties, namely its bolometric albedo, thermal inertia, and surface roughness (e.g. Pluto, see Lellouch et al. 2016,and references therein). In a more intuitive way, in the case of positively correlated thermal and optical light curves, it is also possible to constrain the thermal inertia using the delay between the thermal and optical light curves and their relative amplitudes.

Until recently, only a few thermal light curves of outer Solar System minor bodies (or dwarf planets) have been obtained. Pluto’s thermal light curve was detected with ISO/PHOT, Spitzer/MIPS, and /IRC at a variety of wavelengths longwards of 20 μm (Lellouch et al. 2000, 2011). More recently, Pluto was also observed with the Herschel Space Observatory (Pilbratt et al. 2010) Photodetector Array Camera and Spectrometer (PACS: Poglitsch et al. 2010) and with the Spectral and Photometric Imaging Receiver (SPIRE: Griffin et al. 2010), which provides thermal light curves at six wavelengths: 70, 100, 160, 250, 350, and 500 μm (Lellouch et al. 2016). A marginal thermal light curve of dwarf planet Haumea was reported with Spitzer-MIPS (Lockwood & Brown 2009; Lockwood et al. 2014). The definite detection of Haumea’s thermal light curve was obtained with Herschel/PACS within the Science Demonstration Phase (SDP, Lellouch et al. 2010). Other tentative thermal light curves of TNOs/Centaurs that were observed with Herschel/PACS (i.e. Varuna, Quaoar, Chiron, Eris) were also presented for the first time (Santos-Sanz et al. 2014) and will be published separately (e.g. Kiss et al., in prep.).

Here, we present thermal time series photometry of the dwarf planet Haumea and the Plutinos 2003 VS₂ and 2003 AZ₈₄ that were taken with Herschel/PACS using its 3-filter bands, which are centred at 70, 100, and 160 μm (hereafter blue, green, and red bands, respectively). In the case of Haumea, we present additional and improved data and we merge them with the SDP observations that were originally presented in Lellouch et al. (2010). The thermal time series of 2003 VS₂ and 2003 AZ₈₄ are presented here for the first time. The Herschel/PACS observations are described in Sect. 2, the data reduction and photometry techniques applied are detailed in Sect. 3. Data are analyzed, modeled, and interpreted in Sect. 4. Finally, the major conclusions of this work are summarized and discussed in Sect. 5.

2. Observations

The observations presented here are part of the project “TNOs are Cool”: a survey of the trans-Neptunian region, a Herschel Space Observatory open time key programme (Müller et al. 2009). This programme used ~372 h of Herschel time (plus ~30 h within the SDP) to observe 130 TNOs/Centaurs, plus two giant planet satellites (Phoebe and Sycorax), with Herschel/PACS; 11 of these objects were also observed with Herschel/SPIRE at 250, 350, and 500 μm (see Fornasier et al. 2013), with the main goal of obtaining sizes, albedos, and thermophysical properties of a large set of objects that are representative of the different dynamical populations within the TNOs. For PACS measurements, we used a range of observation durations from about 40 to 230 min based on flux estimates for each object. Results of the PACS and SPIRE measurements to date have been published in Müller et al. (2010), Lellouch et al. (2010), Lim et al. (2010), Santos-Sanz et al. (2012), Mommert et al. (2012), Vilenius et al. (2012), Pál et al. (2012), Fornasier et al. (2013), Lellouch et al. (2013), Vilenius et al. (2014), Duffard et al. (2014b)¹. In addition, four bright objects (Haumea, 2003 VS₂, 2003 AZ₈₄ and Varuna) were re-observed long enough to search for thermal emission variability related to rotation (i.e. thermal light curve). This paper presents results for the first three. Other objects (Pluto, Eris, Quaoar, Chiron) were also observed outside of the “TNOs are Cool” programme to search for their thermal light curve.

The general strategy we used to detect a thermal light curve with Herschel/PACS was to perform a long observation covering most of the expected light curve duration, followed by a shorter follow-on observation, which enabled us to clean the images’ backgrounds, as explained in Sect. 3.1.

Dwarf planet (136108) Haumea was observed twice with Herschel/PACS in mini scan maps mode at 100 and 160 μm. The first visit was performed on 23 December 2009, which covers 86% of its 3.92 h rotational period followed by a shorter follow-on observation on 25 December 2009. The second one was obtained on 20 June 2010 (follow-on observations on 21 June 2010) using the same detector and bands and covering 110% of its rotational period.

The Plutino (84922) 2003 VS₂ was observed with Herschel/PACS in mini scan-maps mode at 70 and 160 μm on 10 August 2010, covering 106% of its 7.42 h rotational period. Follow-on observations were performed on 11 August 2010.

The binary Plutino (208996) 2003 AZ₈₄ was observed with Herschel/PACS in mini scan maps mode at 100 and 160 μm on 26−27 September 2010 (with follow-on observations on 28 September 2010). The observation lasted ~7.4 h, i.e 110% of an assumed single-peak rotational period of 6.79 h (or 55% of a double-peaked period of 13.58 h).

All observations were made using only one scanning direction (see e.g. Santos-Sanz et al.2012, for a detailed description of the mini scan maps mode in the case of TNOs²).

Table 1 shows the orbital parameters, B − R colors, absolute magnitudes, rotational properties, taxonomy, and dynamical classification of the observed objects. This table also includes the radiometric solutions of these three objects (equivalent diameter for an equal-area sphere, geometric albedo, beaming factor) previously published as part of the “TNOs are Cool” project. Table 2 shows the observational circumstances of each one of these TNOs.

Owing to the spatial resolution of Herschel/PACS, the satellites of the binaries/multiple systems (i.e. Haumea and 2003 AZ₈₄) are not resolved, and their thermal fluxes are merged with the thermal flux of the main body.

3. Data reduction and photometry

3.1. Data reduction

PACS images obtained from the Herschel Space Observatory were processed using the Herschel Interactive Processing Environment (HIPE³) and our own adapted pipelines developed within the “TNOs are Cool” project. The application of the pipeline provides individual or single maps, each one of these single images covering ~4.7 min. The re-sampled pixel scale of the single maps is 1.1″/pixel, 1.4″/pixel, and 2.1″/pixel for the 70 μm (blue), 100 μm (green), and 160 μm (red) bands, respectively. Apparent motion over the duration of an individual map is negligible compared to the PACS PSF (FWHM in radius is 5.2″/7.7″/12″ in blue/green/red bands, respectively) and does not need to be corrected. These single maps are combined afterwards, using ephemeris-based recentering processes within HIPE, to obtain enough signal-to-noise ratio (S/N) to perform a good photometry, while keeping enough time resolution to resolve the thermal light curve, in a similar way to Lellouch et al. (2010). The exact number of individual and combined maps for each target is detailed in Sect. 3.2.

To minimize at best contamination by background sources, all light curve data are associated with complementary observations acquired one or few days later (follow-on observations), where the target has moved enough that a background map can be determined and subtracted from each combined map of the light curve. The method was demonstrated in Spitzer/MIPS TNOs/Centaurs observations (Stansberry et al. 2008; Brucker et al. 2009). However, this technique to remove background sources fails when trying to remove some background features in the 2003 AZ₈₄ images. In this case another technique, known as double-differential background subtraction, is applied. A complete and detailed description of the data reduction process, the background subtraction and the double-differential techniques applied to the Herschel/PACS images can be found in Santos-Sanz et al. (2012) and in Kiss et al. (2014). Figures 1, 7, and 10 illustrate the advantages of these background-removing techniques for the three targets respectively.

3.2. Photometry

As indicated above, single maps obtained with a time resolution of 4.7 min were combined to improve image quality for photometry. Typically, the number of single images to be combined was larger at 160 μm than at 70/100 μm, owing to lower S/N and larger sky residuals. The details can be found below:

Haumea data were taken in two epochs, each time using the green (100 μm)/red (160 μm) filter combination. For the first (resp. second) epoch, the total number of single images is 40 (resp. 55) per filter. These images were grouped by 4 (18.8 min time resolution) in the green and by 6 (28.2 min per data point) in the red. The choice of this particular grouping of single images for the green (by 4) and the red (by 6) for the Haumea data is based on a compromise between obtaining enough S/N to extract a reliable photometry and having enough time resolution to properly sample the light curve, as mentioned above. Usually, more single images must be grouped for the red band because those images are normally noisier than images at other shorter wavelengths (even after the application of background-removing techniques).

The different grouping elections for 2003 VS₂ and 2003 AZ₈₄ are based on the same described compromise between S/N and time resolution. After removing clear outliers in the Haumea data, 37 images remain for the first epoch (resp. 53 for the second epoch) at green band, and 35 images for the first epoch (resp. 50) at red band (see Table A.1).

2003 VS₂ was observed in blue/red combination, with 96 single maps for each color. We grouped these single images by 5 for the blue (23.5 min per data point) and by 10 for the red band (47.0 min per data point). Each blue data point is independent of the previous and following point (no time overlap) while, for the red, each data point has a time overlap of 23.5 min with the previous and following point. Consecutive points for both bands have a separation of ~0.05 in rotational phase, clear data outliers have been removed. At the end, 18 data points remain for the blue band, and 18 for the red one (see Table A.2).

Similarly, 95 single maps of 2003 AZ₈₄ were acquired in green/red combination. They were grouped by 6 in the green without time overlap between previous and following point, and with a separation between consecutive points of ~0.07 in phase. Final data points for the green are 15 (see Table A.3). Red images were discarded for a thermal light curve analysis because they are very contaminated by background sources (even after the application of background removing techniques) and the final photometry on these images is very noisy. We still use the mean value of the red band flux (see Table 4) for thermophysical model analysis.

Photometry was performed on the combined maps. The flux of the objects was obtained using DAOPHOT (Stetson 1987) routines adapted to IDL⁴ to perform synthetic-aperture photometry on the final images. The object is usually located at the centre or very close to the centre of the images. Since our targets are bright enough, we do not need to use ephemeris coordinates to find them: photocentre routines are then used to obtain the best coordinates to place the centre of the circular aperture. Once the photocentre is obtained, we performed aperture photometry for radii that span from 1 to 15 pixels. We applied the aperture correction method (Howell 1989) for each aperture radius using the tabulated encircled energy fraction for a point-source observed with PACS⁵. Uncertainties on the fluxes are estimated by means of a Monte-Carlo technique, in which artificial sources are randomly implanted on the images. We obtain and correct by aperture the fluxes of these artificial sources using a median optimum aperture radius. Uncertainties are computed as the standard deviation of these implanted fluxes. These photometric techniques and uncertainty estimations used to extract the PACS fluxes are further described in Santos-Sanz et al. (2012) and Kiss et al. (2014).

In addition to the random photometric errors, the data may suffer from a systematic flux calibration uncertainty, which is estimated to be ~5%. As the latter affects all points of a given light curve in an identical way, it is not included in the individual error bars. Color corrections, which are ~1−2% for Herschel/PACS data, were applied to the fluxes (see caption in Tables A.1−A.3 for the exact value of the color-correction factors applied).

Finally, time-phasing of all the images was computed using the preferred rotational periods (see Table 1) of the objects. For Haumea, the relative phasing of data taken at two epochs separated by six months did not pose any problem, thanks to the highly accurate knowledge of the period⁶. A running mean of the Haumea thermal fluxes in each filter was applied with a bin of 0.05 in rotational phase (=11.75 min of time), which finally obtained 20 points at green band and 19 at red band. The final thermal light curves are shown in Figs. 2, 8, and 11 for Haumea, 2003 VS₂ and 2003 AZ₈₄, respectively.

4. Results and analysis

4.1. (136108) Haumea

Fig. 1 Images (top line: 100 μm; bottom line: 160 μm) extracted from the 20 Jun. 2010 light curve of Haumea. Left: original images. Right: background-subtracted images. The Field of View (FOV) is 2.5′× 2.5′. Haumea is at centre.

Haumea’s optical light curve is one of the best studied among TNOs (Rabinowitz et al. 2006; Lacerda & Jewitt 2007; Lacerda et al. 2008; Thirouin et al. 2010). Besides its strong amplitude (0.28 mag), its most remarkable feature is its asymmetric character, exhibiting two unequal brightness maxima, which cannot be explained by a pure shape effect and is interpreted as being due to the presence of a darker (and redder) region on the object’s surface (Lacerda et al. 2008).

Figure 2 shows the thermal data for Haumea at 100 and 160 μm, respectively, as a function of rotational phase, using a period of 3.915341 h. The zero phase epoch is JD = 2 455 188.720000 (uncorrected for light-time) and phases are calculated using light-time corrected Julian dates and light-time corrected zero date. In addition to the thermal fluxes, Fig. 2 displays the Scaled optical LC, which represents the optical fluxes rescaled by some multiplicative factor to match the mean thermal flux level. The part of the optical light curve that is affected by the dark spot is outlined in yellow. The overall positive correlation between the thermal and optical light curves indicate that both are mostly caused by the object elongated shape, as already noted in Lellouch et al. (2010, hereafter Paper I). However, the highest quality 100 μm data further indicates an asymmetry in the two thermal flux maxima, with the strongest occurring near phase ~0.75, i.e. in the part of the optical light curve that is affected by the dark spot. The possible influence of the spot could not be discerned in Paper I, and the present data are of higher quality. On the other hand, the 160 μm data do not show evidence for an enhanced thermal flux associated with the dark spot.

A Fourier fit of the thermal data permits us to determine the amplitude of the thermal light curve (defined as the difference between maximum and minimum fluxes in the Fourier fit), as well as its phasing relative to the optical light curve. In both filters, the thermal light curve amplitude is larger than its optical counterpart and diminishes with increasing wavelength; these behaviors are in accordance with thermophysical model expectations. However, the two thermal filters do not give fully consistent information of the phase shift between the thermal and optical data: 100 μm data appear well in phase with the optical light curve, while 160 μm data appear shifted by 0.06 in phase (i.e. by about 21 degrees: see Table 3).

Finally, we perform a consistency check of the fluxes obtained in the thermal light curves. To do this, we run a Near Earth Asteroid Thermal Model (NEATM, Harris 1998) for the green and red fluxes at the minimum and maximum of the thermal light curves. We assume H_V = 0.43 ± 0.01 mag as in Fornasier et al. (2013) and the geometric albedo and beaming factor derived in that work (p_V = 80.4%, η = 0.95). Under these assumptions we estimate the area-equivalent diameter from the NEATM for the minimum of the thermal light curves, using an absolute magnitude of (H_V + 0.28 / 2) mag, where 0.21/2 is the semi-amplitude from the optical light curve, obtaining D_min = 1173 km. Running a NEATM in the same way for the maximum, for an absolute magnitude (H_V − 0.21 / 2) mag, we obtain D_max = 1292 km. These diameters are consistent, within error bars, with the best equivalent diameter obtained in Fornasier et al. (2013).

Fig. 2 Haumea’s thermal light curve at 100 μm (top) and 160 μm (bottom), combining data from the two visits (see text for details). The black dashed curve is the “scaled optical light curve” obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The part of this curve that is affected by Haumea dark spot, according to the Lacerda et al. (2008) preferred model, is outlined in yellow. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 455 188.720000 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.1 for further details). The uncertainty in the rotational phase is ±0.001.

Following our work in Lellouch et al. (2010), modeling of the thermal light curve was performed using OASIS (Optimized Astrophysical Simulator for Imaging Systems: Jorda et al. 2010). OASIS is a versatile tool in which an object is described by triangular facets. The orientation of each facet with respect to pole orientation, Sun direction, observer direction, and time as the object rotates, is calculated. OASIS therefore requires a shape model and an assumed aspect angle (i.e., pole orientation). For the shape of Haumea, we used an ellipsoid made of 5120 triangles. For the aspect angle, the large amplitude of Haumea’s optical light curve favors a large angle and, here, we assumed an equator-on geometry (aspect angle θ = 90°). A spectral and bolometric emissivity of 0.9 in all filters is assumed as well.

In Paper I, two shape models for Haumea (defined by the a, b, and c semi-major axes of the ellipsoid) were used, based on optical light curves observations by Lacerda et al. (2008) and assuming a Jacobi hydrostatic equilibrium figure (a>b>c). The two shape models were derived by considering two different scattering properties for Haumea’s surface (Lambertian reflectivity, shape model 1, and Lommel-Seelinger reflectance properties, shape model 2), leading to slightly different values of b/a, c/a, and the object density. For a given shape model, knowledge of the object mass (Ragozzine & Brown 2009) provided the absolute values of a, b, and c, which in turn provided the object mean geometric albedo, based on its H_v magnitude. All these parameters were then implemented in a NEATM thermal model (Harris 1998), and the only free parameter in fitting the thermal light curve was the so-called beaming factor, η. In this process, a phase integral q = 0.7 was adopted; this value is reasonable for a high albedo object (see Lellouch et al.2000, Fig. 7; and Brucker et al.2000) but admittedly uncertain. Considering mostly an aspect angle θ = 90°, the main conclusion of Paper I was that η = 1.15 satisfies the mean thermal flux constraint for both shape models, but matches the light curve amplitude only for model 2, which was therefore favored. The relatively low η value (for an object at this distance from the Sun) pointed to a generally low thermal inertia for the surface and significant surface roughness effects (see Lellouch et al. 2013). Paper I also briefly explored the effect of aspect angle by considering the case θ = 75°, and found that such a model could be valid, but using a slightly larger η value (e.g. η = 1.35 instead of η = 1.15). However, the modest quality of the thermal light curve in Paper I did not warrant the use of more elaborate models.

Given the improved data quality in the current work, including the apparent detection of increased thermal emission at the expected location of the dark spot, we now improve these early models by (i) considering thermophysical models (TPM); (ii) exploring in some detail the effect of a surface spot. The essential physical parameter to constrain is now the surface thermal inertia, Γ. To make allowances for possible surface roughness effects, however, the TPM can also include an η factor, but which in this formalism is by definition ≤1 (see e.g. Groussin et al.2004, Eq. (3); Lellouch et al.2011, Eq. (2)). Unlike NEATM, which for a uniform (constant albedo) elliptical surface, calculates (by construction) a thermal light curve in phase with the optical light curve, the TPM approach enables us to investigate temperature lags owing to thermal inertia. Thus, in principle, the thermal inertia can be derived by investigating the relative phase of the thermal and optical light curves, as constrained by the observations. Once Γ is determined, η and q may be adjusted so as to match the mean flux level and the amplitude of the light curve. However, the problem may be underconstrained, i.e. η and q cannot necessarily be determined separately. If the surface includes a spot of known albedo and spatial extent, the parameters (i.e. Γ, and η and/or q) may be adjusted separately in the spot region and outside. In what follows, and given the low thermal inertias we inferred (see below), we found that the observed flux levels did not require to be enhanced by surface roughness, so we simply assumed η = 1, recognizing that some degeneracy exists between η and q.

Lacerda et al. (2008) show that the asymmetry in Haumea’s optical light curve can be interpreted with different spot models, characterized by the albedo contrast of the spot with respect of its surroundings and its spatial extent. In all cases, the spot is assumed to be centred on Haumea’s equator and to lead one of the semi-major axes by 45°. Possible models range from a very localized (6% of Haumea maximum cross section) and low albedo spot (30% of the non-spot albedo) to much more extended spot (hemispheric) and subdued in contrast (95% of non-spot albedo). In a brief study of the spot effect on the thermal light curve, Paper I considered a spot covering 1/4 of Haumea’s maximum projected cross section, with an albedo contrast (about 80%) of the non-spot albedo, as prescribed by the Lacerda et al. (2008) results. The same spot description was adopted here. Being relatively limited in extent, the spot has negligible effect on about half of the thermal light curve (from phase ~0.0 to ~0.5 with the adopted phase convention). Therefore, as a first step, we focus on the part of the light curve that is not affected by the spot.

Fig. 3 Comparison of the Haumea 100 μm light curve with homogeneous models (no spot) having various values of the thermal inertia Γ in the range 0.0 to 0.5 MKS. No surface roughness effects are included (i.e. η = 1). For each Γ value, the phase integral is adjusted to match the mean flux level and light curve amplitude. The lowest Γ values (0.0−0.2 MKS) provide the best fit to the data.

Figure 3 shows the comparison of the observed 100 μm light curve with several homogeneous (no spot) models differing by their surface thermal inertia (Γ = 0.0 to 0.5 MKS by steps of 0.1). Here, and throughout the following, shape model 2 is adopted (this shape model is also favored by the analysis of Lockwood et al. 2014) with an aspect angle θ = 90°, giving a = 961 km, b = 768 km, c = 499 km, and p_v = 0.71. These a, b, c values lead to a mean area-equivalent diameter D_equiv = 2 a^{1 / 4}b^{1 / 4}c^{1 / 2} = 1309 km, within the error bars of the area-equivalent diameters obtained from radiometric techniques for this object (1324 ± 167 km from Lellouch et al.2010; 1240$\begin{array}{c}\mathrm{+}\mathrm{69}\\ -58\end{array}$ km from Fornasier et al.2013). Using a shape model that is consistent with the optical light curve is preferable to a radiometric solution that may include measurements at different light curve phases. For each value of Γ, the phase integral (q) is adjusted to provide the best fit to the mean flux level and amplitude of the light curve outside the spot region (i.e. at phases 0.0–0.5). For Γ = 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5 MKS, respectively, the required values of q are 0.851, 0.778, 0.712, 0.660, 0.613, and 0.578, respectively. As can be seen in Fig. 3, Γ = 0.0 provides the best fit to the part of the light curve not affected by the spot, while larger Γ values progressively lead to larger delays of the thermal emission, which are not observed in the 100 μm data. Detailed comparisons of the data show that the reduced χ² is minimum for Γ = 0.0 MKS and larger but still reasonable for Γ = 0.2 MKS. From this we conclude that Γ values in the range 0.0–0.2 MKS are consistent with the data. The associated phase integral values are in the range 0.851–0.712. These are generally consistent but still somewhat larger than the 0.7 value assumed in Paper I. We note that using the Brucker et al. (2009) empirical relationship between geometric albedo and phase integral, phase integrals of 0.851 (obtained for Γ = 0 MKS) – 0.712 (for Γ = 0.2 MKS) would imply geometric albedos in the range 1.05–0.64. The latter value is more consistent with p_v = 0.71, as indicated by shape model 2. From this point of view, Γ = 0.1 − 0.2 MKS would seem a more plausible solution but, as indicated above, the associated fits are worse than with Γ = 0 MKS. We also note that these q values hold for our assumption of η = 1 (no surface roughness), and that for a given Γ, any finite surface roughness would require an even larger value of q.

Fig. 4 Models with variable phase integral in the spot region. In these models, the spot thermal inertia is Γspot = 0, as in the no-spot region. The red curve (having q = 1.00 in the spot) represents the first best fit model for the 100 μm data.

We now turn to constraints on the dark spot. Continuing with the assumption of no surface roughness, we investigate the effect of a specific phase integral or thermal inertia in the spot. In all models, the spot has an equivalent geometric albedo of 79% of its surroundings (p_v = 0.71), i.e. p_v = 0.56. Figure 4 shows the effect of changing the phase integral in the spot (q_spot) by steps of 0.05, maintaining its thermal inertia to Γ_spot = 0 as in the best fit solution of the no-spot region (Fig. 3). As is clear from Fig. 4, mild changes in q_spot, e.g. from ~0.85 to ~1.05 produce flux variations at the adequate level. The flux excess near phase 0.75 would point to q_spot = 0.85, essentially identical to the non-spot region (meaning that the flux excess is purely an effect of the spot darker albedo). However this model (blue curve in Fig. 4) clearly overpredicts the fluxes beyond the peak phase, where data are best fit with q_spot = 0.95–1.05. The best overall fit of the light curve using an χ² criterion is achieved with q_spot = 1.00, and the model light curve shows only a very weak flux excess in the spot region (meaning that the effect of the darker spot albedo is essentially compensated for by the effect of a larger phase integral). A phase integral of 1.00 seems very high (and not in line with the Brucker et al.2009, relationship), but it has been observed in other solar system objects (see Hillier et al. 1991, for Triton). However, while the fit is satisfactory (first best-fit model in Fig. 4), we note that it might not be a very realistic solution, since there is a general positive correlation between albedo and phase integral on airless bodies (Lellouch et al.2000, Fig. 7; Brucker et al.2009), in contradiction with the dark (low albedo) spot, which has a larger phase integral (higher albedo) than its surroundings.

Fig. 5 Models with variable thermal inertia in the spot region. In these models, the spot phase integral is kept at the same value as in the no-spot region. The green, black, and blue curves are for Γ = 0.1 MKS and q = 0.778 outside of the spot, and Γ = 0.0, 0.1, and 0.2 MKS in the spot. The red curve represents the second best fit model at 100 μm with Γ = 0.0 MKS and q = 0.851 outside of the spot, and Γspot = 0.05 MKS in the spot.

Figure 5 shows the effect of changing the spot thermal inertia, maintaining its phase integral to the value of its surroundings. For this exercise, we first consider Γ = 0.1 MKS in the no-spot region, associated with q = 0.778. We maintain the q_spot at this value and show models with Γ_spot = 0.0, 0.1, and 0.2 MKS. These models show that variations of Γ within this range produce flux changes at the appropriate level. In this family of models, the best fit is achieved with Γ_spot = 0.2 MKS (blue curve in the figure), which again suggest that, to avoid too large fluxes past the peak, the thermal effect of darker albedo needs to be partly compensated for by a larger thermal inertia. As discussed above, however, Γ = 0.1 MKS is not an optimum fit of the no-spot part of the light curve. The red curve in Fig. 5 shows a model with Γ = 0.0 MKS and q = 0.851 outside of the spot, and Γ_spot = 0.05 MKS (and still q_spot = 0.851) in the spot (second best fit model). This model fit is worse than the first best fit model (Γ_spot = 0.0 MKS, q_spot = 1.00) and, in fact, no better than the no-spot model: the fit improves in the region of flux maximum around phase 0.75, but the non-zero thermal inertia in the spot tends to delay the model too much in the 0.85–1.0 phase region. Thus the models of Fig. 5 indicate that, at most, the thermal inertia in the spot region is slightly higher than in its surroundings. This type of behaviour might be expected. A statistical study of TNO thermal properties (Lellouch et al. 2013) suggests that the highest albedo TNOs generally exhibit particularly low thermal inertia. A plausible explanation is that in addition to a specific composition (e.g. pure ice), a higher albedo may reflect a smaller grain size. In fact, the generally low thermal inertia of TNOs points to a regime where radiative conductivity (i.e. in surface pores) is important in the overall heat transfer, thus high albedo objects might plausibly be associated with low thermal inertia, and the association of lower albedo with larger thermal inertia is also plausible.

Fig. 6 The two best fit models from the 100 μm light curve (see Figs. 4 and 5) applied to the 160 μm light curve (red and blue curves, respectively). The green curve is the best-fit model in terms of least squares adjusted to match the 160 μm phase, having Γ = 0.5 MKS and q = 0.68 throughout Haumea’s surface.

Figure 6 shows the application of the above two best fit models to the Haumea’s 160 μm light curve. In both cases, the agreement with the observed 160 μm mean flux, and especially light curve amplitude is reasonable, but the models fail to match the data in two resspects: (i) the modeled flux levels are on average too high (by ~10%); (ii) the models are somewhat out of phase with the observations. The first problem may suggest a calibration error in the 160 μm data (more subject to sky contamination). The second one is related to the fact that the 160 μm data appear to be shifted by about 0.06 ± 0.01 in rotational phase (~21 degrees or ~13 min of time, see Table 3) with respect to the 100 μm data. Optimum phasing of the model with respect to the 160 μm data would require a thermal inertia Γ~ 0.5 MKS, and fitting the mean flux levels would then require a phase integral q = 0.68. This model tailored to the 160 μm data does not include any specific values of Γ or q in the dark spot region, as this is not required by these data. Overall a simultaneous fit of the 100 and 160 μm data within error bars does not appear possible. Although so-called compromise parameters could be formally found by performing an χ² minimization on both datasets simultaneously, we feel that separate modeling is preferable since it provides a handle on model limitations and realistic range of solution parameters. Giving more weight to the higher quality 100 μm data, we favor Γ = 0.0–0.2 MKS but, based on 160 μm modeling, we regard Γ = 0.5 MKS as an upper limit to Haumea’s thermal inertia.

In summary, thermophysical modeling of the Haumea light curves indicates that: (i) the object’s thermal inertia Γ is extremely small (less than 0.5 MKS and probably less than 0.2 MKS) and its phase integral is high (at least 0.8 for Γ = 0.1 MKS and probably even higher if surface roughness is important); (ii) only small changes in the surface properties of the dark spot (e.g. changes in the thermal inertia by ~0.1 MKS or in the phase integral by ~0.1) are required to significantly affect the emitted fluxes on the hemisphere where the spot resides; larger changes are excluded by the data; (iii) the most plausible scenario may invoke a slightly higher thermal inertia in the dark spot compared to its surroundings, but a fully consistent picture is still not found, since the ~21 degrees shift in phase between the 100 and 160 μm data is difficult to understand.

Finally, we note that in all of our Haumea models we ignored the possible contribution of its satellites Hi’iaka and Namaka. Although their albedos are unknown, these moons are thought to have been formed by a catastrophic impact that excavated them from the proto-Haumea ice mantle and led to the Haumea family (Brown et al. 2007) or from rotational fission (Ortiz et al. 2012). As such, their albedos are probably comparable to, or even higher than, Haumea’s itself. Assuming a 0.70 geometric albedo, Hi’iaka and Namaka diameters are ~320 and ~160 km, respectively (Thirouin et al. 2016). Furthermore, assuming identical thermophysical properties, they would contribute in proportion of their projected surfaces, i.e. 6% and 1.5% of Haumea’s thermal flux. Although this is not negligible, we did not include this contribution owing to its uncertain character. Should the Haumea’s thermal fluxes to be modeled decrease by ~20 mJy × 7.5% ~1.5 mJy, this could be taken care of in the models by a slight increase of the phase integral for a given thermal inertia, without any changes to the conclusions on the object thermal inertia and the dark spot properties.

4.2. (84922) 2003 VS₂

2003 VS₂ is a Plutino without known satellites. Near infrared spectra of this body shows the presence of exposed water ice (Barkume et al. 2008), which probably increases the geometric albedo of this Plutino (~15%, according with Mommert et al. 2012), compared with the mean albedo of TNOs without water ice. This object presents an optical light curve with moderately large peak-to-peak amplitude ~0.21 ± 0.01 mag and a double-peaked rotational period ~7.42 h (Ortiz et al. 2006; Sheppard 2007; Thirouin et al. 2010). To fold the Herschel data to the rotation period with enough precision, we refined the knowledge of the rotational period. To achieve this goal, we used optical observations taken on 4−8 September 2010 by means of the 1.5-m telescope at Sierra Nevada Observatory (OSN, Spain) using a 2k × 2k CCD with a FOV of 7.8′× 7.8′ and 2 × 2 binning mode (image scale = 0.46″/pixel). Then we merged these observations with old optical light curves obtained also at OSN on 22, 26, 28 December 2003 and 4, 19−22 January 2004, and with 2003 observations published in Sheppard (2007). No filter was used to perform the OSN observations. We reduced and analysed all these data as described in Thirouin et al. (2010) to finally obtain an accurate rotational period of 7.4175285 ± 0.00001 h (Δm = 0.21 ± 0.01 mag). A further description of the observations and techniques leading to this very accurate rotational period are detailed in Thirouin (2013).

Fig. 7 Images (top line: 70 μm; bottom line: 160 μm) extracted from the thermal time series curve of 2003 VS2. Left: original images. Right: background-subtracted images. The FOV is 2.5′× 2.5′. 2003 VS2 is at centre.

The thermal light curves at 70 μm and 160 μm are not firmly detected with only a 1.7σ and 1.5σ confidence levels respectively (see Table 3 and Fig. 8). A Fourier fit of the 70 μm data indicates a mean flux of 14.16 mJy and an amplitude of 1.73 mJy, which is slightly smaller than the optical light curve amplitude. The same analysis yields a negligible shift in time of the 70 μm data relative to the optical (Fig. 8). The lower quality 160 μm data would suggest a − 0.054 ± 0.102 phase shift (see Table 3). We do not consider this last negative shift significant since it is well below the estimated error bars. In what follows, we pursue with our thermophysical modeling, focusing on the mean thermal flux and light curve amplitude at 70 μm.

Fig. 8 Thermal time series curve at 70 μm (top) and 160 μm (bottom) for 2003 VS2. The black dashed curve is the scaled optical light curve obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 452 992.768380 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.2 for further details). The uncertainty in the rotational phase is ±0.01.

The large amplitude of the optical light curve and its double peaked nature indicates that the main cause for the variability is a triaxial shape (Sheppard 2007; Thirouin et al. 2010). Then, if we assume that the optical light curve is entirely shape-driven, its period and amplitude can be used to derive a shape model under the assumption of hydrostatic equilibrium. As for Haumea, the large object size (~500 km) makes the assumption of a Jacobi hydrostatic equilibrium figure (semi-major axes a>b>c) reasonable. Further assuming an equator-on viewing geometry (aspect angle = 90°), the Δm = 0.21 mag light curve amplitude is related to shape by $\mathrm{\Delta }\mathit{m}\mathrm{=}\mathrm{2.5}\hspace{0.17em}\log \left(\frac{\mathit{a}}{\mathit{b}}\right)\mathrm{·}$(1)Using the rotation period of 7.42 h and the Chandrasekhar figures of equilibrium tables for the Jacobi ellipsoids (Chandrasekhar 1987), we obtain b/a = 0.82, c/a = 0.53 and ρ = 716kg / m³, where ρ is a lower limit to the density because, if the object is not observed equator-on, the true a, b axial ratios may be higher and the implied density would be also higher. Using the area-equivalent radiometric diameter (D_equiv = 523 km) derived from earlier thermal modeling of Herschel and Spitzer data (Mommert et al. 2012; Lellouch et al. 2013), defined as D_equiv = 2 a^{1 / 4}b^{1 / 4}c^{1 / 2}, we obtain the values of the semi-major axes of 2003 VS₂: a = 377 km, b = 310 km, and c = 200 km. We further adopt a phase integral q = 0.53 and a V geometric albedo p_V = 0.147 from the previous papers and, as for Haumea, we do not consider surface roughness (i.e. η = 1). All these values are used as input parameters to the OASIS code for modeling the 70 μm thermal light curve of 2003 VS₂. With this approach, the only free parameter is the thermal inertia Γ. Figure 9 shows model results for Γ = 1.0, 2.0 and 3.0 MKS. In this figure, the phase of the thermal models is determined by requiring that the model with Γ = 1.0 matches the observed phase of the thermal data (i.e. a maximum at phase 0.55, see top panel of Fig. 8). We note that, unlike in the Haumea case (see Fig. 3), models with the various thermal inertias all appear to be approximately in phase. The difference in behaviour is caused by the combination of the longest period and warmer temperatures at 2003 VS₂ versus Haumea. For a given thermal inertia, this causes a much smaller value of the thermal parameter for 2003 VS₂ (see e.g. Lellouch et al. 2013). A thermal inertia of 2.0 MKS matches reasonably well the mean flux levels and the light curve amplitude, while the other two models significantly over- or underestimate the mean flux. Thus, we conclude to a thermal inertia Γ = (2.0 ± 0.5) MKS for 2003 VS₂. This value is fully consistent with the mean thermal inertia for TNOs and centaurs derived statistically from the Herschel “TNOs are Cool” sample (Γ = 2.5 ± 0.5 MKS, Lellouch et al. 2013), but significantly above that for Haumea. As indicated in the above paper, high-albedo objects seem to have preferentially low thermal inertias.

Fig. 9 Thermophysical model fits of the 2003 VS2 70 μm data for various values of the thermal inertia. From top to bottom: Γ = 1.0, 2.0 and 3.0 MKS. Shape model and other physical input parameters are given in the text.

4.3. (208996) 2003 AZ₈₄

2003 AZ₈₄ is a binary Plutino with a shallow optical light curve (Δm ~ 0.07 mag) and with a rotational period ~6.79 h, assuming a single-peaked light curve (Sheppard & Jewitt 2003; Ortiz et al. 2006; Thirouin et al. 2010). Nonetheless, the double-peaked solution, which corresponds to a rotational period ~13.58 h, cannot be totally discarded (Thirouin et al. 2010). Near infrared spectra have detected water ice on its surface (Barkume et al. 2008; Guilbert et al. 2009; Barucci et al. 2011).

Fig. 10 Left panel: one of the 2003 AZ84 images at 100 μm (green) with the original background. Right panel is the same map after the application of the double-differential background subtraction technique, as described in Kiss et al. (2014). 2003 AZ84 is the source at the center of the images. In the right panel, the black source at the bottom left of 2003 AZ84 is the negative image of 2003 AZ84 resulting from the application of this background-removing technique. The FOV of each image is 2.5′× 2.5′.

As for the other two TNOs, we acquired additional time series images of 2003 AZ₈₄ on 4−5 February 2011 with the 1.23-m telescope at Calar Alto Observatory (CAHA, Spain), equipped with a 2k × 2k CCD camera in the R filter. These data are merged with old CAHA and OSN data from 2003 and 2004 to refine the rotational period, obtaining P = 6.7874 ± 0.0002 h (Δm = 0.07 ± 0.01 mag), which we nominally use to fold the Herschel/PACS data and compare them with the visible light curve. A more detailed description of the observations and techniques of analysis leading to this rotation period are included in Thirouin (2013).

The thermal light curve of 2003 AZ₈₄ is not firmly detected in the PACS data at 100 μm (see Fig. 11). While a Fourier fit to the thermal data formally provides a best fit amplitude of 1.97 ± 1.40 mJy at 100 μm, its significance is thus at the 1.4σ level. By making a visual comparison of the Fourier fit to the thermal data in Fig. 11 with the Fourier fit to the optical data (shown with a dashed line), it appears that there could be a weak anticorrelation of the 100 μm data with the visible data. This would give confidence to the interpretation that the thermal light curve could be generated by enhanced thermal emission in the darker spots or darker terrains that give rise to the optical light curve, in the same way as the dark spot in Haumea generates enhanced thermal emission. However a Spearman test to analyze a possible anticorrelation of the thermal data with the optical light curve gave a non-significant result. Moreover, in the regime of low albedo (~10% for 2003 AZ₈₄), the thermal emission is essentially albedo-independent, so that a thermal light curve resulting from optical markings would have an undetectable amplitude, barely above 0.05 mJy.

Fig. 11 Thermal time series curve at 100 μm for 2003 AZ84. The black dashed curve is the scaled optical light curve obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 453 026.546400 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.3 for further details). The uncertainty in the rotational phase is ±0.01.

4.3.1. Analysis of 2003 AZ₈₄ results

2003 AZ₈₄ is a large enough TNO (D_equiv = 727 km according Mommert et al. 2012) so that it is very likely to be in hydrostatic equilibrium. This means that the expected 3D shape of this TNO should be a figure of equilibrium: either a rotationally symmetric Maclaurin spheroid or a triaxial Jacobi body. 2003 AZ₈₄ has a low light curve amplitude in the visible, which means that either this TNO is a Maclaurin object with small albedo variability on its surface, or it is seen nearly pole on, or both. Recent results on stellar occultations by 2003 AZ₈₄ (Dias-Oliveira et al. 2016) have shown an equal-area diameter of D_equiv = 766 km and a small projected flattening of only 0.05. The small flattening has two possible extreme explanations: the object has a typical density of a TNO of its size but it is seen nearly pole-on so that the large flattening of the body becomes a small projected flattening, or 2003 AZ₈₄ has an exceptionally high density for its size. The density required for a Maclaurin body with flattening of 0.05 and a rotation period of 6.79 h is 5500 kg/m³. This huge density is not feasible in the trans-Neptunian region. Assuming that 2003 AZ₈₄ could have a density of ~2500 kg/m³, which is already too high for a TNO of its size, its true oblateness would be 0.12. This would require an aspect angle <45 degrees for the Maclaurin spheroid to give rise to the projected oblateness of 0.05 seen in the occultation. We note that densities of around 2500 kg/m³ have only been measured for the very largest TNOs, such as Pluto, Eris and Haumea whose internal pressures do not allow for the macroporosity that can exist in bodies of smaller size (see e.g. Jewitt & Sheppard 2002). So it is extremely unlikely that 2003 AZ₈₄ could have such a high density of 2 500 kg/m³. Hence we are confident that the aspect angle of 2003 AZ₈₄ must be smaller than 45 degrees. Therefore the low light curve amplitude in the visible, the small thermal variability, and the occultation results are reasons to believe that 2003 AZ₈₄ could be close to pole-on (have a small aspect angle).

To further constrain the spin axis orientation, we run a thermophysical model (TPM: Lagerros 1996, 1997, 1998). The model takes into account the thermal conduction and surface roughness for objects of arbitrary shapes and spin properties. The model was extensively tested against thermal observations of near Earth asteroids (NEAs), Main Belt asteroids (MBAs), and TNOs over the last two decades. Recent works with this TPM code (e.g. Pál et al. 2015, 2016; Schindler et al. 2017; Fornasier et al. 2013) have shown that a combined analysis of Herschel and Spitzer thermal measurements enable us to constrain the spin-axis orientation of TNOs. Following up on this expertise, we also applied this code in the case of 2003 AZ₈₄, assuming a spherical shape model and a constant emissivity of 0.9 at all MIPS and PACS wavelengths (emissivity effects are only expected at longer wavelengths, see Fornasier et al. 2013). We used all available Herschel/PACS and Spitzer/MIPS thermal data (except the Spitzer/MIPS data point at 71.42 μm, which is affected by a background source). The Spitzer MIPS data were presented in Stansberry et al. (2008). We re-analysed the two MIPS observations of 2003 AZ₈₄. The 71.42 μm measurements are problematic owing to contaminating background sources, but the 23.68 μm points are clean and the object’s point-spread-function is as expected. Table 4 shows the thermal measurements used in the thermophysical modeling. Using H_V = 3.78 ± 0.05 mag and the stellar occultation size D = 766 ± 16 km (Dias-Oliveira et al. 2016), we check models with a range of thermal inertias from 0.0 to 100 MKS and a range of different levels of surface roughness (rms of surface slopes of 0.1, 0.3, 0.5, 0.7, and 0.9). The biggest issue with the radiometric analysis is that the MIPS and the PACS data do not match very well. A standard D − p_V radiometric analysis for only the PACS data favors the pole-on solutions (combined with low surface roughness − rms of surface slopes <0.5) and provide the correct occultation size (760−790 km), for the two possible rotation periods (P = 6.7874 h and P = 13.5748 h) and almost independent of thermal inertia. The overall fit to only PACS observations is excellent. Equator-on solutions can also provide correct sizes, but only under the assumptions of extremely low thermal inertias far below 1.0 MKS and combined with extremely high surface roughness (rms of surface slopes >0.7), which looks very unrealistic. The very best solutions are found for a spin-axis orientation 30° away from pole-on, intermediate levels of surface roughness (i.e. realistics values), and acceptable values for the thermal inertia (Γ = 0.5 − 3.0 MKS). Overall, the pole-on ±30° configuration explains very well the PACS fluxes, but slightly overestimates the MIPS 24 μm within the 2σ level (see Fig. 12). This difference between MIPS flux and model could be due to some light curve effects at the moment of the Spitzer/MIPS observations. The equator-on geometry only works when using very extreme settings, which seem very unrealistic. Summarizing the combined thermal and occultation analysis, we find that its spin-axis is very likely close to pole-on (±30°).

Fig. 12 Absolute PACS and MIPS fluxes for 2003 AZ84 with various thermophysical models overplotted: pole-on solution, pole-on +30°, pole-on +60°, equator-on. All the models use the occultation size and reasonable assumptions for the thermal properties.

5. Summary and brief discussion

Time series thermal data of three bright TNOs (Haumea, 2003 VS₂ and 2003 AZ₈₄) have been acquired with Herschel/PACS in search of thermal light curves, with the following main results:

• The thermal light curve of Haumea is clearly detected at 100 and 160 μm, superseding the early results of Lellouch et al. (2010) with Herschel/PACS. Both light curves are correlated with the optical one, implying primarily shape-driven light curves. Nonetheless, the 100-μm data indicates a small extra flux at rotational phases affected by the optical dark spot.

• Thermophysical modeling of the Haumea thermal light curves indicates an overall surface with an extremely small thermal inertia (Γ < 0.5 MKS and probably Γ < 0.2 MKS) and high phase integral (q ~ 0.8 for Γ = 0.1 MKS and no surface roughness), which will be even higher if surface roughness is present.

• The energetic and thermophysical properties of Haumea’s dark spot appear to be only modestly different from the rest of the surfaces, with changes of only ~+0.05–0.1 MKS in thermal inertia or ~+0.1 in phase integral. We favor the case for a small increase of thermal inertia in the dark region.

• The thermal light curve of 2003 VS₂ is not firmly detected at 70 μm and at 160 μm. However, Fourier fits to the thermal data are correlated with the optical light curve. The amplitude and mean flux of 2003 VS₂’s 70 μm light curve indicate a thermal inertia Γ = (2.0 ± 0.5) MKS.

• The thermal light curve of 2003 AZ₈₄ at 100 μm is not firmly detected. A thermophysical model applied to the mean thermal light curve fluxes and to all the Herschel/PACS and Spitzer/MIPS thermal data favors a close to pole-on (±30°) orientation.

• Our conclusion of extremely small thermal inertias for 2003 VS₂ and even smaller for Haumea statistically nicely matches inferences on the TNO/Centaurs population based on Spitzer/Herschel radiometry (Lellouch et al. 2013), including an albedo dependence of the thermal inertia. These authors interpreted their results in terms of highly porous surfaces, in which the heat transfer efficiency is affected by radiative conductivity within pores and increases with depth in the subsurface. For heat conduction dominated by radiation, the thermal inertia is essentially proportional to ${\mathit{r}}_{\mathrm{h}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{4}}$ (Delbo et al. 2015), or to ${\mathit{r}}_{\mathrm{h}}^{\mathrm{-}\mathrm{(}\mathrm{0.9}\mathrm{-}\mathrm{1.0}\mathrm{)}}$ if the temperature dependence of the specific heat of ice is taken into account (Lellouch et al. 2013). Our thermal inertia for the three objects (2.0 ± 0.5 MKS for 2003 VS₂ at r_h = 36.5 AU, ~0.2 MKS for Haumea at r_h = 51 AU, and 0.5−3.0 MKS for 2003 AZ₈₄ at r_h = 45 AU) convert into Γ = 10 − 135 MKS and Γ = 4 − 50 MKS at 1 AU for the two temperature-dependence cases, respectively. While somewhat even lower, these numbers compare generally well with the thermal inertias of large (>100 km) asteroids (10−300 MKS, corrected to 1 AU, see Fig. 9 in Delbo et al. 2015), where the smallest values are indicative of fine grain regolith. Recently, Ferrari & Lucas (2016) re-addressed the general issue of low thermal inertias in outer solar system bodies (including icy satellites), and pointed out several other important factors, in addition to surface porosity, affecting surface effective thermal inertias. One such factor is the quality of grain contact (i.e. tight or loose) in determining solid-state conductivity. For H₂O-ice covered surfaces, another factor, already recognized by Lellouch et al. (2013), is the amorphous vs. crystalline state of water, as the two states are associated with different bulk conductivities (and different temperature dependence thereof). On the basis of a detailed physical model of conductivity, including radiative conductivity, Ferrari & Lucas (2016) were able to reproduce the order of magnitude and heliocentric dependence of the thermal inertias measured by Lellouch et al. (2013), by invoking loose contacts in a moderately porous regolith of sub-cm-sized grains made of amorphous ice. Since water ice, when detected on the surface of TNOs, is usually in cristalline form, this scenario implies the presence of amorphous ice at cm depths below a thin layer of crystalline ice.

Acknowledgments

Herschel is an ESA space observatory with science instruments provided by European led Principal Investigator consortia and with important participation from NASA. Herschel data presented in this work were processed using HIPE, a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia. This research is partially based on observations collected at Centro Astronómico Hispano Alemán (CAHA) at Calar Alto, operated jointly by the Max-Planck Institut fur Astronomie and the Instituto de Astrofísica de Andalucía (CSIC). This research is partially based as well on observations carried out at the Observatorio de Sierra Nevada (OSN) operated by Instituto de Astrofísica de Andalucía (CSIC). The research leading to these results has received funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement No. 687378. P. Santos-Sanz and J. L. Ortiz would like to acknowledge financial support by the Spanish grant AYA-2014-56637-C2-1-P and the Proyecto de Excelencia de la Junta de Andalucía J.A. 2012-FQM1776. C. Kiss acknowledges financial support from NKFIH grant GINOP-2.3.2-15-2016-00003. E. Vilenius was supported by the German DLR project number 50 OR 1108. R. Duffard acknowledges financial support from the MINECO for his Ramon y Cajal Contract.

References

Appendix A: Tables

In this appendix, we include the light curve photometric results obtained from the Herschel/PACS observations of Haumea, 2003 VS₂ and 2003 AZ₈₄.

All Tables

All Figures

	Fig. 1 Images (top line: 100 μm; bottom line: 160 μm) extracted from the 20 Jun. 2010 light curve of Haumea. Left: original images. Right: background-subtracted images. The Field of View (FOV) is 2.5′× 2.5′. Haumea is at centre.
In the text

Fig. 2 Haumea’s thermal light curve at 100 μm (top) and 160 μm (bottom), combining data from the two visits (see text for details). The black dashed curve is the “scaled optical light curve” obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The part of this curve that is affected by Haumea dark spot, according to the Lacerda et al. (2008) preferred model, is outlined in yellow. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 455 188.720000 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.1 for further details). The uncertainty in the rotational phase is ±0.001.

In the text

	Fig. 3 Comparison of the Haumea 100 μm light curve with homogeneous models (no spot) having various values of the thermal inertia Γ in the range 0.0 to 0.5 MKS. No surface roughness effects are included (i.e. η = 1). For each Γ value, the phase integral is adjusted to match the mean flux level and light curve amplitude. The lowest Γ values (0.0−0.2 MKS) provide the best fit to the data.
In the text

	Fig. 4 Models with variable phase integral in the spot region. In these models, the spot thermal inertia is Γspot = 0, as in the no-spot region. The red curve (having q = 1.00 in the spot) represents the first best fit model for the 100 μm data.
In the text

Fig. 5 Models with variable thermal inertia in the spot region. In these models, the spot phase integral is kept at the same value as in the no-spot region. The green, black, and blue curves are for Γ = 0.1 MKS and q = 0.778 outside of the spot, and Γ = 0.0, 0.1, and 0.2 MKS in the spot. The red curve represents the second best fit model at 100 μm with Γ = 0.0 MKS and q = 0.851 outside of the spot, and Γspot = 0.05 MKS in the spot.

In the text

	Fig. 6 The two best fit models from the 100 μm light curve (see Figs. 4 and 5) applied to the 160 μm light curve (red and blue curves, respectively). The green curve is the best-fit model in terms of least squares adjusted to match the 160 μm phase, having Γ = 0.5 MKS and q = 0.68 throughout Haumea’s surface.
In the text

	Fig. 7 Images (top line: 70 μm; bottom line: 160 μm) extracted from the thermal time series curve of 2003 VS2. Left: original images. Right: background-subtracted images. The FOV is 2.5′× 2.5′. 2003 VS2 is at centre.
In the text

Fig. 8 Thermal time series curve at 70 μm (top) and 160 μm (bottom) for 2003 VS2. The black dashed curve is the scaled optical light curve obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 452 992.768380 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.2 for further details). The uncertainty in the rotational phase is ±0.01.

In the text

	Fig. 9 Thermophysical model fits of the 2003 VS2 70 μm data for various values of the thermal inertia. From top to bottom: Γ = 1.0, 2.0 and 3.0 MKS. Shape model and other physical input parameters are given in the text.
In the text

Fig. 10 Left panel: one of the 2003 AZ84 images at 100 μm (green) with the original background. Right panel is the same map after the application of the double-differential background subtraction technique, as described in Kiss et al. (2014). 2003 AZ84 is the source at the center of the images. In the right panel, the black source at the bottom left of 2003 AZ84 is the negative image of 2003 AZ84 resulting from the application of this background-removing technique. The FOV of each image is 2.5′× 2.5′.

In the text

Fig. 11 Thermal time series curve at 100 μm for 2003 AZ84. The black dashed curve is the scaled optical light curve obtained by rescaling the optical brightness to match the mean of the thermal fluxes. The grey solid curve is a Fourier fit to the thermal data. The reference for zero phase is at JD = 2 453 026.546400 days, uncorrected for light-time. Rotational phases have been computed using light-time corrected Julian dates and light-time corrected zero date (see caption of Table A.3 for further details). The uncertainty in the rotational phase is ±0.01.

In the text

	Fig. 12 Absolute PACS and MIPS fluxes for 2003 AZ84 with various thermophysical models overplotted: pole-on solution, pole-on +30°, pole-on +60°, equator-on. All the models use the occultation size and reasonable assumptions for the thermal properties.
In the text