Photometry of Outer Solar System Objects from the Dark Energy Survey. I. Photometric Methods, Light-curve Distributions, and Trans-Neptunian Binaries

The modification for a binary TNO requires changing Equation (2) to include a second point-source term:

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{M}}}_{{ij}}^{\nu }=\displaystyle \sum _{u,v}{\boldsymbol{P}}{\boldsymbol{S}}{{\boldsymbol{F}}}_{\nu }(i-u,j-v){{\boldsymbol{P}}}_{{uv}}+{{\boldsymbol{b}}}_{\nu }\\ \,+{f}_{1}{\boldsymbol{P}}{\boldsymbol{S}}{{\boldsymbol{F}}}_{\nu }(i+{\rm{\Delta }}{x}_{1},j+{\rm{\Delta }}{y}_{1})\\ \,+{f}_{2}{\boldsymbol{P}}{\boldsymbol{S}}{{\boldsymbol{F}}}_{\nu }(i+{\rm{\Delta }}{x}_{2},j+{\rm{\Delta }}{y}_{2}),\end{array}\end{eqnarray} \tag{ 5 }$$

where the indices 1 and 2 correspond to the primary and secondary sources, respectively. For a fixed set of {Δx₁, Δy₁, Δx₂, Δy₂}, the linearity of the least-squares procedure in Equation (3) is preserved. ⁴⁵ The minimization of the χ² of the binary model proceeds by first obtaining the shifted PSFs, and then fitting for {f₁, f₂} by solving the linear least-squares procedure as in Equation (3).

To determine whether the binary model is preferable to the single PSF model, we check the difference in χ² between the two, ${\rm{\Delta }}{\chi }^{2}={\chi }_{\mathrm{single}}^{2}-{\chi }_{\mathrm{binary}}^{2}$. As we have an additional five parameters (the flux of the secondary source and the position shifts), we expect Δχ² to also follow a χ² distribution with 5 degrees of freedom. However, we note that there are other reasons (other than a binary object) for the χ² to improve, such as a poorly modeled background feature (such as a star), a miscentering in the PSF of the main source (as the position is fixed, and derived from the orbit fitting in our single PSF photometry), or another transient in the same image (such as an asteroid). We visually inspect all cases where Δχ² ≥ 9, where the probability of the null hypothesis being true (i.e., a single source) is 0.1%. Figure 3 shows two examples for the high-confidence binaries identified in this data.

Zoom In Zoom Out Reset image size

The application of this technique to all ≈25,000 of our detections in the griz bands leads to two objects where several images led to both an improvement in the χ² and S/N ≥ 5 for the fainter source in multiple images, as shown in Figure 4. We also note that Eris's satellite Dysnomia is not resolved in the DES images (Bernstein et al. 2023). While 2014 LQ₂₈ had been previously identified as a binary (Thirouin & Sheppard 2019), 2013 RJ₁₂₄ is a new binary discovery. The objects (612620) 2003 SQ₃₁₇, 2014 QL₄₄₁, and 2016 TT₉₄ also show two images each with Δχ² > 10 and S/N ≥ 5 with no associated transient or artifact, but these images are not enough to confirm their binary status, and require follow-up observations. We note that 2003 SQ₃₁₇ has been imaged as a single point source with HST in Noll (2007), which indicates a probable false positive from our analysis.

Zoom In Zoom Out Reset image size

Our two binaries have multiple resolved observations across several years of the survey, so we can attempt to determine a mutual orbit for these systems. In this particular case, we can relax our requirement for a resolved observation to Δχ² ≥ 8: this is justified, as now we have a strong "prior" that this is a binary object, and the probability of the null hypothesis is still low (0.7%). We can determine the uncertainty in the shifts from the center of mass position by inverting the Hessian of the χ² of the model given by Equation (5).

These position shifts (and their corresponding uncertainties) can be readily transformed into a separation vector r and position angle ϕ ∈ [0°, 180°], avoiding any degeneracy in the determination of the primary source. We fit to each object a Keplerian orbit, where the six orbital parameters $({a}_{m},{e}_{m},{i}_{m},{{\rm{\Omega }}}_{m},{\omega }_{m},{{ \mathcal M }}_{m})$ and the period P_m are used to derive distance and light-travel time corrected sky-plane projections $(\hat{r},\hat{\phi })$. Here, the angles refer to the equatorial plane. Indexing the resolved images by μ, we use a Markov Chain Monte Carlo (MCMC) to sample this seven-dimensional parameter space with the likelihood

