The Density of Mid-sized Kuiper Belt Objects from ALMA Thermal Observations

For small solar system bodies like Kuiper Belt objects (KBOs) about which little is known, density can be one of the key characteristics that informs our understanding of bulk composition, formation, and physical state. For most small bodies, measurement of density requires separate accurate measurement of the mass and of the volume of the system. While mass can often be measured to high accuracy by fitting a Keplerian orbit to the motion of a small satellite, measurement of the size of the body is often more difficult and less accurate. While stellar occultations have given extremely accurate estimates of a limited number of larger bodies (Sicardy et al. 2006, 2011, 2016; Braga-Ribas et al. 2013; Brown 2013b), most size estimates for these small bodies rely on the much less accurate technique of thermal radiometry.

In thermal radiometry, the thermal emission of a body at a known distance from the Sun is used to estimate the emitting area—thus the diameter—of the body. To date, most measurements of thermal emission have been made from the Spitzer Space Telescope (at 24 and 70 μm; Stansberry et al. 2008, pp. 161–179) and with the Herschel Space Telescope (at 70–500 μm; Müller et al. 2010). Such measurements have the advantage of spanning both the Wien and Rayleigh–Jeans side of the thermal emission, which typically peaks around ∼100 μm for bodies at these temperatures. The disadvantage of these measurements, however, is the moderately low signal-to-noise ratio (S/N) of the thermal emission at these wavelengths with these telescopes and also the sensitivity of the size measurements to model assumptions such as the thermal properties, spin period, and pole position of the objects

The Atacama Large Millimeter Array (ALMA) radio observatory provides the potential for thermal observations of KBOs with a complementary set of advantages and disadvantages. With the massive collecting area of ALMA, thermal observations can be obtained at significantly higher S/N. Unfortunately, these observations are at longer wavelength and thus exclusively on the Rayleigh–Jeans tail of the thermal emission, preventing an accurate constraint on the temperature and temperature distribution of the KBO surface. Interestingly, however, the Rayleigh–Jeans portion of the emission is only linearly dependent on temperature, instead of exponentially. More uncertainly, there is the possibility that at these long wavelengths emissivity-like effects depress the thermal emission in unknown ways (Müller & Lagerros 1998; Fornasier et al. 2013; Lellouch et al. 2016), making prediction of thermal flux densities (or of radii from thermal measurements) unreliable. We note that a similar depression has been seen in the flux densities of asteroids at long wavelengths; this is not due to a true emissivity effect, but rather the properties of the subsurfaces of these bodies (i.e., when a proper radiative transfer model is used with realistic surface and subsurface properties, no depressed emissivity is needed to fit the observations (Keihm et al. 2013)). To date, with no available Spitzer/Herschel measurements along with high signal-to-noise longer wavelength measurements, we cannot assess the effect of such measurements on our ability to measure size—and thus density—of KBOs.

We explore the utility of ALMA measurements of KBO size for a carefully selected group of objects. As noted in Brown (2013a), the apparent increase in density with diameter of known KBOs is difficult to explain in conventional accretion scenarios. The smallest bodies have densities below 1 g cm⁻³, suggesting both significant porosity and significant ice fraction, while the largest bodies have densities greater than 2 g cm⁻³, implying a substantial rock fraction. The difficulty of creating the large bodies by accretion of the small bodies could be overcome if the small bodies have similar rock fractions to the larger bodies but are significantly more porous. The mid-sized KBOs provide the key to answering this question. At diameters greater than about 600 km little porosity can be sustained inside of KBOs (Yasui & Arakawa 2010), so the density measured comes close to reflecting a true ice-to-rock ratio. Brown (2013a) found that 2002 UX25, with a diameter of 692 ± 23 km, had a density of 0.82 ± 0.11 g cm⁻³, suggesting that, indeed, the smallest bodies were incapable of accreting to form the larger bodies. Here we attempt to improve the density measurements of a sample of four mid-sized KBOs in order to more carefully explore this transition from small-to-large bodies and to determine if the results found from 2002 UX25 carry through to this larger sample.

All observations were undertaken with the 12 m array of the ALMA. This synthesis array is a collection of radio antennas, each 12 m in diameter, spread out on the Altiplano in the high northern Chilean Andes. Each of the pairs of antennas acts as a two-element interferometer, and the combination of all of these individual interferometers allows for the reconstruction of the full sky brightness distribution, in both dimensions (Thompson et al. 2001).

