© ESO, 2012

1. Introduction

Since the detection of the first trans-Neptunian object (TNO) 20 years ago – 1992 QB₁ – (Jewitt & Luu 1992) besides Pluto, about 1550 TNOs/Centaurs have been discovered (as of October 2011). Different dynamical classes, according to their orbital parameters, have been identified within this distant icy population: classical objects with moderate eccentricity orbits; resonant objects trapped in mean motion resonances with Neptune (Plutinos, in the 3:2 resonance with Neptune, represent the most populated resonant group, Jewitt & Luu 2000); scattered disk objects (SDOs) with high eccentricities and sometimes high inclinations, resulting from close encounters with Neptune; detached objects, whose perihelia are well outside Neptune’s orbit (Gladman et al. 2002), and which likely constitute a large population in the outer Kuiper Belt; centaurs, which reside between the orbits of Jupiter and Neptune, represent TNOs scattered towards the inner solar system, and so are usually considered as Kuiper Belt objects (KBOs). Although it is convenient to classify objects by “populations”, the limits between some of these populations (i.e. between SDOs, centaurs, high inclination/high eccentricity classicals, detached, etc.) are often “blurred” and may ultimately rely on definition (Gladman et al. 2008). The present paper deals with 15 SDO or detached objects, for which we first specify the adopted definition.

The terminology “scattered disk” was originally used for TNOs with orbits characterized by both a high eccentricity and a perihelion close to Neptune. It is believed that the majority of SDOs have unstable orbits due to close encounters with Neptune, but some of them could have survived since the origin of the solar system (Duncan & Levison 1997). Alternatively, the origin of the scattered disk could be related to other effects, including sweeping resonances that move the objects to this region (Malhotra 1995), passing star (or stars) (Morbidelli & Levison 2004, and references therein). Several approaches to the classification of objects as SDOs numerically integrate their orbits backward in time over dozens of millions of years (Elliot et al. 2005; Gladman et al. 2008). Here we adopt the Gladman et al. (2008) classification, which defines SDOs as those objects, independently of their origin, that are currently scattering actively off Neptune. SDOs are characterized by the rapid variation of semi-major axis in the orbit integrations. Specifically, Gladman et al. (2008) adopt a criterion similar to Morbidelli et al. (2004): the object is classified as an SDO if its semimajor-axis presents a variation  ≥1.5 AU during a 10-m.y. backward integration. The exact value used (1–2 AU variation) makes little difference, as SDOs undergo large semimajor-axis changes in short times. Then, SDOs in this nomenclature might be called “scattering” objects instead of “scattered” objects. Following this classification SDOs exist over a large semi-major axis range, not only restricted to a > 50 AU; rather the SDO population extends down to the Centaur region (beginning at a = 30 AU). At the other edge the inner Oort cloud would begin at very large semimajor-axes (a > 2000 AU, Gladman et al. 2008).

Detached objects (Delsanti & Jewitt 2006) are also known as “extended scattered disk objects” (ESDOs), though it is preferable to use the term “detached” because not all these objects have necessarily been placed in this region through scattering. Detached objects have orbits with perihelion distances unaffected (“decoupled”) from the giant planets, in particular Neptune. This population may include inner Oort cloud bodies or transitional objects between SDOs and inner Oort cloud objects. The limits in parameter space are unknown, and are expanding with the new discoveries of TNOs at larger and larger semimajor axes – e.g. a = 965 AU for 2006 SQ₃₇₂ – (Gladman et al. 2002). Elliot et al. (2005) adopted the Tisserand parameter (T) relative to Neptune to define the detached population. A problem with this definition is the large number of high inclination TNOs within the classical population. In the Gladman et al. (2008) scheme, numerical integrations separate SDOs from detached TNOs. A remaining issue is to decide where the detached population ends towards low eccentricities; Gladman et al. (2008) adopt a lower bound of e = 0.24. Possible processes that placed the detached objects on these orbits are currently being investigated (e.g. Gomes et al. 2008; Kenyon et al. 2008; Duncan et al. 2008, etc.).

Determining KBO sizes and albedos is not a straightforward task. Measurements from direct imaging are available for only a handful of them (Brown & Trujillo 2004; Brown et al. 2006). Most recently, sizes have become available from stellar occultations (Elliot et al. 2010; Sicardy et al. 2011). Obtaining size/albedo values for a larger sample is best achieved from thermal radiometry, whereby thermal and optical observations are combined within a thermal model. After the pioneering measurements on two objects from ISO (Thomas et al. 2000) and ground-based observations on a few others (see Stansberry et al. 2008, and references therein), Spitzer provided the first large such dataset, with about 60 Centaurs/KBOs measured at 24 and/or 70 μm (Cruikshank et al. 2005; Grundy et al. 2005, 2007; Stansberry et al. 2006, 2008; Brucker et al. 2009). This approach (multi-wavelength photometry) also constrains the thermal regime of the object, phenomenologically described by the “beaming factor” (Jones & Morrison 1974; Lebofsky & Spencer 1989; Harris 1998). The launch of Herschel Space Observatory in 2009 offered a new tool for this purpose, providing photometry over 6 bands from 70–500 μm. Early results on a handful of objects were presented in Müller et al. (2010), Lim et al. (2010), and Lellouch et al. (2010). Beyond their intrinsic scientific value, albedos are needed to properly interpret reflectance spectra. Furthermore, correlations between size, albedo, color, composition and orbital parameters may provide diagnostics on the dynamical, collisional and physical history of KBOs. The size distribution of TNOs is also a diagnostic relating to that history, but in the absence of measured sizes the distribution has typically been estimated from the luminosity function and an assumed albedo, an uncertain step given that the average albedos seem to be different for the different dynamical populations (Grundy et al. 2005; Brucker et al. 2009). Instead, and to constrain Kuiper belt formation/evolution models it is desirable to directly determine the size distribution of the different populations, and infer their individual total masses. Large bodies retain information about the angular momentum of the presolar nebula and their size distribution is dominated by accretion. On the other hand, small bodies – in the 50 to 100 km range – are the result of collisional evolution (Petit et al. 2008).

A large program (370 h of telescope time), named “TNOs are Cool: A survey of the trans-Neptunian region” (Müller et al. 2009), has been awarded on Herschel to measure the thermal flux of about 140 KBOs/Centaurs. We present here results on 15 SDO and detached objects. These observations lead to new sizes/albedos for 9 objects and improved results on 6 others. Results on 19 classical and 18 Plutinos are presented in companion papers by Vilenius et al. (2012) and Mommert et al. (2012), respectively.

2. Observations

The Herschel Space Observatory (Pilbratt et al. 2010) is a 3.5 m telescope launched on 14 May 2009, covering the far-infrared range by means of three focal plane instruments. The “TNOs are Cool” program consists of photometry of  ~140 objects, mostly with the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010).

The PACS instrument has 3 photometric bands centered at 70, 100 and 160 μm. These bands (termed “blue”, “green” and “red”, respectively) cover the 60−85 μm, 85−130 μm, 130−210 μm ranges. With two bolometer arrays, PACS can observe a field in two bands at a time: blue/red or green/red. The red detector is a 32    ×    16 pixel array, while the blue/green is 64    ×    32 pixels. The 3 bands image the same field of view of 3.50′    ×    1.75′. Projected onto the sky, pixel sizes are 6.4″    ×    6.4″ (red band) and 3.2″    ×    3.2″ (blue/green bands), respectively.

The observations discussed here were performed using the scan map mode. This mode was selected due to its better overall performance compared to the point-source (chop-nod) mode (Müller et al. 2010). Scan maps are performed by slewing the spacecraft at a constant speed along parallel lines, or “legs”. Our observations were typically performed at medium speed (20″/s). We use 10 scan legs in each map, separated by 4″. The length of each leg is 3.0′. Each such scan-map is repeated two to five times, depending on the expected brightness of the objects (see below). Within a single scan map, either the blue-red or the green-red channel combination is selected. For any band selection, two individual maps are made using two scanning directions with respect to the detector array: 70° and 110°, allowing to fit the scan pattern on the rectangular array in an optimum way.