$$\begin{eqnarray}&&\begin{array}{c}{ \mathcal L }\propto \\ \quad \displaystyle \prod _{\mu }\exp \left[-\frac{\left({r}_{\mu }-\hat{r}({t}_{\mu }),\phi -{\hat{\phi }}_{\mu }\right){{\rm{\Sigma }}}_{\mu }^{-1}{\left({r}_{\mu }-\hat{r}({t}_{\mu }),\phi -{\hat{\phi }}_{\mu }\right)}^{\top }}{2}\right].\end{array}\end{eqnarray} \tag{ 6 }$$

A uniform prior is also applied to restrict these parameters to physically reasonable values for a bound orbit (a_m , P_m > 0, $0\leqslant {{\rm{\Omega }}}_{m},{\omega }_{m},{{ \mathcal M }}_{m}\lt 2\pi$, 0 ≤ i_m ≤ π, 0 ≤ e_m ≤ 1). We present the results of these fits in Table 1. In particular, we note that only our semimajor axes, eccentricities, and periods are well constrained for each system, as indicated by the large permissible ranges of $({i}_{m},{{\rm{\Omega }}}_{m},{\omega }_{m},{{ \mathcal M }}_{m})$. That is, these large error bars indicate that the orientation of the orbit is poorly constrained. From the orbital elements, we can also derive the masses of these systems (using Kepler's third law), as well as the ratio between the mutual semimajor axis and the system's Hill radius (see, e.g., Parker et al. 2011).

Table 1. Mutual Orbits of the Trans-Neptunian Binaries

Object	am	em	im	Ωm	ωm	${{ \mathcal M }}_{m}$	Pm	Epoch	Mass	Hill Radius
	(106 m)		(deg)	(deg)	(deg)	(deg)	(days)	(MJD)	(1018 kg)	(107 m)
2014 LQ28	${32.9}_{-3.5}^{+7.1}$	${0.64}_{-0.16}^{+0.14}$	${43}_{-2}^{+2}$	${5}_{-3}^{+5}$	${256}_{-6}^{+4}$	${268}_{-20}^{+18}$	${1470}_{-230}^{+390}$	57,615.258	${1.4}_{-0.5}^{+0.7}$	${52}_{-7}^{+7}$
2013 RJ124	${34.9}_{-4.0}^{+6.6}$	${0.14}_{-0.09}^{+0.13}$	${90}_{-15}^{+14}$	${103}_{-40}^{+49}$	${152}_{-43}^{+162}$	${38}_{-22}^{+121}$	${1850}_{-250}^{+290}$	56,545.355	${1.0}_{-0.2}^{+0.3}$	${49}_{-4}^{+5}$

Note. All values correspond to the 68% limits of the posterior distribution marginalized over all other parameters.

Download table as:  ASCIITypeset image

In additional to astrometric data, we can also determine flux ratios. We define the magnitude difference $\delta m\,\equiv | 2.5{\mathrm{log}}_{10}({f}_{1}/{f}_{2})|$, where we take the absolute value to avoid ambiguity in determining which object is the primary in the PSF-fitting procedure. We also have a few back-to-back observations (where the orientation of the system does not change) in multiple bands, and we can immediately determine colors without any ambiguity. These results are included in the data release (see Section 6.4).

2014 LQ₂₈ has several observation pairs, presented in Figure 5. All color pairs are less than 3σ away from the colors being equal, and so the results are consistent with Benecchi et al. (2009), with the colors of the primary being statistically indistinguishable from the secondary, implying similarities in the surface composition of each member of the system.

Zoom In Zoom Out Reset image size

To make physical inferences about some selected subpopulation of TNOs, we will assume that their LCAs are drawn independently from some distribution q(A∣θ_A ), where θ_A is one or more parameter(s) characterizing the distribution. Physically, the observed q(A) will be some convolution of a distribution of intrinsic shapes (or surface variations) with a distribution of obliquity angles of the rotation axes to the line of sight. An excellent overview of the derivation of these distributions is in Showalter et al. (2021). Previous analyses of TNO variability have, however, been severely hampered by selection effects in the published values of LCAs—selections both on which objects were targeted, and which were published. Nondetections of variation or periods often go unpublished, and even published upper limits on LCAs have not been incorporated into population analyses in a rigorous way. An advantage of the DES TNO sample is that we have multiepoch photometry for all of the targets, and we also know that the discovery rate is essentially independent of LCA (Bernardinelli et al. 2022), so there are no LCA-dependent selection effects.

Once a selection on orbital characteristics or H is made, we have the posterior distribution p(A_i ∣{f_ij }) of source i available in the form of samples from its MCMC chain, whether or not this posterior indicates a clear detection of variability. We can now form a posterior likelihood for the distribution of the parameters θ_A of the LCA distribution for this class as

$$\begin{eqnarray}&&p({\theta }_{A}| D)\propto p(D| {\theta }_{A})p({\theta }_{A})\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,p({\theta }_{A})\displaystyle \prod _{i}p(\{{f}_{{ij}}\}| {\theta }_{A})\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\,=\,p({\theta }_{A})\displaystyle \prod _{i}\int {{dA}}_{i}\,p(\{{f}_{{ij}}\}| {A}_{i})q({A}_{i}| {\theta }_{A})\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\propto \,p({\theta }_{A})\displaystyle \prod _{i}\left[\displaystyle \frac{1}{{N}_{\mathrm{samp},i}}\displaystyle \sum _{k=1}^{{N}_{\mathrm{samp},i}}q({A}_{{ik}}| {\theta }_{A})\right].\end{eqnarray} \tag{ 14 }$$

This follows from the use of Bayes' theorem, plus the fact that the samples A_ik from the MCMC chain of object i are drawn from a distribution assuming a flat prior on A_i . We can thus use Equation (14) to assign a likelihood to any postulated LCA distribution for a population by summing over all of the MCMC outputs, regardless of whether they indicate a decisive detection of variability.

In a similar vein, we can derive an estimate of the posterior distribution of color for any individual object by applying a chosen prior q(A_i ∣θ_A ) to the MCMC outputs, e.g.,

$$\begin{eqnarray}&&\begin{array}{c}p(b-{b}^{{\prime} }| \{{f}_{{ik}}\},{\theta }_{A})\\ \propto \displaystyle \sum _{k=1}^{{N}_{\mathrm{samp},i}}q({A}_{{ik}}| {\theta }_{A})\delta (b-{b}^{{\prime} }+2.5{\mathrm{log}}_{10}{\bar{f}}_{{ibk}}/{\bar{f}}_{{{ib}}^{^{\prime} }k}).\end{array}\end{eqnarray} \tag{ 15 }$$

One could in principle attempt to constrain the physical parameters of the TNO subpopulations by propagating such parameters through to a functional form for q. Such modeling involves specifying the geometry of the objects, their surface scattering properties, and distribution of rotation axes distributions, as described in detail by Showalter et al. (2021). This is before one even considers the nature of albedo variations across the object's surface. The available information is far short of what would be necessary to constrain the large number of free parameters in any realistic physical model.

We opt instead to adopt a simple heuristic parametric form for q(A), which we will fit to various subpopulations with the intention of looking at the distinctions across the trans-Neptunian region, rather than attempting to extract physical parameters. In other words, we will use q(A) as a kind of genetic marker for the relationships among the subpopulations.

A flexible family of distributions normalized over the range 0 ≤ A < 1 is the β distribution, whose probability distribution function depends on two shape parameters $\alpha ^{\prime}$ and $\beta ^{\prime}$ and is usually written as

$$\begin{eqnarray}&&q{(A| \alpha ^{\prime} ,\beta ^{\prime} )\propto {A}^{\alpha ^{\prime} -1}(1-A)}^{\beta ^{\prime} -1}.\end{eqnarray} \tag{ 16 }$$

We transform the parameters slightly, using first the mean of the distribution, $\bar{A}=\alpha ^{\prime} /(\alpha ^{\prime} +\beta ^{\prime} ),$ and a second parameter $s\equiv {\mathrm{log}}_{10}(\alpha ^{\prime} +\beta ^{\prime} ),$ a "sharpness" parameter: at fixed $\bar{A},$ the β distribution becomes narrower about the mean as s increases. Thus, we will have

$$\begin{eqnarray}&&q{(A| {\theta }_{A})=q(A| \bar{A},s)\propto {A}^{\bar{A}{10}^{s}-1}(1-A)}^{(1-\bar{A}){10}^{s}-1}.\end{eqnarray} \tag{ 17 }$$

The β distributions are a superset of the power-law distributions considered in early modeling of the variability of TNOs (Lacerda & Luu 2006). Since Bernardinelli et al. (2022) showed that the selection function of the DES survey is nearly independent of A at fixed mean apparent magnitude, we do not need to introduce selection terms into our posterior probability for the parameters of q(A).

For a chosen sample of DES TNOs, we calculate the posterior probability of $(\bar{A},s)$ using Bayes' theorem, $p(\bar{A},s| D)\propto p(D| \bar{A},s)p(\bar{A},s),$ where the data D consist of the MCMC chains sampling the A_i distribution for TNO i within the sample. This can be evaluated across the full $(\bar{A},s)$ space numerically using Equation (14), and assuming uniform priors over $\bar{A}$ and s. We will also marginalize over s to obtain $p(\bar{A}| D)$ for different subpopulations. The overlap between these distributions for two different subpopulations becomes a measure of their morphological similarity.

Table 2 defines several disjoint samples for which we calculate $p(\bar{A},s)$, and Figure 7 plots the posterior probability in the $(\bar{A},s)$ space, and marginalized down to $p(\bar{A})$ for each. It is expected that more massive TNOs will be more spherical due to self-gravity—and hence have q(A) distributions weighted toward lower values—so it is important to control for size when testing differences. All but one of the TNO samples we define therefore are restricted to 6 < H_r < 8.2 (and have 〈H_r 〉 = 7.4 ± 0.2), a range in which most of the dynamical families have a useful fraction of their DES-detected TNOs. The "Big" sample contains all DES TNOs with H_r < 6. Aside from Eris, this sample's members have 3.3 < H_r < 6.0, and only two of them meet the cold-classical criterion, so the "Big" subset is essentially composed of dynamically excited TNOs.

Zoom In Zoom Out Reset image size

Table 2. TNO Sample Definitions

Subset	N	Size	Median (Hr )	Dynamical Class	Inclination	$\bar{A}$ 68% CL
Cold classical (CC)	95	6 < Hr < 8.2	7.18	Classical	ifree < 5°	0.149–0.165
Hot classical (HC)	261	6 < Hr < 8.2	7.38	Classical	ifree > 5°	0.124–0.128
Resonant	170	6 < Hr < 8.2	7.56	Resonant	⋯	0.127–0.136
Detached	133	6 < Hr < 8.2	7.39	Detached	⋯	0.114–0.115
Scattering	34	6 < Hr < 8.2	7.57	Scattering	⋯	0.148–0.176
Big	36	Hr < 6	5.43	⋯	⋯	0.081–0.095

Note. Each row gives the conditions defining one subsample of the DES TNO sample, and the resultant number and median absolute magnitude H_r of the sample. Dynamical classifications are taken from Bernardinelli et al. (2022), and the free inclinations i_free of all classical objects are provided by Huang et al. (2022). The final column gives the 68% confidence interval for the mean LCA $\bar{A}$ of the distribution of the subset's members' A values. The dynamical subpopulations follow the scheme of Gladman et al. (2008)—we refer the reader to Bernardinelli et al. (2022) for a full description of the classification procedure for these objects.

Download table as:  ASCIITypeset image

Figure 7 immediately indicates that the "Big" H_r < 6 sample is indeed less variable (lower $\bar{A}$) than the 6 < H_r < 8.2 members of any of the dynamical classes, as might be expected from the effects of self-gravity. We wish to determine whether two subsamples D_A and D_B could be drawn from the same parent distribution of parameters. Hypothesis ${{ \mathcal H }}_{1}$ is that they are from a single β distribution. Hypothesis ${{ \mathcal H }}_{2}$ is that they are from distinct distributions with two pairs of β-distribution parameters. We give two ways of determining the strength of a conclusion in favor of ${{ \mathcal H }}_{2},$ distinct distributions. A frequentist-oriented statistic is

$$\begin{eqnarray}&&{\rm{\Delta }}{\chi }^{2}\equiv 2\mathrm{log}\displaystyle \frac{\left[{\max }_{\bar{A},s}p(\bar{A},s| {D}_{A})\right]\left[{\max }_{\bar{A},s}p(\bar{A},s| {D}_{B})\right]}{{\max }_{\bar{A},s}p(\bar{A},s| {D}_{A},{D}_{B})}.\end{eqnarray} \tag{ 18 }$$

This statistic compares the maximum posterior likelihood of ${{ \mathcal H }}_{2}$ to ${{ \mathcal H }}_{1}.$ Under the assumption that the parameters are Gaussian distributed, this statistic would have a χ² distribution with ν = 2 degrees of freedom, if the two populations were indeed drawn from the same model. Then the probability of Δχ² > 4.6 would be 10% and Δχ² > 9.2 would be 1%. We will designate these as "indicative" and "decisive" evidence of differences in the underlying q(A) distribution of two samples. In practice, we do not meet the assumption of Gaussianity, so these probability ranges should be seen as an approximation. By this criterion, the "Big" sample is indeed decisively distinct from all of the higher-H_r subpopulations, except for being indicative of the detached TNOs.

A more Bayesian statistic is the (log) evidence ratio, defined as

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal R }=\mathrm{log}\displaystyle \frac{p(D| {{ \mathcal H }}_{2})}{p(D| {{ \mathcal H }}_{1})}\\ \quad =\,\mathrm{log}\displaystyle \frac{\int d\bar{A}\,{ds}\,p({D}_{A}| \bar{A},s)p(\bar{A},s)\times \int d\bar{A}\,{ds}\,p({D}_{B}| \bar{A},s)p(\bar{A},s)}{\int d\bar{A}\,{ds}\,p({D}_{B},{D}_{A}| \bar{A},s)p(\bar{A},s)}.\end{array}\end{eqnarray} \tag{ 19 }$$