ALMA is tunable in 7 discrete frequency bands, from ∼90 to ∼950 GHz. All observations in this paper were taken in Bands 6 and 7, near 230 and 350 GHz, in the "continuum" (or "TDM") mode, with the standard frequency tunings. For band 6, this yields four spectral windows in the frequency ranges: 220–222 GHz; 222–224 GHz; 236–238 GHz; and 238–240 GHz. For band 7 the frequency ranges are: 337–339 GHz; 339–341 GHz; 349–351 GHz; and 351–353 GHz. In the final data analysis we average over the entire frequency range in both bands, and use 230 GHz and 345 GHz as the effective frequencies in our modeling.

All of these observations are in dual-linear polarization; in the end we combine these into a measurement of Stokes I. While we expect polarized emission from the surfaces, it is weak and in an unresolved image averages to zero.

Table 1 shows the information for all of our observations; most of which were executed in 2013 November (with a single observation in 2014 March). The number of antennas included in the observation was between 24 and 35, but in all but the 2014 March observation of Quaoar at band 7 there were three antennas that were significantly distant from the others and provided no useful signal so had to be removed from the processing. The array was mostly in the C32-3 configuration, which yields a resolution of ∼09 in band 6 and ∼06 in band 7—the actual resolution for each of the observations is shown in Table 1. Each observation was of order 1 hr in duration, including all calibration overheads, which resulted in roughly 20–30 minutes on source. Neptune and Titan were used as the absolute flux density scale calibrators for all observations (Butler 2012). Nearby point-like radio source calibrators were used to calibrate the phase of the atmosphere and antennas as a function of time.

Initial calibration of the data was provided by the ALMA observatory, and is done in the CASA reduction package via the ALMA pipeline (Muders et al. 2014). The actual measured quantity of a complex interferometer like ALMA is a sampling of the complex visibility function at the positions of the baselines between each of its antennas. The visibility function is the two-dimensional Fourier transform of the sky brightness distribution. The individual samples of the visibility function are referred to as visibilities, and are complex quantities (real and imaginary, or amplitude and phase). After the initial calibration, the data product was a set of visibilities for each of the observing dates.

At this point we exported the data from CASA and continued the data reduction in the AIPS package (www.aips.nrao.edu/CookHTML/CookBook.html). We imaged each of the objects at each band with natural weighting. The result is shown in Figures 1 (band 7) and 2 (band 6). The final step of the data analysis was to estimate the observed flux density for each body in both bands. We obtained this value in a number of ways, to check for consistency: flux density in the image; flux density in the CLEAN components; fitting a Gaussian in the image; and fitting the visibilities directly. We found good agreement (better than the final 1σ uncertainty) for all of these techniques. We take the visibility fit value as the best value, as it avoids the biases of fitting in the image plane in the presence of correlated noise (Greisen 2004). As a byproduct of the fitting process, an offset from the phase center is also derived. Table 2 shows the final fitted flux densities and offsets for our observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Background

3.2. Multi-parameter Markov Chain Monte Carlo Thermal Modeling

We model thermal emission from these KBOs using a standard thermal model with a beaming factor which is adjusted to account for the combined effects of thermal inertia, pole position, and surface roughness. Our thermal model takes as input the radius, R, beaming factor, η, and Bond albedo A. The surface temperature at angle Θ away from the sub-solar point is calculated as

$$\begin{eqnarray*}&&T={\left[\displaystyle \frac{S\cos ({\rm{\Theta }})(1-A)}{\sigma \epsilon \eta }\right]}^{1/4},\end{eqnarray*}$$

where is the bolometric emissivity—discussed below—σ is the Stefan–Boltzmann constant, and S is the solar insolation at the distance of the object. The thermal emission from the observer-facing area is then integrated to calculate the total emission at each wavelength. The possible range of bolometric emissivity, , for KBOs, is unclear. Bolometric emissivities of Pluto and Charon have been estimated to be between 0.83 and 0.93 (Lellouch et al. 2016). In well-measured asteroids, emissivities vary from about 0.8 to 0.9 (Müller & Lagerros 1998). We will take values from 0.8 to 1.0 to encompass the range of possible uncertainty and allow our emissivity to vary between these values. Typical KBO thermal models assume a fixed emissivity of 0.9.

The model also calculates the expected absolute visible magnitude, H_V, of the object from its radius and visible geometric albedo, p_V, as