The “TNOs are cool” target list was constructed by predicting the thermal fluxes in the PACS bands for the known TNOs using the “hybrid standard thermal model” (STM – see Sect. 4). The flux predictions were based on absolute visual magnitudes (H) from the Minor Planet Center (MPC), an assumed visual geometric albedo of 8% (based on previous studies using Spitzer), and geometric circumstances appropriate to the expected observing epoch (c. 2010). For objects with known albedos (for this study those were Eris: 70%, Typhon: 5%, 1996 TL₆₆: 4%, and 2002 PN₃₄: 5%) we used those values rather than the 8% assumption. Objects with the largest predicted thermal fluxes were selected as potential targets; that list was further refined by rejecting objects with predicted positional uncertainties  <5′′. The objects presented in this paper are a subset of the targets selected in the above manner. Once the targets were selected, integration times were chosen to provide signal-to-noise ratio (SNR) values of  ~10 in the blue and green channels. Integration times for individual observations were tailored by varying the repetition factor for the scan maps described earlier. Note that the values of H from the MPC were used for planning purposes only, and can differ significantly (by up to 0.5 mag) from the values ultimately used in our analysis.

It is important to note that the target selection process introduces a bias towards brighter (low H values) objects. A consequence of the assumed and generally fixed albedo is that the object sample favors large objects at large heliocentric distances (r_h). In addition, the object sample is affected by a discovery bias that: (i) favors high-albedo objects at high r_h, and (ii) introduces a correlation between the perihelion (q) and the heliocentric distance, because objects are most easily discovered near perihelion and haven’t moved much since discovery. These biases and correlations must be kept in mind when studying correlations between size, albedo and other physical and orbital parameters (see Sect. 5). As for the present SDOs/Detached sample, objects are listed in Table 1, including relevant orbital parameters, absolute magnitudes in V-band, rotational periods and amplitudes, taxonomic color class (when known), existence or not of satellites, and dynamical classification following Gladman et al. (2008). Observational circumstances are described in more details in Table 2.

Each target was visited twice, with identical observations repeated in both visits. During the first visit, we grouped 4 maps to be observed in sequence, corresponding to the two scan directions and the two channel combinations. The observation geometry is such that on the final combined map, the central  ~50″ × 50″ region has high coverage. Photometric accuracy is usually limited by instrumental noise in the blue channel and by background confusion noise in the red, while either of the two sources of noise may dominate in the green¹. The timing of the observations, i.e. the selection of the visibility window, was optimized to ensure the lowest far-infrared confusion noise circumstances (Kiss et al. 2005), in such a way that the estimated SNR due to confusion noise has its maximum in the green channel. The timing of the second visit was calculated such that the target moved 30−50′′ between the two visits. This ensures that the target position during the second visit is still within the high-coverage area of the map from the first visit. At the same time, this enables us to determine and subtract the background in the individual maps, so as to obtain a clean combined image where we extract the final flux (see Sect. 3.1).

3. Data reduction and photometry

In this section, we first describe the reduction process applied to the PACS data in order to obtain what we call “single” and “combined” maps, and then the photometric techniques used to determine the target fluxes and associated uncertainties in each band. Further details of the data reduction are given in Appendix A.

3.1. Data reduction

The data reduction of the PACS maps was performed within the Herschel interactive processing environment (HIPE², version 6.0.2044) by means of adapted standard HIPE scripts. Since our targets move slowly (a few ″/h), and the total observation time per scanning direction in one visit is only 10–25 min, we do not correct for the apparent motion.

After the application of the HIPE scripts we obtain one “single” map per visit, filter, and scan direction (i.e. in total 8 maps per object in the red, and 4 maps in the blue/green). The final re-sampled pixel scale of the single maps are 1.1″/pixel, 1.4″/pixel, and 2.1″/pixel for the blue (70 μm), green (100 μm) and red (160 μm) channels respectively.

We use these single maps to generate final “combined” maps on which the photometry will be performed. To generate the combined maps, we combine the two visits to the target to determine the background map and subtract it from each single map. Finally we co-add all the background-subtracted maps in the co-moving frame of the target. Similar methods were used to subtract background sources in Spitzer TNO observations, and are described in Stansberry et al. (2008) and Brucker et al. (2009). An alternate, simpler method, to obtain final maps is to co-add directly the original single maps (i.e. not background-corrected) in the target frame. This method is obviously less optimal than the previous one in terms of SNR, but provides a useful test to demonstrate that the background subtraction is not introducing any spurious effects. A more complete and detailed description about the data reduction process can be found in Appendix A, and will be published in Kiss et al. (in prep.-a).

3.2. Photometry

We use IRAF/Daophot³ to perform standard synthetic-aperture photometry on the final sky-subtracted, co-added maps for our targets. The target is typically obvious and located very near the center of the maps. For very faint targets and/or ones with some residual structure from the background subtraction (as can happen in the green and red channels), the ephemeris position of the target is used to position the photometric aperture.

Once the target is identified, we measure the flux at the photocenter position for aperture radii ranging from 1 to 15 pixels. We perform an aperture correction technique (Howell 1989) for each aperture radius using the encircled energy fraction for a point–source for the PACS instrument⁴. We construct an aperture-corrected curve-of-growth using these fluxes (red line of Fig. A.4). Given the optimized reduction technique above, the fluxes used to construct these aperture-corrected curves-of-growth need not be corrected for background contribution.

After inspection of the curves-of-growth, the optimum synthetic aperture to measure the target flux is selected, and is typically  ~1.0–1.25 times the point spread function (PSF) full width at half maximum (FWHM) in radius (PSF FWHM in radius is 5.2″/7.7″/12″ in blue/green/red bands respectively). Usually, the optimum aperture radius lies in the “plateau” zone of the curves-of-growth (see Fig. A.4). In cases in which there are other sources or artifacts very close to the target, we use a smaller aperture radius to avoid contamination of the target flux.

Uncertainties on the flux measurements are estimated by means of a Monte-Carlo technique, in which 200 artificial sources are implanted on the actual background-subtracted final maps. These artificial sources have the structure of the PSF for the specific (blue, green, or red) channel. Sources are randomly implanted on a square region of 50″    ×    50″ around the target photocenter, with an exclusion zone corresponding to a circular region around the target photocenter with a radius = 2    ×  PSF FWHM. We measure and aperture-correct the fluxes of these artificial sources using an aperture radius of 5 re-sampled pixels, which is the median optimum aperture radius. Uncertainties are computed as the standard deviation of these 200 fluxes, and finally multiplied by a factor . The latter factor is due to the fact that given the way the final maps are generated, the remaining background is measured only once in the immediate vicinity of the real target but two times in the rest of the image (i.e. for the artificial sources). Note that except for the background subtraction, all the above reduction and photometric steps were validated on standard stars with known flux.

The fluxes and associated uncertainties are finally color-corrected for each band. Flux densities within the PACS photometric system are defined as those that a source with a flat spectrum would have at the PACS reference wavelengths of 70, 100 and 160 μm (Poglitsch et al. 2010). Converting the in-band flux to monochromatic flux densities therefore depends on the spectral energy distribution of the source and the transmission profile for each PACS band. Color-corrections for each filter are tabulated as a function of color temperature⁵ but are usually within a few percent except at the very low end (T < 30 K) of the temperature range. For the color temperature, we used an approximation of the disk-averaged temperature, i.e. T_color = T_ss/2^1/4, where T_ss, the subsolar temperature (see Eq. (2) ), was obtained from a first run of the thermal model described below, using at this step the in-band (i.e. color-uncorrected) fluxes. T_color values for our sample range from 25 K for Eris to 77 K for 2002 PN₃₄. The mean color correction factors for all our sample are: 0.991  ±  0.028 (blue band), 0.989  ±  0.006 (green band), and 1.011  ±  0.017 (red band). The impact of these factors is smaller than the uncertainties of our measurements.