Adopting the scale of Jeffreys (1961), values of ${ \mathcal R }\gt 2.30,$ 3.45, and 4.61 are labeled as "strong," "very strong," and "decisive" in favor of ${{ \mathcal H }}_{2}.$ The first and last of these correspond to ≈10% and ≈1% chances of H₁ being correct, if we assign equal prior probability to ${{ \mathcal H }}_{1}$ and ${{ \mathcal H }}_{2}.$

For our problem, we require a normalized prior $p(\bar{A},s).$ We choose a prior that is uniform over the range $0.06\lt \bar{A}\lt 0.19,$ 0.7 < s < 1.9. The evidence ratio is indeed decisively in favor of the Big size distribution being distinct from all other subsets, except the detached population, for which ${ \mathcal R }=0.54$ is indecisive about ${{ \mathcal H }}_{1}$ vs ${{ \mathcal H }}_{2}.$

Also conclusive from the Δχ² test and "very strong" at ${ \mathcal R }=4.33$ in the evidence ratio is the hypothesis that the CCs at i_free < 5° are distinct, namely more variable, than the HCs (i_free > 5°). The CCs are believed to be the most dynamically pristine, and are known to have redder colors and a higher rate of binarity (Stephens & Noll 2006; Fraser & Brown 2012; Nesvorný & Vokrouhlický 2019). Our results demonstrate that variability is another physical distinction between higher- and lower-inclination classical TNOs.