$$\begin{eqnarray*}&&{H}_{v}=-5{\mathrm{log}}_{10}({{Dp}}_{v}^{1/2}/1330),\end{eqnarray*}$$

where D is the diameter, in kilometers. The visible geometric albedo, p_v, which determines the zero-phase reflected sunlight, is connected to the Bond albedo, which determines the energy absorbed, through the phase integral, q, as A = qp_v. No phase integrals have been measured for KBOs, as high phase observations have been unavailable, but Brucker et al. (2009) show that the icy Saturnian and Uranian satellites roughly follow a linear function given by q = 0.336p_v + 0.479. Significant outliers occur, however, with Phoebe—perhaps a good analog for darker KBOs—being a factor of two below this value, and Europa being a factor of two above. We take this factor of two below or above variation to represent the uncertainty in our knowledge of the phase integral and allow our phase integral to vary by this factor from the Brucker fit.

The models are compared to the observations of the flux density as a function of frequency and to the measured absolute visible magnitude using a Markov chain Monte Carlo (MCMC) model in which the free parameters are diameter, albedo, and beaming factor, with variations in emissivity and phase integral and the uncertainties in measured absolute magnitude also included as nuisance parameters. The individual thermal observations were obtained when the object was at different heliocentric and geocentric distances, so, although the differences are small, each flux density is modeled for the individual distance that the observation was made at. We use the Python package emcee (Foreman-Mackey et al. 2013), which provides a convenient and parallel implementation of the Hou et al. (2012) affine invariant ensemble sampler for MCMC. In each of the modeling cases that we discuss below, we assign uniform priors to the parameters and run an ensemble of 100 chains through 10⁴ steps after a 10³ step initialization ("burn-in") period. The chains converge with no obvious memory of initial condition. We examine marginalized distributions of all of the parameters. The distribution of diameter is nearly Gaussian, thus we report the median and the middle 68.2% range to represent the 1σ uncertainties.

First, we explore the effects of the model uncertainties on the previously estimated diameters of our four target KBOs. Fornasier et al. (2013) suggest that for many objects effective emissivity drops perhaps longward of 350 μm, so we limit the initial modeling to shorter wavelengths to explore this effect using only the data from Spitzer and the shorter Herschel bands. Figure 3 shows samples of the thermal model from the Markov chain for the four targets, while Table 3 gives derived diameters and albedos. Note that Fornasier et al. (2013) report detections that are less than 1σ as upper limits at the 1σ flux level, which makes these fluxes impossible to use in any statistical way; as we have no additional knowledge of the actual measured flux level, we run a simple Monte Carlo model to estimate that the correct 1σ upper limit to a randomly distributed signal is measured to be somewhere between 0 and 1σ is 1.86σ; we accordingly use those values for the limits, though we do not incorporate them into the fits.

Zoom In Zoom Out Reset image size

The results from Table 3 show that while using a more realistic estimate of model and parameter uncertainties does not change the best fit values of the diameters and albedos, it does increase the uncertainties in some of these values. The uncertainties in the diameters are increased, on average, by a factor of about 25%, albeit with considerable variation. While this value is not extreme, it does lead to a 75% increase in the uncertainty of the measurement of the volume and thus the density of these objects. We find that the uncertainty increase is dominated by the expanded range of parameters such as the emissivity and phase integral. The models do not give substantially different best-fit results, but the range of diameters over which adequate fits can be obtained is clearly higher.

While appropriate use of the upper limits makes it less clear that the longest Herschel wavelengths suffer from depressed emission due to decreasing effective emissivities, Figure 4, which shows an expanded view of these Spitzer/Herschel-only fits to the ALMA observations, shows a clear discrepancy. The ALMA flux densities are systematically about 30% lower than predicted from the Spitzer/Herschel observations. With the high signal-to-noise ratio of the ALMA observations, these discrepancies are highly significant. In the next section we explore whether these higher signal-to-noise observations can improve the estimates of the diameters.

Zoom In Zoom Out Reset image size

With the significantly smaller uncertainties, the ALMA measurements have the possibility of significantly improving the diameter uncertainties. We use our MCMC model to fit the full suite of Spitzer, Herschel, and ALMA data (with the exception of the 350 and 500 μm Herschel data, which we still exclude because of the uncertainty about the uncertainties). Figure 5 shows the fits to this full data set. In most cases, the small uncertainties of the ALMA results drive the model fits to significantly under predict the Spitzer and Herschel data near the peak of the thermal emission, while still over predicting the ALMA flux densities. Best fit diameters are decreased by ∼15%, and while the formal uncertainties decrease by as much as 25%, the actual fits appear quite poor.