The absolute flux calibration of PACS is based on standard stars and large main belt asteroids and has uncertainties of 3% for the blue/green channel, and 5% for the red channel⁶. The final flux uncertainties we use in our modeling are the root-sum-square combination of these absolute calibration uncertainties and the uncertainties we obtain from the photometric measurements for individual targets. A summary of the final, color-corrected and absolute-calibrated, fluxes of our SDOs/Detached sample is given in Table 3.

3.3. Spitzer-MIPS observations

About 102 TNOs/Centaurs have been observed with the 0.85 m Spitzer Space Telescope (Werner et al. 2004) using the Multiband Imaging Photometer (MIPS) instrument (Rieke et al. 2004). The MIPS 24 μm channel consist of a 128 × 128 pixel detector and the 70 μm channel of a 32 × 32 pixel detector. For these 2 channels the telescope-limited resolution is 6″ and 18″ respectively. Fluxes for 60 objects have been published in Stansberry et al. (2008) and Brucker et al. (2009) but others are unpublished. Eight targets of our sample have available (published or unpublished) Spitzer-MIPS flux density measurements at 24 and/or at 70 μm (see Table 4). Unpublished data have been processed in a similar way as in Stansberry et al. (2008) and Brucker et al. (2009), using the calibration (Gordon et al. 2007; Engelbracht et al. 2007) and color-corrections (Stansberry et al. 2007) described in the references. As in Stansberry et al. (2008), we adopt absolute calibration uncertainties of 3% and 6% for the 24 and 70 μm observations (50% larger than the uncertainties derived for observations of stellar calibrators). The larger uncertainties allow for effects from the background-subtraction technique, the faintness of the targets, and uncertainties in the color corrections. The background subtraction techniques used for our Herschel data (see Sect. 3.1) are in fact derived from techniques originally developed for the MIPS data reductions. The previously unpublished fluxes presented here are based on new reductions of the MIPS data, utilizing updated ephemeris positions for the targets. This allows for more precise masking of the target when generating the image of the sky background, and for more precise placement of the photometric aperture. In most cases the new fluxes are very similar to the previously published values for any given target, but in a few cases significant improvements in the measured flux and SNR were achieved.

As discussed below, the combination of Spitzer-MIPS and Herschel-PACS fluxes significantly improves results from the thermal modeling. The Spitzer 70 μm band is similar to the Herschel 70 μm channel, providing a consistency check on data from the two observatories. The inclusion of the Spitzer 24 μm flux, when available, strongly improves the ability to constrain the surface temperature distribution, i.e. in general, allows us to determine the beaming factor (see Sect. 4). Published, re-analyzed, and unpublished Spitzer-MIPS flux density measurements available for our sample are given in Table 4. Note that the effective wavelengths for these fluxes are 23.68 and 71.42 μm.

4. Thermal model

4.1. General formulation

The combination of thermal-infrared and optical data makes it possible to infer an objects’ diameter (D), geometric albedo (p_V) and constrain its thermal properties. This requires the use of a thermal model to calculate the disk-integrated thermal emission. Here we use the so-called hybrid standard thermal model (hybrid STM, Müller et al. 2010; Stansberry et al. 2008). Except for the fact that the hybrid STM assumes zero phase angle (which is an excellent approximation for TNOs), it is identical to the Near-Earth Asteroid Thermal Model (NEATM, Harris 1998). NEATM was originally developed for Near Earth Objects (NEOs) but is applicable to all atmosphereless bodies and is being used in companion papers by Mommert et al. (2012) and Vilenius et al. (2012). In this model, the temperature distribution follows instantaneous equilibrium of a smooth spherical surface with solar input, but is modified by a factor η (beaming factor), which represents empirically the combined effects of thermal inertia and surface roughness. As shown by Harris (1998) and Stansberry et al. (2008), multi-wavelength measurements can be used to constrain η if the data are of sufficient quality and encompass both the short-wavelength and the peak and/or long-wavelength portions of the spectral energy distribution (SED) of the target. In this case the data give a direct measurement of the color-temperature of the emission, which is sensitive to the value of η. Our preferred approach was to simultaneously model Herschel and Spitzer measurements (of course accounting for their different observer-centric, Δ, and heliocentric, r_h, distances) and fit an η factor to the combined dataset. We call this the “floating-η” approach. When Spitzer data were unavailable (or when this method led to implausible η values; see Mommert et al. 2012) we instead used a “fixed-η” approach, with η = 1.20 ± 0.35, the mean value derived by Stansberry et al. (2008) from Spitzer observations of TNOs.

The model therefore fits the Herschel (3 data points) or Herschel + Spitzer (5 data points) fluxes in terms of two (D and p_V) or three (D, p_V, and η) parameters. The model also requires the choice on a phase integral q, that we assumed to be related to p_V through q = 0.336p_V + 0.479 (Brucker et al. 2009). Diameter and geometric albedo are further related through: (1)where \begin{lxirformule}$H$\end{lxirformule} is the absolute magnitude (nominally in V band), and V_⊙ = −26.76 ± 0.02 (Bessell et al. 1998) the V magnitude of the Sun. A grey emissivity (ϵ = 0.9) was assumed when calculating the local temperatures and monochromatic fluxes. With these definitions the subsolar temperature is given by: (2)where \begin{lxirformule}$F_{\odot}$\end{lxirformule} is the solar flux at 1 AU. The solution for the D, p_V (and η) parameters are determined from a χ² minimization of the deviation of the model from the observed fluxes, accounting for their individual 1-σ error bars. Note that when only upper limits were available, they were treated as non-detections (i.e. zero flux) with a 1-σ uncertainty as the upper limit. In the case of multiple systems D is the area-equivalent diameter, i.e., that of the sphere with the same total surface area, because none of the binaries are resolved in our data.

4.2. H magnitudes

Equation (1) above outlines the role of absolute magnitude (H) in the thermal modeling. Determining accurate values of H, preferably in V band, is therefore an essential task for this project. We describe hereafter our adopted approaches.

In the best situation, photometrically calibrated observations are available at a variety of phase angles. In this case, correcting the observed V magnitudes for the heliocentric (r_h) and geocentric (Δ) distances provides m_V(1,1) = V−5·log (r·Δ). A linear fit of the m_V(1,1) magnitudes vs. phase angle α (m(1,1) = H_V + β·α), taking into consideration measurement errors, then provides the phase coefficient β (mag/deg) and the absolute magnitude H_V, and their associated error bars. Whenever data are unsufficient in V band but adequate in R band, we similarly and alternatively infer H_R, which by means of the thermal model will yield the red geometric albedo (p_R).

In the case when not enough data in either V or R are available, or if the above approach fails due to poor quality data, we adopt a mean β = 0.10 ± 0.04 mag/deg which is the average value (excluding Pluto) from Belskaya et al. (2008). We then derive H magnitudes for each individual photometric measurement and use the weighted average result.

When photometrically-calibrated magnitudes are not available in the literature, we use astrometric observations compiled by MPC, as those include rough estimates of the optical magnitudes (most often in R band). The above approaches are then applied to this data, bearing in mind that individual V or R magnitudes reported in MPC data may have large (and systematically biased) uncertainties up to 0.5 mag (Romanishin & Tegler 2005).