Before proceeding further, we check whether the β distributions are adequate descriptions of the measured size distributions. In the left panel of Figure 8, we plot in solid lines the q(A) β-distributions that maximize the posterior probability for each of the Big, HC, and CC samples. Both the differential and cumulative distributions are given. The dotted curves plot the posterior distributions for the observed data, i.e., we average the posterior distributions p(A∣D_i ) = p(D_i ∣A)q(A)/∫dA p(D_i ∣A)q(A) over all members of the labeled subset. ⁴⁷ The agreement is good, with deviations that are larger for the subsets with fewer members.

Zoom In Zoom Out Reset image size

Next we examine the size-controlled samples of the other dynamically excited populations, the resonant, scattering, and detached TNOs. Neither the resonant nor scattering TNO populations can be confidently distinguished from the HC's or CC's in terms of their q(A) parameters. Indeed the evidence ratio ${ \mathcal R }=-3.53$ for the HC-resonant pair indicates very strong preference for the hypothesis of a shared distribution. By the Δχ² criterion, there is indicative evidence that the resonants are distinct from the CCs, and that the scattering TNOs are distinct from the HCs. We can summarize this by saying that the variabilities of the resonant and scattering populations are fully consistent with being somewhere between the HCs and CCs.

The last dynamically excited population, the detached TNOs, lies between the HC and Big loci, and is decisively distinct from CCs in variability by all criteria. The evidence ratio indicates that the detached TNOs could be from the same distribution as either the Big or HC, while the Δχ² offers indicative and decisive distinction from the Big and HC populations, respectively. A surprising result is that the detached and scattering TNOs are decisively inconsistent by both criteria.