Zoom In Zoom Out Reset image size

We can conceive of only two likely explanations for this effect. Either there is a serious difference in calibration between ALMA and the infrared space telescopes, or effective emissivity effects at longer wavelengths are indeed suppressing the longer wavelength thermal emission on these bodies.

A calibration uncertainty of this magnitude in the ALMA data is unlikely, given the care with which these data are calibrated. Flux densities measured by ALMA are specified to have an accuracy of ±5%, and all tests of this claim to date have verified this accuracy. There have been some recent questions regarding the use of monitored quasars and asteroids to set the flux density scale for ALMA observations, but in our case only Neptune and Titan were used for this, and confidence is high that the models for these two bodies are good (Butler 2012).

A suppressed millimeter emissivity seems the most likely cause of the low ALMA flux densities. Fornasier et al. (2013) discuss decreasing emissivity at 350 and 500 μm in the Herschel data, but with the large uncertainties at these wavelengths the true effect is difficult to determine. With the ALMA data, the discrepancy becomes unmistakable. In Figure 6, we use the Spitzer/Herschel only fits of the previous section to determine the effective emissivity as a function of wavelength.

Zoom In Zoom Out Reset image size

Note that Lellouch et al. (2016) have also found a depressed emissivity at these long wavelengths in observations of the Pluto/Charon system. They interpret this depression as due to a combination of surface and subsurface dielectric and particle/volume scattering rather than from surface roughness effects (which have often been invoked for these and other bodies at the shorter thermal emission wavelengths, and as described above are actually part of the genesis of the beaming factor). A similar depression on the four KBOs we have observed here may be due to the same effect.

Some satellites are potentially large enough to have a significant effect on the total thermal emission. If the primary and the satellite have identical surface properties (and pole positions), their emission spectra are identical and the individual contributions simply scale as the squares of their diameters. In this simplest case the diameter of the primary can be obtained from the equivalent diameter, D_e, calculated above, and the difference in magnitude between the primary and satellite ΔV, as ${D}^{2}={D}_{e}^{2}{(1+{10}^{-{\rm{\Delta }}V/2.5})}^{-1}$.

KBO primaries and their satellites have highly correlated surface colors (Benecchi et al. 2009), leading to the reasonable expectation that the surfaces are similar and have similar enough properties that this simplest case prevails. The Orcus–Vanth system clearly deviates from this simple expectation. Orcus has significant water ice absorption that Vanth is lacking, and they have differing optical colors (Brown et al. 2010). The derived system albedo of 0.23 ± 0.02 is high for typical KBOs but not unusual for large ones with visible water ice. All of these considerations strongly suggest that Vanth is darker than Orcus. In such a case Vanth will be warmer than Orcus and a simple scaling cannot be used to estimate their sizes. We run an additional MCMC model to simultaneously fit the emission from a potentially darker Vanth. In this model we force Orcus and Vanth to have the same beaming factor, bolometric emissivity, and phase integral; these assumptions could be suspect if the two bodies have very different surfaces. We fix Vanth to be 2.61 mag fainter than Orcus (Brown et al. 2010), and add the constraint that the albedo of Vanth must be less than or equal to the albedo of Orcus. The MCMC model fits these models to all of the data out to 250 μm, appropriately accounting for thermal emission from, for example, a large darker warmer Vanth, or a small bright colder Vanth. The median of the marginalized posterior distribution of the diameter is ${885}_{-80}^{+55}\,\mathrm{km}$ for Orcus with an albedo of ${0.25}_{-.03}^{+.05}$ and a diameter of ${370}_{-70}^{+160}\,\mathrm{km}$ for Vanth and an albedo of 0.13 ± 0.06. This model fit should not be considered a detection of Vanth, but merely a statement of the best limits to the size of Vanth with the stated assumptions and observations. The retrieved albedo of Vanth remains higher than typical values for objects of that size, and, indeed, a significant tail to the distribution extends to albedos as low as 0.04 (below that the thermal emission from such a large Vanth becomes too significant), with corresponding larger diameters for Vanth. We will retain our formal error bars as calculated from the MCMC model but not discount the possibility that Vanth could still be darker and larger.