When full lightcurve information is available (i.e. the amplitude is known and the lightcurve can be temporally phased to the Herschel measurements), the H magnitude is corrected accordingly. When the lightcurve amplitude (or an upper limit) is known, but the period is unknown or phasing is impossible, 88% of half peak-to-peak amplitude (i.e. 1σ = 68% of the values of a sinusoid are within this range) is quadratically added to the uncertainty on H. Finally, if there is no lightcurve information, we assume peak-to-peak amplitude less than 0.2 mag, as based on Duffard et al. (2009),  ~70% of the objects have a peak-to-peak amplitude less than this. In this case we quadratically add 0.09 mag (i.e. 88% of 0.2/2) to the H uncertainty. Absolute magnitudes quoted in Table 1 have been obtained as just described.

Note finally that in case we have to use R magnitudes (i.e. for 2007 OR₁₀, 2007 RW₁₀ and 1999 KR₁₆) modeling provides the R geometric albedo (p_R) instead of V albedo (p_V). p_R can be converted to p_V through: (3)where (V − R)_⊙ = 0.36 mag is the Sun’s V − R, and (V − R)_obj that of the object. We assume (V − R)_obj = 0.55    ±    0.11, the mean value from Hainaut & Delsanti (2002) for 66 SDOs/Detached objects, except for 2007 OR₁₀ (see discussion for this object in Sect. 5.3).

4.3. Uncertainties

Within the framework of the thermal model, uncertainties on the fitted (D, p_V, η) parameters were obtained using a Monte-Carlo approach (see Mueller et al. 2011) in which 1000 synthetic datasets were randomly generated using the uncertainties in the measured fluxes. Note that the uncertainty in the H magnitude (and the  ± 0.35 uncertainty in η, in the case of the “fixed-η” approach) was also accounted for in the Monte-Carlo approach. Naturally, we eliminated those Monte-Carlo simulations which generated negative fluxes (or negative η). The diameter, albedo (and η) of each of the 1000 synthetic objects were determined, and their distributions were used to determine the error bars on these parameters. As outlined in Mueller et al. (2011), these distributions are generally not Gaussian (especially the albedo distribution). Therefore, we adopted the median of the Monte-Carlo results as the nominal value, and asymmetric error bars to include the 68.2% of the results.

Figure 1 shows best fits for all objects in our sample using only the Herschel-PACS measurements. All these fits are “fixed-η”. For the 8 objects that also have Spitzer fluxes, Fig. 2 shows best fits (“floating-η”, except for Eris, which is not detected at 24 μm) for the combined Herschel and Spitzer measurements. In that Figure, the same best fit model is plotted twice, once (solid line) for the Herschel (r_h, Δ) conditions and once for the Spitzer (dashed line) conditions. In general, this shows that the model appropriately accounts for rather different Herschel vs. Spitzer fluxes associated to different r_h and Δ (e.g. Ceto, 2002 PN₃₄). Many of the fits shown in Figs. 1 and 2 are satisfactory, in that the model fits all data points within or close to the 1-σ error bars. However, in a number of cases, the fits are obviously bad (e.g. 1999 KR₁₆, 2005 QU₁₈₂, Typhon); in those cases, the data points are mutually inconsistent and the model is no more than a mere compromise between them. Given that the PACS green vs. blue observations (and the Spitzer vs. Herschel data) are not taken simultaneously, such behaviour could in principle result from rotational variability of the objects. An alternate physical explanation would be that the fluxes are affected by spectral features, e.g. depressing the 100-μm flux with respect to the 70-μm and 160-μm values, or more generally that our simple STM-like approach does not capture the complexity of the thermal emission. However, the discrepancies are often large and these explanations may not be quantitatively plausible. In those cases, we conclude that the error bars we have determined on the fluxes are probably underestimates of the true error bars, though we are unable to track down the reasons for that.

In the situation where the formal measurement errors are too small (or similarly that the model is not completely adequate), the above described Monte-Carlo method will underestimate the uncertainties on the solution parameters. To handle this problem, we adopted a “rescaled error bar” approach described in detail in Appendix B (Sect. B.1). Essentially, it consists of “inflating” the measurement error bars by a suitable factor before running the Monte-Carlo simulation. This approach was applied to the six objects for which the radiometric fit is poor (see Figs. 1 and 2): 1999 KR₁₆, 2002 PN₃₄, 2002 XU₉₃, 2005 QU₁₈₂, 2007 UK₁₂₆, and Typhon, with rescaling factors varying from 1.45 to 1.99. Fluxes for the remaining nine objects are satisfactorily fit within the 1-σ flux uncertainties and did not require any adjustment of the measurement errors. We finally stress that the diameter, albedo, and beaming factor values and their uncertainties reported in the next section are obtained in the framework of the hybrid STM, which may have its own limitations (see Harris 2006, for a discussion in the case of NEOs). More physical models (e.g. TPM) may lead to somewhat different values and error bars, as illustrated by Müller et al. (2010) in the case of a few objects observed in the Science Demonstration Phase (SDP) of the “TNOs are cool” programme.

5. Results and discussion

5.1. General results

Results of the model fits in terms of the objects diameter, albedo and beaming factor are gathered in Table 5. When Spitzer-MIPS measurements are available (8 objects), the model was run twice, once considering the PACS fluxes only, and once combining PACS and MIPS. As illustrated in Table 5, these two fitting options give consistent solutions within error bars. The PACS+MIPS solution is always adopted as the preferred one (outlined in bold face in Table 5), as it (i) provides reduced error bars (ii) permits the determination of the beaming factor.

Table 5 includes first size/albedo determination for 9 objects (4 SDO and 5 detached objects). Six of them are based on PACS-only measurements (2007 OR₁₀, 2007 RW₁₀, 1999 KR₁₆, 2005 QU₁₈₂, 2007 OC₁₀, and 2007 UK₁₂₆), and the remaining three make use of PACS plus unpublished-MIPS data (2001 FP₁₈₅, 2002 XU₉₃, and 2003 FY₁₂₈). Published sizes and albedos exist for the other 6 objects of our sample (1996 TL₆₆, 2002 PN₃₄, Ceto, Typhon, 2005 TB₁₉₀ and Eris) but we here present fluxes for these objects at wavelengths not observed before (i.e. 100 μm and 160 μm), and the combination with earlier measurements leads to improved estimations of sizes, albedos, and beaming factors (see Table 5).

The V geometric albedos for our targets vary from 3.8% (2002 XU₉₃) to 84.5% (Eris), with an unweighted mean value of 11.2  ±  7.6% excluding Eris. If weighted by the relative errors (i.e. 1/(σ/p_V)²), the mean albedo is 6.4  ±  6.3% (for σ we used the mean of upward and downward uncertainties on p_V). This error bar on the mean albedo includes the dispersion of the albedo measurements and their individual RMS errors. It is dominated by the dispersion between individual measurements, so it has a limited significance. This variation of geometric albedos can be seen in Fig. 3-left. For SDOs, there is a lack of albedos  >20% and a clustering in the [0–5%] range. In contrast, detached objects in our sample all have albedos  >5%, with a peak in the [10–20%] range. The unweighted/weighted mean is 6.9/5.2% for the SDOs and 17.0/12.3% for detached objects excluding Eris.

Diameters vary from 112 km (2002 PN₃₄) to 2454 km (Eris), with a remarkably uniform distribution (Fig. 3, right). However, the distribution of diameters is seemingly different between SDOs and detached objects. SDO diameters are smaller than 400 km, except for an outlier at 1280 km (2007 OR₁₀), while all detached objects are larger than 200 km. This apparent tendency of detached objects to be larger than SDOs is, however, probably an effect of the discovery bias (see discussion in Sect. 2), as detached objects are typically discovered at larger heliocentric distances than SDOs.