The results of both forms of model comparison test can be easily summarized: if the 95% confidence level (CL) ellipses of two populations do not overlap in the upper panel of Figure 7, then they are decisively distinct from each other by one or both of the more formal criteria.

Further division of the DES TNOs into smaller subsets than listed in Table 2 did not yield any ability to distinguish statistically significant differences. In particular, we have tested for dependence on H_r within the HC sample by splitting at the median H_r ; neither hypothesis test indicates a preference for distinct q(A) between the halves. Similarly we split the resonant class into high-a and low-a halves; high-i and low-i halves; and Plutinos vs. others. In no case did the splits yield significantly better modeling. Of course this does not preclude the existence of such dependence, but the quantity and quality of our data are insufficient to detect any.

Comparison of these conclusions with previous investigations, and implications for the origin of the populations, are discussed in Section 6.

This paper describes a series of techniques to obtain optimal measures of fluxes, colors, binarity, and variability of moving sources like TNOs, with principled estimation of uncertainties. We apply these techniques to the >800 TNOs detected by the DES, yielding accurate colors for each, discovery of several likely new wide binaries, and the ability to make inferences about the relative levels of variability of TNO populations with more rigor and confidence than previous studies of TNO variability. Aside from the ≈10 × larger sample of TNOs than preceding studies, this analysis addresses several methodological issues that affect some investigations:

• 1.  
  The photometry methods obtain the optimal signal-to-noise level on TNO fluxes from every frame by using PSF-fitting techniques, and also subtract any background sources that may be behind the TNO—even those that may be undetectable at the noise level of the target image.
• 2.  
  The estimates of mean H_r and colors make use of all of the exposures of a given TNO, avoiding biases that would result from using only the images with S/N sufficient for detection. The color estimates and their uncertainties include the effects of potential variability of the source while exploiting any information from near-simultaneous colors.
• 3.  
  Estimates of the variability yield a full probability distribution for A, rather than arbitrary decisions about and treatment of nondetections. This is true whether or not a period is detectable. This enables a correct quantitative comparison to candidate distributions q(A), even for TNOs with poorly constrained A.
• 4.  
  The analysis includes all of the TNOs detected by DES that meet selection criteria on dynamics or H. In other words the selection function for the input sample is well defined and independent of variability. Most previous studies have relied on samples drawn from the very heterogeneous literature, which are subject to unquantifiable selection biases. Most pernicious might be a publication bias, whereby many authors might only publish amplitudes of objects if they have detected variation, leaving out nondetections.
• 5.  
  The 3–5 yr temporal spread of the DES photometry is far longer than any TNO rotation period, insuring that light-curve phases are randomly sampled and our distance-corrected variability estimates are independent of the period.
• 6.  
  The DES sample is large enough to make meaningful comparisons between dynamical groups at a fixed and limited range of H so we can separate dynamical-state dependence from size dependence.

These techniques should be of great value for future, even larger-scale, sky surveys such as the Legacy Survey of Space and Time, which will measure at least an order-of-magnitude more TNOs than DES has found. The vast majority of TNO measurements in any survey have S/N levels too low to allow for traditional measurements of light curves, but there is still a great deal of information available from their colors and variability.