Our main interest in calculating diameters and more rigorous uncertainties for this set of KBOs is to calculate the densities of these mid-sized objects. One final uncertainty that we must consider is the possibility of a difference in the density of the primary and the secondary. The two KBOs that have estimates for satellite densities show that the range can be extreme: Charon has a density nearly equal to that of Pluto (Stern et al. 2015), while the two satellites of Haumea have densities several times smaller than that of Haumea (Ragozzine & Brown 2009). We will make the assumption that satellite densities can range from 0.5 g cm⁻³ up to the density of the primary. For Salacia, Quaoar, and 2002 UX25, we simply calculate the range of satellite masses, subtract that from the system mass, and determine the density of the primary. For Orcus, we take the Orcus and Vanth sizes from our Markov chain and randomly assign Vanth a density between 0.5 g cm⁻³ and the density of Orcus. We arrive at an Orcus density of ${1.65}_{-0.24}^{+0.34}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. The derived sizes and densities of all objects are shown in Table 4.

Table 4.  Derived Diameters and Densities

Body	ΔV	Primary Diameter	Satellite Diameter	Primary Density	Occultation Diameter
	(mag)	(km)	(km)	(g cm−3)	(km)
2002 UX25	2.28 ± 0.06	659 ± 38	230 ± 19	0.80 ± 0.13	⋯
Orcus	2.54 ± 0.01	${885}_{-80}^{+55}$	${370}_{-70}^{+170}$	${1.65}_{-.24}^{+.35}$	⋯
Salacia	2.37 ± 0.06	866 ± 37	290 ± 21	1.26 ± 0.16	⋯
Quaoar	5.6 ± 0.2	1079 ± 50	82 ± 17	2.13 ± 0.29	1110 ± 5a

Note.

^aBraga-Ribas et al. (2013).

Download table as:  ASCIITypeset image

The high sensitivity of ALMA at millimeter wavelengths allows precise measurements of thermal emission from objects in the outer solar system. Unfortunately, the lack of independent knowledge of millimeter emissivities of these objects prevents these less uncertain measurements from decreasing the uncertainties of diameter measurements of these objects. Indeed, when parameter uncertainties are more carefully included, we find that the current estimates of the uncertainties in the sizes of KBOs are approximately a factor of 25% too small.

Currently, no observatories are capable of observing the ∼100 μm thermal peak of KBOs, so any current size or albedo measurements will have to rely on these longer wavelength observations. To explore how well ALMA-only measurements would recover the sizes of our four observed objects, we assume that the effective emissivity of each object is the average of that measured from the other three objects. We then rerun our thermal emission MCMC fitting only the ALMA data scaled by this effective emissivity (but retaining the 0.8–1.0 bolometric emissivity for the emission peak). We find diameters as follows: 2002 UX25 is estimated to be ${742}_{-109}^{+76}\,\mathrm{km}$, Orcus is ${1075}_{-156}^{+121}\,\mathrm{km}$, Salacia is ${900}_{-140}^{+95}\,\mathrm{km}$, and Quaoar is ${1057}_{-168}^{+107}\,\mathrm{km}$. In all cases the ALMA-only diameters are within 1-σ of the Spitzer/Herschel diameters, albeit with uncertainties two to three times larger.

As a second check, we derive a diameter for Charon, which has been measured to have a flux of 7.0 ± 0.07 mJy at 840 μm (Butler et al. 2015). We use our ALMA-only model with an assumed emissivity of 0.685—the average of the values for our four KBOs—and find a diameter of 1355 ± 110 km, which is within 1.3σ of the measured value of 1212 ± 6 km (Stern et al. 2015).

The larger uncertainties in KBO diameters inferred from the analysis here renders the original goal—better constraining the behavior of density versus size for these objects—impossible. Assuming that most KBO diameters measured with radiometry (rather than occultations) have similarly underestimated uncertainties does not substantially change the interpretation of Brown (2013a), but it makes further progress difficult.

The densities of mid-sized KBOs remain a key constraint for understanding the accretional history of the solar system. Because of the moderate model and parameter uncertainties associated with diameter measurement from thermal radiometry combined with a factor-of-three greater effect on volume and density measurements, we conclude that this technique is unlikely to yield the precision necessary to further constrain these densities. A high priority should be placed on obtaining occultation measurements of satellite-bearing mid-sized objects.