Except for Eris, where the “floating η” fit fails, the combination of Spitzer and Herschel data allows us to determine the beaming factor for 7 objects (6 SDO and 1 detached – see Table 5 –). The (weighted by absolute errors) mean-η obtained from these results is 1.14  ±  0.15. This is fully consistent with the mean 1.20  ±  0.35 η-value derived from the Spitzer dataset (Stansberry et al. 2008) and justifies the use of this mean value for those cases where η cannot be determined. Beaming factors only barely larger than  ~1 (as opposed to η ≫ 1) indicate that KBOs in general and our sample objects in particular are close to the slow rotator regime (i.e. STM), which is a remarkable result given the very low temperatures prevailing on these objects and their usually short rotation periods. Except in the unlikely situation where they are all seen pole-on, this implies very low thermal inertia and/or very large surface roughness. Lellouch et al. (2011) determine the thermal inertia of Charon’s surface to be Γ = 10−20 J m^-2 s^−1/2 K^-1 (MKS) and find this to be equivalent to η = 1.21−1.40. Yet Charon is a slow rotator (6.4 day) and a rough scaling indicates that the same η range, applied e.g. to Haumea or Orcus, would be representative of even much smaller thermal inertias (Γ = 0.5−3 MKS). Hence a general conclusion is that unless they are very rough, KBO surfaces must have very low thermal conductivity, suggesting a high degree of surface porosity. Similar thermal inertias and η’s for KBOs were obtained in Müller et al. (2010).

5.2. Comparison with earlier results

Six objects of our sample (1996 TL₆₆, 2002 PN₃₄, Ceto, Typhon, 2005 TB₁₉₀ and Eris) have previously published sizes/albedos. We here briefly compare results for the first five, and discuss Eris in more detail in next subsection. Note that 2007 UK₁₂₆ has been reported as binary system (Noll et al. 2009), but there is no mass estimation for this object, for this reason so we do not discuss it here.

(15874) 1996 TL₆₆. Our preferred solution for this SDO (D = 339    ±    20 km, %, ) is at odds with the Stansberry et al. 2008 results ( km, %, η = 1.8 ± 0.3) based on Spitzer-MIPS fluxes at 24 and 70 μm. Our reprocessing of these data leads to very different fluxes at 24 μm (0.32 mJy vs. 0.38 mJy in Stansberry et al. 2008) and especially at 70 μm (10.1 mJy vs. 22 mJy). Fitting these revised Spitzer fluxes alone gives a best fit (D = 328 km, p_V = 11.8%, η = 1.12) in excellent agreement with our preferred solution.

(73480) 2002 PN₃₄. Although the 70 μm fluxes from Spitzer and Herschel are very different (Fig. 2), this can be ascribed to the different observing distances. Our preferred results are: D = 112 ± 7 km, p_V = 4.9 ± 0.6%, , in a good agreement with Stansberry et al. (2008) results: (D = 120 ± 10 km, %, ). This object is the smallest one in our SDOs/Detached sample.

(65489) Ceto/Phorcys is a binary SDO, previously observed with Spitzer-MIPS (Grundy et al. 2007; Stansberry et al. 2008). Model fits indicate that the optimum model (D = 281 ± 11 km, p_V = 5.6 ± 0.6%, η = 1.04 ± 0.05) matches the 24 μm flux but somewhat overestimates the Spitzer 70 μm flux, indicative of a colder object than the Spitzer data would suggest. It is therefore not surprising that the Spitzer-only solution of Stansberry et al. (2008) has a lower η (), lower D ( km) and higher p_V (%).

Grundy et al. (2007) estimate a mass of (5.41 ± 0.42) × 10¹⁸ kg for the Ceto/Phorcys system. Using an equivalent diameter of  km, they obtain a bulk density of  g cm^-3. We derive a smaller bulk density of  g cm^-3 due to our larger equivalent diameter (D = 281 km which gives D_Ceto = 223    ±    10 km, and D_Phorcys = 171    ±    10 km, assuming the same albedo for the two components and Δm = 0.58    ±    0.03 mag). This low bulk density is compatible with a porous surface, which is also supported by the very low η = 1.04 value.

(42355) Typhon/Echidna is a binary SDO. Our preferred results are D = 185 ± 7 km, p_V = 4.4 ± 0.3%, η = 1.48 ± 0.07, consistent with the Stansberry et al. (2008) solution obtained by applying a fixed η = 1.2    ±    0.35 model to the 24 μm flux ( km, ). In contrast, these results are inconsistent with the preliminary Herschel results (D = 138 ± 9 km and p_V = 8.0 ± 1.0%), based on SDP data (Müller et al. 2010). The latter observations were performed in chop-nod mode (less sensitive than scan-map for point sources) and led to lower 70 and 100 μm fluxes and no detection at 160 μm.

Grundy et al. (2008) obtain a mass of (9.49    ±    0.52) × 10¹⁷ kg for this binary system. Using an equivalent diameter of  km, they obtain a bulk density of  g cm^-3. We derive a nominally even smaller (but compatible within error bars) bulk density of  g cm^-3 from a slightly larger equivalent diameter (D = 185 km, i.e. D_Typhon = 162    ±    7 km, and D_Echidna = 89    ±    6 km, assuming same albedo and Δm = 1.30    ±    0.06 mag). Again, this extremely low bulk density points to a porous body.

(145480) 2005 TB₁₉₀. This detached object was not observed by Spitzer-MIPS, but was already measured in chop-nod mode by Herschel. Our new results, using a fixed-η mode (D = 464    ±    62 km, ) are reasonably consistent (to within 1-σ) with the preliminary values published by Müller et al. (2010) (D = 375    ±    45 km, p_V = 19.0    ±    5.0%).

5.3. Two large objects, 2007 OR₁₀ and Eris

2007 OR₁₀. Its H_R magnitude of  ~2 indicates that this very distant object (currently at 86.3 AU) must be one of the largest TNOs, unless it exhibits a very high albedo. Brown et al. (2011) recently reported that its spectrum is markedly red in the near-infrared and shows features due to water ice. This is an unusual combination as most KBOs with water ice absorptions have neutral colors, but this characteristic is shared by Quaoar. Although their data are unsufficient to show a detection of methane ice, Brown et al. (2011) used the analogy with Quaoar, where methane ice is present (Schaller & Brown 2007a), to suggest that the red color of 2007 OR₁₀ is due to irradiation of small amounts of methane ice. They further assumed an albedo equal to that of Quaoar (nominally 0.18 for an adopted Quaoar radius of 960 km, but allowing for 50% uncertainties), and showed, based on the Schaller & Brown (2007b) sublimation model, that 2007 OR₁₀ is in a similar volatile retention regime as Quaoar, with the main volatile species (N₂, CO, CH₄) being marginally stable over the solar system age.

2007 OR₁₀ is clearly detected in the three PACS channels. In the lack of V measurements, we assumed H_R = 1.96    ±    0.16 mag based on the fit of eight good quality R-measurements from MPC data (see Table 1). A fixed-η fit is satisfactory (Fig. 1), providing D = 1280 ± 210 km, . This puts the object in the Charon, Makemake, Haumea, etc. size class, i.e. a plausible candidate for dwarf planet status. To convert p_R to p_V we use Eq. (3) from Sect. 4.2, but use Quaoar’s V − R value = 0.67  ±  0.02 (Fornasier et al. 2004). This yields , nicely consistent with the value assumed by Brown et al. (2011). 2007 OR₁₀ thus indeed appears as a sibling to Quaoar, probably of even larger size.