The level of variability in our "Big" sample (3.3 < H_r < 6.0, plus Eris) is decisively lower than that of any of the dynamical classes' 6.0 < H_r < 8.2 members. A relative paucity of variability at Δm > 0.3 mag among H ≲ 6 objects has been noted in past works, and is readily ascribed to the effects of greater surface gravity in forcing objects toward sphericity. Benecchi & Sheppard (2013), for example, found a correlation of LCA versus H_V at 3σ significance in a sample comprising 32 objects at H_V < 6 that they measure, combined with 96 objects at 0 < H_V < 12 with Δm measures drawn from the extant literature, mixing together all dynamical classes. As noted by Showalter et al. (2021) in reviewing variability statistics, we should be alert to selection biases and publication biases in samples drawn from the literature. They opt to nonetheless "proceed on the assumption that our sample is unbiased simply because we have no clear alternative." With the DES sample, we no longer need to make this leap. The measurements by Alexandersen et al. (2019) also have a selection criterion that should be independent of A, since they target all 63 OSSOS TNO detections at 6 < H_r < 11 within a few chosen fields of view of the HyperSuprimeCam imager. Within this sample, they report a correlation between σ_mag (the observed magnitude dispersion per object) with H_r at p = 0.013 false-positive rate in a Spearman rank correlation test. The correlation persists at p = 0.015 when the sample is restricted to dynamically excited populations (i.e., excluding CCs). ⁴⁸ Our data can put to rest any lingering doubt about the reality of a variability-H connection, since we use samples guaranteed to be free of variability biases, and can show that the 3.3 < H_r < 6 sample is distinctly less variable from the smaller objects in every dynamical class (with the possible exception of detached TNOs). Our distributions for H_r ≤ 6 are also in general agreement with the targeted sample of primarily large objects by Kecskeméthy et al. (2023), who showed that ≈55% of their sample has Δm ≤ 0.2 mag (compare to the upper-left panel of Figure 8).

Another decisive result of our variability study is a higher level of variability in "cold" classical TNOs (i_free < 5°) than in "hot" classicals (i_free > 5°). Correlations of variability with inclination, or higher variability in CCs than in non-CCs, have been noted in other samples. Benecchi & Sheppard (2013) reported differences at the ≈95% CL between the variabilities of classical TNOs and those of resonant and/or scattering objects, and detected a correlation between Δm and i at similar confidence across the full TNO sample. The decisive difference between HC and CC populations' variabilities seen in our data reinforce the suspicions based on dynamical and color differences that the CCs have a distinct physical origin, and we should therefore split these in any comparison to other dynamical families. We should also realize that the simple i_free cut is probably not the definitive division between these two physical populations, so our HC and CC samples are a bit mixed, which means we are potentially underestimating their distinction.

Thirouin & Sheppard (2019) compared their own measures of CC variability against other populations' data in the literature, reporting that 36% of CCs have Δm < 0.2 mag—quite close to our value on the blue curve in the upper left of Figure 8—compared to 65% of a mixed bag of "other" TNOs. Our data agree with this trend in the sense of CCs being the most variable of all of the dynamical classes, at fixed H range. Similarly, Strauss et al. (2023) measured a sample of 24 objects (that they suggested are primarily cold Classicals) with high cadence observations during 4 hr, finding that roughly half their sample had Δm ≤ 0.4 magnitudes, similar to our CC distribution. Alexandersen et al. (2019) also reported a correlation between i and H_r within the CC class at 94% CL. We do not confirm this trend (nor falsify it): splitting the CC sample at the median i_free does not yield statistically significant distinctions in q(A)–neither does splitting the much larger HC sample at its median i_free of 15°. Our data support only the dichotomy of variability between HC and CC dynamical types, and no further evidence exists for a continuous gradient in i_free.