(136199) Eris. Although Eris was observed by Spitzer, it was not detected at 24 μm (with a 0.005 mJy upper limit), preventing us from determining its beaming factor. Assuming therefore η = 1.20    ±    0.35, we obtain a best fit solution for D = 2454 ± 117 km, p_V = 84.5 ± 8.8%. Using the mass determined by Brown et al. (2007) – recomputed in Grundy et al. (2011) – for the Eris/Dysnomia system ((1.688  ±  0.035)  ×  10²² kg) we obtain a bulk density of  g cm^-3, assuming the same albedo for Eris and Dysnomia (D_Eris = 2434    ±    117 km, and D_Dysnomia = 316    ±    23 km). It is probably more realistic to assume that Dysnomia’s surface is much darker than Eris’. Assuming for definiteness a 5 times lower albedo, we obtain D_Eris = 2356    ±    117 km and D_Dysnomia = 685    ±    50 km, yielding a global density of  g cm^-3, similar to the Brown et al. (2007) value (2.3   ±   0.3 g cm^-3).

Near-IR spectroscopy (Tegler et al. 2010; Merlin et al. 2009; Brown et al. 2005) indicates that Eris’ surface composition is very similar to Pluto’s, with the presence of methane ice along with another dominant ice, presumably nitrogen, suggestive of a high albedo. Prior to our measurements, Eris’ size and albedo have been determined by several groups and methods. The first measurement was achieved by Bertoldi et al. (2006) who used IRAM thermal observations at 1200 μm and obtained  km and p_V = 60 ± 10%. Subsequently, Brown et al. (2006) estimated D = 2400 ± 100 km, p_V = 86 ± 7% from direct imaging with the Hubble Space Telescope. Based on the Spitzer 70 μm flux, Stansberry et al. 2008 obtained  km, from Spitzer-MIPS and hybrid-STM fit. Within uncertainties, our results agree with all these values except with those of Bertoldi et al. (2006), which overestimate the size and underestimate the albedo. Most recently, a stellar occultation by Eris on November 6th 2010 provided a very accurate determination of Eris’ diameter and albedo: D = 2326 ± 12 km, % (Sicardy et al. 2011), assuming a spherical shape. (Note however that allowing for elliptical shape, Sicardy et al. find a broader range of solutions, 68.3% (1-σ) of which have an effective diameter of 2330  ±  180 km.) Eris’ diameter is thus indistinguishable from Pluto’s given error bars on both, supporting the idea that Eris is a Pluto twin, with an even brighter surface, presumably due to a collapsed atmosphere at aphelion.

At face value and given error bars, our result on Eris’ diameter agrees nicely with that of Sicardy et al. (2011), especially if we assume that Dysnomia is intrinsically much darker than Eris. Similarly, our density for the system is in line with theirs (2.52  ±  0.05 g cm^-3). However, the consistency is actually not complete. Sicardy et al. (2011) reanalyzed the Spitzer 70 μm and IRAM 1200 μm fluxes in the light of their determined diameter and albedo. They considered both a standard thermal model (STM) and an isothermal latitude model (ILM), introducing an “η” (=η_R) parameter describing exclusively surface roughness effects (and ranging from η_R = 1.0 – no roughness – to η_R = 0.7 – large roughness –). In this framework, the only free parameter of their thermal model was the phase integral q. They found that the thermal measurements implied q = 0.49–0.66 for the STM case (fully consistent with the value of Saturn’s brightest satellites), but an implausible q = 0−0.24 for the ILM. On this basis, the STM was strongly favored, implying either a pole-on orientation or a negligible thermal inertia. The inconsistency with our analysis resides in the fact that we assumed η = 1.20    ±    0.35 and that our p_V solution implies q = 0.75, according to the assumed Brucker et al. (2009) expression. Another way of seeing this is to note that if we now assume the two limiting STM cases considered by Sicardy et al. 2011 (i.e. q = 0.49, η_R = 1.0, and q = 0.66, η_R = 0.7) and re-fit the PACS + MIPS data, we obtain best fit diameters of 2034 km and 2139 km respectively, and V-albedos larger than 100% in the two cases, inconsistent with the occultation results. Ultimately, the discrepancy can be traced to the fact that Sicardy et al. (2011) included the IRAM 1200 μm flux in their model, which as they noted tends to indicate warmer temperatures than the Spitzer 70 μm flux, while on the contrary our best fit model tends to underestimate the latter. The resolution of this problem will require a simultaneous fits of all datasets, considering in particular (i) non-spherical shapes (ii) separate thermal models for Eris and Dysnomia and (iii) more elaborate thermophysical models in which the thermal inertia, surface roughness and spin orientation are free parameters. This study is deferred to a subsequent paper (Kiss et al., in prep.-b). Note finally that the interpretation of thermal observations of a volatile-covered body such as Eris is expected to be even more complex than presented here. An additional complication is that volatile transport will result in surface variegation, i.e. multiple surface albedos. This is demonstrated in the case of Pluto, which does show large albedo spatial variations and whose thermal data require three-terrain model (Lellouch et al. 2011; Mommert et al. 2012). Another example is Makemake, the visual appearance of which is unknown but for which a two-terrain model is required to fit the radiometric data (Lim et al. 2010).

5.4. Search for correlations

5.4.1. Description of the method

Using Spearman-rank correlation coefficient (Spearman 1904), also known as Spearman-ρ, we have examined our results for correlations between albedo (p_v), beaming factor (η), diameter (D), absolute magnitude (H), color (B − R), perihelion distance (q), orbital inclination (i), semi-major axis (a), and heliocentric distance at the time of observation (r_h).

Spearman-ρ is distribution-free and less sensitive to outliers than most other methods but, like others, treats data-points as “exact” and does not take into account their possible error bars. Any variations in the measured data points within their error bars, may change the correlation coefficient, however. Furthermore, each correlation coefficient has its own confidence interval, which depends on the number of data points and on the magnitude of the correlation. We use procedures described in Appendix B (Sect. B.2) and in Peixinho et al. (2004) to deal with these effects to obtain Spearman-ρ values that take into account the error bars.

We call correlations “strong” when . Efron & Tibshirani (1993) analyze the evidence for correlation with the follow criteria (where CL is the confidence level defined in Appendix B): (): very strong evidence of correlation; and ( 0.003), the “standard” 3σ significance: clear correlation. We restrict ourselves to the cases with ( 0.005), i.e. at least 2.8σ confidence.

5.4.2. Correlations analysis

Based on the previously described methods, the strongest and/or more interesting correlations are shown in Table 6. These are given for the whole sample (including or excluding Eris), and for SDOs and detached objects separately. Before discussing them, we stress that results must be taken with care because of small number statistics. In general, more significant correlations are found for SDOs than for detached objects.

“Clear correlations” (99.7%) in Table 6 are seen between D and H (negative correlation, 3.0σ), p_V and H (negative, 4.2σ), D and r_h (positive, 4.9σ), p_V and r_h (positive, 3.8σ). The first two anticorrelations are a simple consequence of Eq. (1), with no physical implication. The positive correlations between D and r_h, and between p_V and r_h, are most probably biases due to object selection criterion and discovery bias, as discussed in Sect. 2.

Additional “clear” or “very strong evidence” (99%) for positive correlations are found between D and q (3.4σ) and p_v and q (2.9σ) (Fig. 4). There is no a priori relationship between diameter, albedo and perihelion distance, hence these correlations could plausibly be real, as discussed below. However, they could also reflect biases due to the correlation between q and r_h (at 3.9σ for our sample, see Table 6), i.e. the discovery bias discussed in Sect. 2. On the other hand, it is noteworthy that no high albedo objects (e.g. p_v > 10%) are seen at small heliocentric or perihelion distance (≤ 30 AU), which should not be a bias. Note that when SDOs and detached objects are treated separately, the D vs. q correlation remains significant at 2σ only, while the p_v and q correlation essentially vanishes. Correlations also decrease (marginally) if Eris is excluded from the sample.

Based on Spitzer measurements of about 50 objects, Stansberry et al. (2008) noted the apparent positive correlation between p_v and q (3.5σ in their case). While recognizing the above discovery bias, they proposed it to be real, offering two explanations. First, objects near the Sun tend to experience higher temperatures, hence may be subject to increased sublimation of their intrinsically bright ices. Conversely, increased UV-photolysis and solar wind radiolysis associated with the low heliocentric/perihelion distances may preferentially darken the surfaces of these objects. It seems more difficult to find a physical explanation for an hypothetical correlation between diameter and perihelion distance, especially because the current orbits of the SDOs/Detached objects tell us little if anything about the formation region. This apparent correlation is most probably a bias, as discussed above. Note also that unlike Stansberry et al. (2008) we do not find any (even tentative) correlation between p_v and the semi-major axis a (ρ = 0.21, 0.8σ, see Table 6).

We find a very strong correlation between albedo and diameter (p_V vs. D) when SDOs are included, in particular for the whole sample at 2.9σ significance (see Fig. 4). This result was obtained by Stansberry et al. (2008) for KBOs (defined there as any object but Centaurs) but was not robust against small changes in the classification. This correlation implies a trend for larger objects having higher albedos. An explanation could invoke the fact that larger objects can retain bright ices (released either by collision or intrinsic activity) more easily than small objects do. Eris follows this tendency, although its extremely high albedo puts it as an outlier in the plot of Fig. 4. An important aspect here is that there are no small objects (≤300 km) with high albedo (≥10%). This cannot be a bias, as such objects would have been targetted in our sample, yet poorly detected. Note that if objects with albedos  ≥10% are to be associated with the presence of surface ices, then such ices should be only moderately volatile (i.e. CH₄ or less), as volatile retention models (Schaller & Brown 2007b; Brown et al. 2011) show that only the largest (≥1000−1500 km) bodies can retain the most volatile ices such as N₂ or CO over the solar system age.

Finally, we do not detect any correlation between diameter/albedo and color, or inclination, nor between the beaming factor η with any parameter (though the sample is limited to 7 objects in the latter case).

In summary, in spite of the relatively small sample covered in this paper, we find three possible correlations (i.e. p_V vs. D, p_V vs. q, and D vs. q), of which at least the first one – large objects tend to have brighter surfaces than small objects – is real and may be physically explained. More subtle correlations may be found, but with lower significance, and larger datasets will be required for confirmation.

6. Summary

We have determined albedos/sizes of 8 SDOs and 7 detached objects using Herschel-PACS fluxes measurements at 70, 100 and 160 μm and Spitzer-MIPS fluxes at 23.68 and 71.42 μm when available. The main results, which include nine first determinations, are summarized hereafter:

• Diameters range from 100 to 2400 km. The overall diameter distribution in our sample is remarkably uniform, but apparently different for SDOs vs. detached objects. This difference might however just reflect a sample/discovery bias, as our sample favors large objects at large heliocentric distances, and detached objects are typically discovered at larger distances than SDOs.

• V-geometric albedos range from 3.8 to 84.5%. Unweighted mean albedo for the whole sample is 11.2% and 6.9/17.0% for the SDOs/Detached (excluding Eris) respectively. There is a lack of albedos  >20% with a peak in the [0–5%] range for SDOs. Detached objects have albedos  >5%, with a peak in the [10–20%] range.

• We re-determine bulk densities of three binary systems, based on our estimated sizes and published masses: Ceto/Phorcys ( g cm^-3), Typhon/Echidna ( g cm^-3), and Eris/Dysnomia ( g cm^-3).

• Beaming factor (η), determined from thermal fits of 7 objects, range from 1.02 to 1.48, with a weighted mean obtained of 1.14  ±  0.15, consistent with the mean value derived from Spitzer, 1.20  ±  0.35 (Stansberry et al. 2008). These low η values are compatible with low thermal inertias, presumably due to high porosity of regolith-like surfaces.

• We obtain for the first time the albedo and size of the very large SDO 2007 OR₁₀: D = 1280    ±    210 km, %. Along with its surface composition, this indicates that the object is Quaoar-like.

• Our size/albedo estimations for the Eris/Dysnomia system (D = 2454    ±    117 km, p_V = 84.5    ±    8.8%) are compatible within error bars with occultation results by Sicardy et al. (2011), but further analysis is needed.

• We find a significant correlation between albedo and diameter (more reflective objects being bigger), probably due to the fact that large objects can retain bright ices more easily than small objects. Stansberry et al. (2008) obtained a similar result for a larger sample of TNOs, but they noted that the correlation is sensitive to changes in the classification used (i.e. MPC vs. DES; Elliot et al. 2005).

• We also find positive correlations between albedo and perihelion distance, and between diameter and perihelion distance (brighter and bigger objects having larger perihelia). The first correlation has been seen before based on Spitzer measurements (Stansberry et al. 2008) and interpreted in terms of increased ice sublimation and/or increased space weathering at low heliocentric distances. It seems more difficult to find a physical explanation for the second correlation.

Online material

Appendix A: Data reduction

In this appendix, we describe the reduction process applied to the PACS data in order to obtain what we call “single” and “combined” maps.

A.1. General PACS photometer data reduction

As explained in the main text (see Sect. 3.1) the data reduction of the maps was performed within HIPE (version 6.0.2044). We modified and adapted the standard HIPE scripts for the needs of our programme. Thus, the PACS photometer data is reduced from Level-0 to Level-2 using a modified version of the standard mini-scanmap pipeline. A major difference of our reduction chain compared with the standard pipeline is the application of frame selection on scan speed rather than on building block identifier (also known as BBID, which is a number that identifies the building block – a consistent part – of an observation). It has been widely tested that the application of the scan speed selection increases the number of useable frames, and finally increases the SNR of the final maps by 10–30% (depending on the PACS observation requests (AOR) setup and on the band used), which is especially important for our faint targets.

Since our targets move slowly (a few ″/h), and the total observation time per scanning direction in one visit is between 10−25 min we do not correct for the proper motion. Performing this kind of correction would smear the background, and would make it impossible to make a correct background subtraction in the later stages (see Sect. A.2).

We use a two-stage high-pass filtering procedure, as in the standard mini scan-map pipeline. In the first stage a “naive” map is created. This map is used in the second stage to mask the neighborhood of the source and the high flux pixels. In the second stage the vicinity of the source has been masked (2 times the FWHM of the beam of the actual band), as well as all pixels with fluxes above 2.5 ×  the standard deviation of the map flux values. We used the standard high-pass filter width parameters of 15, 20 and 35. The final maps are created by the PhotProject() task, using the default dropsize parameter.

These first steps lead to the production of one “single” map per visit, filter, and scan direction (i.e. in total 8 maps per object in the red, and 4 maps in the blue/green created using the MosaicTask() in HIPE). The sampling of the single maps generated with these scripts are 1.1″/pixel, 1.4″/pixel, and 2.1″/pixel for the blue (70 μm), green (100 μm) and red (160 μm) bands respectively.

A.2. Building final maps for photometry

We then use these single maps to generate the final “combined” maps on which the photometry will be performed. Essentially, to generate the combined maps, we determine the background map, subtract it from each single map, and finally co-add all the background-subtracted maps. A detailed description of this process is as follows:

• i)

  We have initially a set of 8 (red) or 4 (blue and green) single maps taken on different dates. Considering the red maps, let us call t1 to t4 (t5 to t8) the dates corresponding to the first (second) visit. t1-t4 and t5-t8 are separated by typically 1–2 days so that the target motion (typically 30–50″ as described in Sect. 2 of the main text) produces a significant change in the observed sky field (see Fig. A.1).

  Fig. A.1 Left four panels: Typhon red band 1st visit single maps (t1 to t4). Right four panels: Typhon 2nd visit single maps (t5 to t8). t1, t3, t5, and t8 were observed with a scan angle of 110°, and t2, t4, t6, and t7 with a scan angle of 70° with respect to the detector array. The green circle marks the object position.