The resonant and scattering TNOs are not distinguishable from either the HC or CC TNOs with regards to their variability levels. This is consistent with any scenario in which resonant or scattering TNOs share a physical origin with any mixture of the classicals' source populations. The surprising result, however, is that the detached and scattering populations have decisively distinct variability distributions (${\rm{\Delta }}{\chi }^{2}=18.2,{ \mathcal R }=6.5$). The detached and resonant populations meet the frequentist criterion for decisive distinction (${\rm{\Delta }}{\chi }^{2}=12.0,{ \mathcal R }=2.4$). This would exclude scenarios in which members of the scattering population have their perihelia raised and become detached by a process that is independent of the physical nature of the TNOs. This might be consistent with the observations if TNOs were substantially physically altered in the detaching process, to lower their mean variability by a factor ≈1.5. An isotropization of spin axes that were initially aligned normal to the ecliptic could lower the mean apparent variability by this much, but one would need to explain why the detached TNOs are isotropized but the scattering TNOs are not. A weaker version of this restriction applies to the detached TNOs arising through the resonant channel. This is a puzzle for models of the origin of the detached TNOs.

We discover one new binary, in addition to the previously known 2014 LQ₂₈, and characterize their colors and sizes. Our ability to identify these binaries is a proof-of-concept of the combination of SMP with binary deblending. These objects also have enough astrometric data over the sparse DES observations that we can determine a mutual orbit for the binary system. All of our objects are tightly bound, with a_m /R_H ≈ 5%. Their masses are in significant agreement with masses of other objects of similar sizes reported in the literature (e.g., Grundy et al. 2011; Parker et al. 2011).

2014 LQ₂₈ and 2013 RJ₁₂₄ are cold Classicals, adding to the now large sample of characterized binary CCs (Noll et al. 2020). Notably, 2014 LQ₂₈ is the second-largest CC in the DES sample (H_r ≈ 5.7), which is in line with the hypothesis that most large CCs are binaries (Thirouin & Sheppard 2019). We can estimate from the differences in H_r and masses that 2013 RJ₁₂₄ is either half as dense or half as reflective as 2014 LQ₂₈. Assuming the same albedo for both objects, the volume ratio between both objects is 2.8, while the mass ratio is ≈1.38. On the other hand, if we assume that 2014 LQ₂₈'s albedo is twice that of 2013 RJ₁₂₄, we obtain roughly the same volume for both objects, and such albedo differences are in line with measurements from thermal data (Vilenius et al. 2012).

Readers should resist the temptation to use our binarity detections to measure the fraction of wide binaries in the TNO populations. While the binary search was exhaustive over all of our images, it was not systematically characterized, meaning that we have not characterized our ability to discover binaries from our images as a function of mutual orbit parameters, seeing, and S/N. The binary discovery fraction in the DES data strongly underestimates the binary fraction in TNO populations, since the majority of the objects detected in our survey are near the detection threshold, where our ability to resolve binaries is inherently diminished. Further observations of each of these binary objects, with even slightly better angular resolution, would greatly improve constraints on their orbits and masses.

The data release presented here substantially increases the number of TNOs with known colors and LCAs in the literature, and is the largest of such a data release from a single survey with a consistent filter set and calibration. Digital resources containing the results of application of these techniques to the >800 TNOs detected in the DES Wide Survey are available on Zenodo at doi:10.5281/zenodo.8231238. Included in the data release are the following:

• 1.  
  A table of all of the valid individual photometric flux measurements of all of the TNOs.
• 2.  
  A table of all photometric measurements of the resolved binaries.
• 3.  
  For the two TNO binaries, we also include the MCMC chains with their derived mutual orbits.
• 4.  
  A Jupyter notebook containing code that can run MCMC chains to sample the space of mean fluxes and variability, given a set of multiband, sparsely sampled fluxes as described in Section 4.
• 5.  
  For each TNO, a table giving the output MCMC chains produced from its photometry.
• 6.  
  An update of the table of TNO properties given by Bernardinelli et al. (2022), augmented with the means and standard deviations of H_r , colors, and LCAs A for each object, and derived from the MCMC chains.
• 7.  
  A Jupyter notebook that derives the plots and statistics of this paper from the MCMC chains, and demonstrates how a user could apply the methodology to other subsets of the DES TNOs or to other data.

We also include the scene modeling photometry software for measuring fluxes of both single and binary sources, as well as a Jupyter notebook with some examples of its usability, in a different Zenodo repository at doi:10.5281/zenodo.8231232 (Bernardinelli et al. 2023).