• ii)

  A background map is generated using the single maps. To do this we “mask” the target in each of the 8 (4) single maps and co-add the maps in the sky coordinate system. This step produces a background map with high SNR without the target (see the right panel in Fig. A.3 for the Typhon case at 160 micron).

• iii)

  The background map is subtracted from the single maps (Fig. A.1), producing 8 (4) single maps with background removed. We call these images “background-subtracted maps” (Fig. A.2).

  Fig. A.2 Same as Fig. A.1 for the background-subtracted maps.

• iv)

  Finally the background-subtracted maps are co-added in the target frame (center panel of Fig. A.3), producing the final combined map on which photometry is performed.

Note that in step ii) before combining the images of the two epochs, we need to calculate the optimal matching offset between the coordinate systems of the two maps since, when we use the original World Coordinate System (WCS) information stored in the file headers, artifacts of bright sources remain in the background subtracted image. The best offset position was found in the following way: we added a small offset – in steps of 0.5″ in both the RA and Dec directions – to the appropriate original WCS parameters in the file headers, and performed the background subtraction. For each of these maps the standard deviation was calculated for pixels with a coverage value of at least 35% of the maximum coverage of the map. The offset values related to the lowest standard deviation correspond to the best matching of the two backgrounds. After this procedure the artifacts of bright sources (e.g. pairs of positive-negative spots) completely disappeared (or, at least, were below detection), and the noise in the background decreased notably. This method was performed on the 160 μm band measurements, where the effect of the background is the largest. However, it has been checked in a series of measurements that the same offsets can be applied to the 70 and 100 μm images as well.

An alternate, simpler method, to obtain final maps is to co-add directly the original single maps (i.e. not background-corrected) in the target frame. We call this the simple coaddition method (left panel of Fig. A.3). This method is obviously less optimal than the previous one in terms of SNR, but provides a useful test to demonstrate that the background subtraction is not introducing any spurious effects.

A more detailed description about the whole PACS data reduction process will be published in Kiss et al. (in prep.-a).

Appendix B: Statistical aspects

B.1. “Rescaled error bar” approach

As indicated in the main text (see Sect. 4.3), error bars on the fitted parameters (diameter, albedo and beaming factor) are determined through a Monte-Carlo approach which essentially consists of randomly building synthetic datasets using uncertainties in the measured fluxes. However, when the measurement errors bars are too small, or similarly when the model is not completely adequate (both cases being indicated by a “poor model fit”), the described Monte-Carlo approach will underestimate the uncertainties on the solution parameters. In this case, we adopted a “rescaled error bar” approach described hereafter.

Quantitatively, the fit quality is determined by the value of the reduced χ²: (B.1)where O, M and σ are the observed, modelled and error bar flux values, and ν is the degree of freedom (ν = N − 1 for fixed-η, and ν = N − 2 for floating-η, where N is the number of thermal wavelengths available). While indicates a good fit, indicates a poor fit. The idea is therefore to empirically “rescale” (uniformly, in the lack of a better choice) all error bars σ by  before running the Monte-Carlo approach. This method leads to much more realistic error bars on the solution parameters. However, the range of for which this approach is warranted depends on the number of data points, as for few observations does not have to be closely equal to 1 to indicate a good fit. Specifically, the variance of the variance estimator of a distribution of N points picked from a Gaussian distribution with dispersion = 1 is itself a Gaussian with a mean equal to unity and a standard deviation equal to . If we had for example N = 1000 observations, 68.2% of the good fits have a reduced χ² between 0.955 and 1.045, i.e. is strongly constrained to 1; for N = 5 (resp. 3) observations in contrast, this range is much broader, 0.368–1.632 (resp. 0.184–1.816). This means for example that for 3 data points (PACS only), a reduced χ² of e.g. 1.6 is acceptable. We considered these N-dependent limits as the thresholds beyond which the fits are statistically bad and the error bars need to be rescaled. Note finally that the rescaling approach is merely an “operative” method to avoid unrealistically low error bars on the fitted diameter, albedo and beaming factor, and that error bars shown in Figs. 1 and 2 of the main text are the original, not rescaled measurement uncertainties.

B.2. Spearman-rank correlation with error bars

As described in the main text the Spearman-ρ correlation test is distribution-free and less sensitive to outliers than most other methods but, like others, treats data-points as “exact” and does not take into account their possible error bars. Any variations in the measured data points within their error bars, may change the correlation coefficient, however. Furthermore, each correlation coefficient has its own confidence interval, which depends on the number of data points and on the magnitude of the correlation. To account for these effects we used the following three procedures:

• 1.

  To estimate the most probable correlation coefficient given the error bars, we generate 1000 samples of data points, building each synthetic dataset from its associated probability function. When errors bars are symmetric this probability function is considered to be Gaussian and the error bar correspond to one standard deviation. When errors bars are asymmetric, we fit a lognormal probability function such that the observed value corresponds to the distribution’s mode, and the interval between the upper and lower limits of the observation corresponds to the shortest interval containing 68.2% of all possible values. The resulting distribution of correlation values is not Gaussian, but its Fisher-transform z = argtanh(ρ) is. So, we can determine the most probable correlation value,  ⟨ ρ ⟩  from our Monte-Carlo simulations, and its upper and lower limits (+ σ_MC/−σ_MC), which are not necessarily symmetric after the reconversion using ρ = tanh(z) (see Peixinho et al. 2004, for more details). It is noticeable (and expected) that observational error bars “degrade” the correlation. The approximate significance p-value, or confidence level (CL), of ρ can be computed from , which follows a t-Student distribution with n − 2 degrees of freedom. However, when n ≲ 15, using t provides a slight overestimation of the significance. We have compared the previous calculation of the significance with exact values tabulated by Ramsey (1989) and computed adjustments required to obtain a more accurate approximation of the ρ’s and p-values in the case of low n.

• 2.

  The confidence level of a given ρ results from testing the hypothesis “the sample is correlated” against the hypothesis “the sample correlation is zero”. Knowing the confidence interval (CI) within which the correlation ρ value of the parent population lies may be more informative. For example, suppose we have ρ = 0.7 with a 3σ confidence level; we would conclude that the sample is correlated. But if the 68% confidence interval of the correlation were say  [0.3,0.9] , we would be unsure if the correlation was very strong or rather weak. To estimate the shortest (+ σ_B/−σ_B) interval containing 68.2% of the population’s ρ (i.e. the equivalent to the Gaussian 1σ interval), we used 1000 bootstrap extractions from the data-sample computing this range as we did for the previous item (e.g. Efron & Tibshirani 1993).

• 3.

  Finally, to be fully correct, one should perform the 1000 bootstraps on each one of our 1000 Monte-Carlo simulations to obtain the true 68% confidence interval of  ⟨ ρ ⟩ , a computationnally heavy task. Fortunately, the combination of the two effects can be obtained by quadratically adding the standard deviations of the Gaussian distributions of the Fisher-transform of the Monte-Carlo simulations and the bootstraps, i.e. , which, after reconversion to ρ, will give our final best estimate for the 68% confidence interval noted as .

---

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. P. Santos-Sanz would like to acknowledge finantial support by the Centre National de la Recherche Scientifique (CNRS), also acknowledge Brett Gladman for his help with the dynamical classification of TNOs, and Daniel Hestroffer and Will Grundy for his help with the bulk density estimation of binary systems. M. Mommert acknowledges support through the DFG Special Priority Program 1385. Part of this work was supported by the German Deutsches Zentrum für Luft- und Raumfahrt, DLR project numbers 50 OR 1108. We are grateful to Danielle Tibi for her unvaluable expertise on statistical aspects (Sect. B.1). J. L. Ortiz acknowledges the spanish grants AYA2008-06202-C03-01, AYA2011-30106-C02-01 and 2007-FQM2998. R. Duffard acknowledges financial support from the MICINN (contract Ramon y Cajal).

References

All Tables

All Figures