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

We report the methods of and initial scientific inferences from the extraction of precision photometric information for the >800 trans-Neptunian objects (TNOs) discovered in the images of the Dark Energy Survey (DES). Scene-modeling photometry is used to obtain shot-noise-limited flux measures for each exposure of each TNO, with background sources subtracted. Comparison of double-source fits to the pixel data with single-source fits are used to identify and characterize two binary TNO systems. A Markov Chain Monte Carlo method samples the joint likelihood of the intrinsic colors of each source as well as the amplitude of its flux variation, given the time series of multiband flux measurements and their uncertainties. A catalog of these colors and light-curve amplitudes A is included with this publication. We show how to assign a likelihood to the distribution q(A) of light-curve amplitudes in any subpopulation. Using this method, we find decisive evidence (i.e., evidence ratio <0.01) that cold classical (CC) TNOs with absolute magnitude 6 < H r < 8.2 are more variable than the hot classical (HC) population of the same H r , reinforcing theories that the former form in situ and the latter arise from a different physical population. Resonant and scattering TNOs in this H r range have variability consistent with either the HCs or CCs. DES TNOs with H r < 6 are seen to be decisively less variable than higher-H r members of any dynamical group, as expected. More surprising is that detached TNOs are decisively less variable than scattering TNOs, which requires them to have distinct source regions or some subsequent differential processing.


INTRODUCTION
The trans-Neptunian region is a distant reservoir of small bodies that trace the formation history of the Solar System (Nesvorny 2018).We currently know of more than 3000 of these objects, with recent surveys capable of discovering several hundreds at a time (e.g.Petit et al. 2011;Bannister et al. 2018;Bernardinelli et al. 2022).The combination of dynamical and physical characterizations of these populations has led to our understanding of several key aspects of the formation of the outer Solar System (see Gladman & Volk 2021, for a recent review).Photometric measurements of these trans-Neptunian objects (TNOs) are of particular interest, and analyses of such data have led to the understanding of the bulk properties of surface shapes (Showalter et al. 2021), the determination of distinct compositional classes from surface colors (e.g.Doressoundiram et al. 2008;Fraser & Brown 2012;Schwamb et al. 2019), characterization of the size distribution (e.g.Bernstein et al. 2004;Fraser et al. 2014;Kavelaars et al. 2021), discovery of a large fraction of binary systems (Stephens & Noll 2006;Parker et al. 2011;Noll et al. 2020), and the identification of a collisional family (Brown et al. 2007).
The Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005Collaboration , 2016) ) received an allocation of 575 nights on the 4m Blanco telescope at Cerro Tololo, using the Dark Energy Camera (DECam, Flaugher et al. 2015) to cover 5000 deg 2 of the sky in the grizY photometric system from 2013 to 2019.The survey's primary objective has been to study the distribution of dark matter and the nature of dark energy (The Dark Energy Survey Collaboration 2022), but the data have enabled the discovery of hundreds of outer Solar System objects (Bernardinelli et al. 2020(Bernardinelli et al. , 2022) ) -we refer the reader to these two publications for a comprehensive presentation of the discovery pipeline.All 814 DES objects have been dynamically classified following Khain et al. (2020), and the DES survey simulator (Bernardinelli et al. 2022) allows us to carefully estimate our detection biases as a function of dynamical population, magnitude, color and light curve amplitude.
The DES observing strategy was optimized for extragalactic science, which means that it is less efficient in terms of TNO discoveries per night of telescope time than surveys designed for TNO discovery (Bannister et al. 2018;Trilling et al. 2023).DES observes a given region too many times, in too many filters, over too long a time span, to be optimal for discovery.This redundancy has the advantage that we can extract significantly more information about each discovered source than an optimized discovery search.Each object has been observed many times: between 6-10 times in each of the grizY bands, on average, over the six years of data collection.This means that each object is imaged on average between 30 and 50 times, depending on whether or not the object was inside the footprint for the entirety of the survey and each region's cadence.This allows estimation of colors in the griz bands (Y is generally too low in signal-to-noise ratio, SNR) and also an estimate of the time variability of each source.Inferring colors and variability simultaneously allows the color estimates to include uncertainties that arise when colors are measured from non-simultaneous exposures in the presence of variability, while also exploiting the occasions when DES targets a given source in multiple filters in quick succession.This paper will present the methods and results of such estimates for the full DES TNO catalog, as well as some initial physical inferences that can be made from the variability information.Section 2 will describe the extraction of fluxes, and Section 4 the extraction of colors and light curve amplitudes (LCA's) from the flux time series.In Section 5 we examine the distributions of light-curve amplitudes in different TNO sub-populations as an indicator of different physical states.Scientific analyses of the color catalog will appear in future publications.
The process of extracting optimal flux estimates from the pixel data can also be generalized to fit two fluxes of a potential binary pair to each exposure of a given TNO.In Section 3 we describe this process and present the resulting binary candidates found among the DES TNOs.
Accompanying this publication is a data release with the ≈ 30, 000 photometric measurements, absolute magnitudes, colors and light curve amplitudes of our objects, as well as some of the software required for the analysis we present here.The data release is available at [forthcoming], and discussed in Section 6.4.

INDIVIDUAL PHOTOMETRIC MEASUREMENTS
To extract unbiased flux measurements of each TNO from the DES images, we use the images in which the object is detected and also images the orbit predicts should contain the object, but there is no detection in the catalog ("non-detections")-similar to the sub-threshold significance (STS) measurement of Bernardinelli et al. (2020).We will collectively call these "observations."We first start by describing the photometry for each individual observation of a TNO.
We determine the TNO positions in each exposure by using the values predicted by the orbit fit, and model each detection using the scene modeling photometry (SMP) technique commonly used in supernovae type Ia cosmology (Brout et al. 2019).We use the DES point spread function (PSF) model of Jarvis et al. (2020), with PSFs derived for the full field of view of each DES exposure.A full description of the DES calibration procedure is presented in Burke et al. (2017) and Sevilla-Noarbe et al. (2020), and comparisons between DES and Gaia show an exquisite calibration with 3 mmag root-mean-square differences between the two surveys (Abbott et al. 2021).
Each location where a TNO is observed on some single night has n − 1 images in the same band at different epochs in which the TNO is not present.We posit that the background is composed of point sources on a m × m square grid, with spacing of roughly the 1σ width of the PSF, so the PSF will blur them into a smooth distribution to represent any extended sources.The grid sources have fluxes P uv , and are centered on the TNO location.Given the ∼ 0.95 FWHM of typical DES imaging, we place sources every 0.35 and define m = 20, so these sources span a 7 × 7 region.These point sources are then mapped into the (u, v) pixel coordinates of each image by inverting the DECam astrometric model (Bernstein et al. 2017).The astrometric solution and the PSF models are both functions of source color; we begin by assuming a nominal g − i = 0.61 color (typical of stellar sources) for each background source.We will let the u, v symbols serve both as indices into the grid of background sources, as well as their exact positions in the pixel coordinate system.The expected signal in pixel (i, j) for exposure µ for this mosaic of sources is modelled as the convolution of each P uv with exposure's PSF derived for its location (u, v); plus some constant background level b µ : In the single exposure ν where the TNO is present, we adopt the same model with an additional point-source term: That is, f TNO represents the integrated flux of the TNO at this epoch and band.Due to the short exposures times (90 seconds for most images), these sources are not trailed, and so corrections such as pill apertures (e.g.Fraser et al. 2016) are not needed.Initially we use a default color for the TNO in evaluating the PSF and its expected pixel position.This model, then, has N = n + m 2 + 1 free parameters, and nk 2 − N degrees of freedom, where k > m is the number of pixels in the postage stamp modeled (we used stamps with 30 × 30 pixels, corresponding to a sky area of 7.8 ×7.8 ).We fit these using a least squares minimization, comparing the model to the measured pixel values Im µ : . (3) The (constant) k µ,ν terms correct for the different zero-points in each exposure, bringing them to a common flux scale (Burke et al. 2017), and σ 2 µ,ij is the noise variance at each pixel from sky background and detector noise.This model is a linear system determined by a design matrix A that contains all PSF realizations and the constant, unitless background terms; parameters X = {P uv , b µ , f TNO }; and the target matrix Y (i.e., Y = AX) that represents all images.Thus, we can solve for the parameters using a standard linear least-squares solution.
The initial estimates of the σ ν,ij account only for shot noise from the uniform sky background, and for detector noise.To account for the additional variance due to shot noise from the sources, we update the weights in exposure ν using the first least-squares solution for M ν ij : where g ν,ij is the amplifier-dependent gain in each pixel and exposure (Bernstein et al. 2018).We refit the model with these new variances, and derive the flux uncertainties from the photometric solution.Once all of the exposures of a given TNO have been measured, we estimate its mean g − i color,1 then repeat the entire measurement process for each TNO apparition while using this g − i value in the color-dependent astrometric solution and PSF for the TNO flux.All ≈30,000 individual observations were visually inspected, and cases where the scene-modelling procedure failed were discarded.We show successful examples in Figure 1 and failures in Figure 2. We note that this methodology can be easily applied to extended sources (for example, comets) and to larger areas of the sky: in Bernardinelli et al. (2021) we applied this methodology to 400 × 400 pixel stamps in order to detect C/2014 UN 271 's extended coma.A generalization to binary sources is presented in the next section.The modification for a binary TNO requires changing Equation 2 to include a second point source term: where the indices 1 and 2 correspond to the primary and secondary sources, respectively.For a fixed set of {∆x 1 , ∆y 1 , ∆x 2 , ∆y 2 }, the linearity of the least squares procedure in Equation 3 is preserved2 .The minimization of the χ 2 of the binary model proceeds by first obtaining the shifted PSFs, and then fitting for {f 1 , f 2 } 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 χ 2 between the two, ∆χ 2 = χ 2 single − χ 2 binary .As we have an additional 5 parameters (the flux of the secondary source and the position shifts), we expect ∆χ 2 to also follow a χ 2 distribution with 5 degrees of freedom.However, we note that there are other reasons (besides a binary object) for the χ 2 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 The first row shows a case where the TNO was near an extended background galaxy.This example is remarkable because, even though there is no usable photometry for the TNO, the model reproduces the spiral structure of the galaxy really well, leading to a near perfect subtraction.The second row shows a similar case of a bright background source near the TNO, but in this case the model failed to reproduce this background.The third row corresponds to a bright source near the edge of the stamp, with its scattered light precluding any photometric measurement nearby.The final row is a case where the least squares procedure of Equation 3 led to an incorrect solution (namely, the model predicts a large negative flux for the central point source).photometry), or another transient in the same image (such as an asteroid).We visually inspect all cases where ∆χ 2 ≥ 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.
The application of this technique to all our ≈ 25, 000 detections in the griz bands leads to 2 objects where several images led to both an improvement in the χ 2 , 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 28 had been previously identified as a binary (Thirouin & Sheppard 2019)  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 317 has been imaged as a single point source with HST in Noll (2007), which indicates a probable false positive from our analysis.

Mutual orbit determination and derived properties
Our two binaries have multiple resolved observations across several years of the survey, 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 ∆χ 2 ≥ 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 χ 2 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 6 orbital parameters The other cases of objects where the binary fit yields a statistically significant improvement in the χ 2 arise when there is a poorly subtracted background object in the single-source fit, or when the binary fitting corrects small errors in the assumed pixel position of the primary object caused by atmospheric turbulence.
(a m , e m , i m , Ω m , ω m , M m ) and the period P m are used to derive distance and light-travel time corrected sky-plane projections (r, φ).Here, the angles refer to the equatorial plane.Indexing the resolved images by µ, we use a Markov Chain Monte Carlo to sample this seven dimensional parameter space with the likelihood A uniform prior is also applied to restrict these parameters to physically reasonable values for a bound orbit (a m , P m > 0, 0 1).We present the results of these fits in Table 1.In particular, we note that only our semi-major axes, eccentricities and periods are well constrained for each system, as indicated by the large permissible ranges of (i m , Ω m , ω m , 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 Note-All values correspond to the 68% limits of the posterior distribution marginalized over all other parameters.
well as the ratio between the mutual semi-major axis and the system's Hill radius (see, e.g., Parker et al. 2011).
In additional to astrometric data, we can also determine flux ratios.We define the magnitude difference δm ≡ |2.5 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 28 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.

EXTRACTING COLORS AND LIGHT CURVE AMPLITUDES
After the scene-modelling photometry, each TNO i ∈ {1, 2, . . ., N TNO } has a series of measured fluxes f ij for j ∈ {1, . . ., N obs,i } observations in bands b ij , with nearly Gaussian uncertainties σ ij .We will assume here that the fluxes have been adjusted for their heliocentric and geocentric distances to represent fluxes that would be observed with both distances at d ref = 30 AU.What are the best estimates of the mean fluxes fib for source i in band b? From these, the best estimated absolute magnitudes H ib and colors (b − b ) i ≡ H ib − H ib can be determined, and their uncertainties.
A key element of the uncertainties is the potential for variability in the TNO fluxes due to rotation coupled with asphericity and/or surface inhomogeneities of the sources.The level of variability is a quantity of interest itself as an indicator of the physical state of the TNO, as well as being a source of noise in flux/color measures, so we would like to have a principled estimate (and uncertainty) of the variability amplitude, as well as estimates of mean flux/color that have been marginalized over the variability.
The traditional means to determining colors would be to obtain high-S/N multiband observations of each source within a time interval that is short compared to the variability period, allowing the measure of "instantaneous" colors.Alternatively one could take enough observations in each band to be able to determine the period, reconstruct the light curve as a function of variability phase φ, and get phase-averaged fluxes fib for each band b.Unfortunately the DES observing cadence does not admit either method.In each band, the source is typically observed 8 times over a five-year span, far too sparse to determine a period, never mind construct a light curve.Furthermore only occasionally are two observations of a given TNO are made on the same night, making instantaneous colors usually unavailable or low-S/N.Consecutive observations in distinct bands do occur, however, so we would like a method that can exploit these events for deriving accurate colors.
We have developed a method to make optimal use of the information that we do have in determining the mean fluxes and the level of variability.Our approach is conceptually similar to that used by Schemel & Brown (2021) to measure photometric variability of Jovian Trojans.We assume a model in which the flux at observation j is assumed to equal Here we have dropped the index i of the TNO, for brevity, since the process is independent for each TNO.A phase function h(φ) is a function of the variability phase φ ∈ [0, 1) with h = 0 when averaged over phase, and max φ h−min φ h = 2.The parameter A then gives the semi-amplitude of the fractional variation of true flux, and fb is the time-averaged mean flux in band b.Note that we define the mean and the light-curve semi-amplitude (LCA) in flux space rather than in magnitude space.
The peak-to-peak magnitude variation would be ∆m ≡ 2.5 log 10 (1 + A)/(1 − A).We assume the convention A ≥ 0, and only A < 1 is physically possible.We make the following critical assumptions: • All bands have the same phase function h and amplitude A, i.e. the variation is achromatic.
• The variability phase φ j is a random unit deviate for each observation, i.e. we know nothing about the period except that it is short enough to leave us with no knowledge of the relative phases of different observations-unless observations j and k are within 1 hour of each other, in which case we assume φ j = φ k but the value is unknown.
• We will assume that the phase function is h(φ) = sin(2πφ).This is likely inaccurate for nonellipsoidal geometries, but our goal is to obtain A values that are indicative of variability level, so as to admit comparison of populations' variabilities, so we don't care that individual values of A are precise LCA's.
With the first assumption in hand we can write the probability of the observations f j with Gaussian uncertainties σ j as To obtain p {f j }|{ fb }, A we need to marginalize Equation ( 8) over all possible phases using the second assumption.We divide all the observations j into sets indexed by s such that j and j are in the same set s if and only if they occur within one hour of each other.With this convention, the marginalization over light-curve phases becomes Numerically, the integral over φ can be executed as a sum over ≈ 20 equally-spaced samples of 0 ≤ φ < 1.Using Bayes' theorem and assuming uniform priors for fb and for 0 ≤ A < 1, we can take We created a straightforward Metropolis-Hastings Markov Chain (MHMC) to sample values of f and A from this posterior probability distribution.We can create samples from the posterior distribution of a TNO color b − b by calculating −2.5 log 10 fb / fb for each step of the chain.This will fully capture the (often non-Gaussian) distribution of the color.Discarding the values of A in the chain is equivalent to marginalizing over variability-or one can apply a prior on A by weighting the color values by p(A).Discarding the fluxes yields the posterior distribution of the LCA values A, which is also often non-Gaussian, whenever A is poorly constrained.
The accuracy of the inferences on A depend critically on having accurate estimates of the measurement errors σ j on individual epochs.The careful characterization of image noise (Bernstein et al. 2018) and photometric calibration (Burke et al. 2017) of the DES data work together with the scene-modeling methods of Section 2 to return reliable uncertainties.We take several further steps to guard against spurious measurements that will inflate the LCA estimate.First, every exposure of every TNO is visually inspected and those with image defects (bad columns, cosmic rays, scattered light, etc.) or poor scene-modelling residuals are excluded from the analysis.After the MHMC runs, any individual photometric data points that lie > 3σ outside of the span of the model light curve at the median A value are clipped, and the MHMC is re-run.We also visually inspect diagnostic plots of the MHMC (as in Figures 6) for outlying measurements or unusual posterior distributions, which results in identification of a handful of additional measurement issues.Figure 6 illustrates this procedure for two objects.
One validation of the LCA estimation is that the brightest observed TNO, Eris, returns a value of A sharply constrained to a 68% confidence interval of 0.020-0.025.Bernstein et al. (2023) determine the period and phased light curve of Eris from the combination of DES photometry with data from 3 other observing campaigns at different observatories, and derive a sinusoidal light curve with an amplitude of A = 0.015 ± 0.001.The effect of spurious sources of fluctuation-such as calibration errors, illumination phase variation, photometric systematic errors, under-estimation of measurement errors-is to induce perhaps an extra 0.01 mag of perceived variability semi-amplitude.
Another cross-check on the determination of the light-curve amplitude is comparison of our inferred A ≈ 0.2 for (612620) 2003 SQ 317 (implying ∆m = 0.44 mag) with the ≈ 0.85 mag peak-to-peak variability reported by Lacerda et al. (2014).As illustrated the topmost panel of Figure 6, it is plausible that DES photometry sampled variations in flux as broad as Lacerda et al. (2014) observed.But the center panel of the second row shows that our estimated A is incompatible with a sinusoidal light curve with ∆m = 0.85 mag.This discrepancy is at least partially attributable to the fact that 2003 SQ 317 has a very non-sinusoidal light curve (as is characteristic of contact binaries).It is also possible that the light curve is varying over time as the system precesses or our viewing angle changes.This object has a rich long-term photometric dataset and is worthy of further investigation.In any case, the discrepancy between our A inference and the higher observed ∆m does not invalidate our methodology.We are not using the sinusoidal-assumption A as a definitive measure of peak-topeak variability; we are using it as an indicative measure of variability that we can use to compare different TNO populations.As long as we use a repeatable, well-defined measure of flux variability, such comparisons remain valid.

CONSTRAINING POPULATION-VARIANCE MODELS 5.1. Probability calculation
To make physical inferences about some selected subpopulation of TNOs, we will assume that their LCA's 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 LCA's-selections both on which objects were targeted, and on which were published.Non-detections of variation or periods often go unpublished, and even published upper limits on LCA's have not been incorporated into population analyses in a rigorous way.An advantage of the DES TNO sample is that we have multi-epoch 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 MHMC 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 This follows from the use of Bayes' Theorem, plus the fact that the samples A ik from the MHMC 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 the MHMC 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 MHMC outputs, e.g. (15)

Application to the DES TNOs
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 modelling 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 objects 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, with a probability distribution function usually written as We transform the parameters slightly, using first the mean of the distribution, Ā = a/(a + b), and a second parameter s ≡ log 10 (a + b), a "sharpness" parameter: at fixed Ā, the β distribution becomes narrower about the mean as s increases.Thus we will have 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) show that the selection  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 light curve amplitude Ā of the distribution of the subset's members' A values.
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 ( Ā, s) using Bayes' Theorem, p( Ā, s|D) ∝ p(D| Ā, s)p( Ā, 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 ( Ā, s) space numerically using Equation ( 14), and assuming uniform priors over Ā and s.We will also marginalize over s to obtain p( Ā|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( Ā, s), and Figure 7 plots the posterior probability in the ( Ā, s) space, and marginalized down to p( Ā) 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.
Figure 7 immediately indicates that the "Big" H r < 6 sample is indeed less variable (lower Ā) 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 H 1 is that they are from a single β distribution.Hypothesis 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 H 2 , distinct distributions.A frequentist-oriented statistic is   2, the top plot shows their 68% and 95% posterior credible regions of the ( Ā, s) distribution parameters as used in Equation ( 17).The maximum posterior likelihood of each is marked with an "X."The bottom panel projects these posterior distributions onto the single parameter Ā that gives the mean of the inferred distribution of noiseless light curve amplitudes for that population.It is immediately clear that larger TNOs ("Big") with H r < 6 are less variable than the smaller members of all the dynamical populations, with the possible exception of detached TNOs.Furthermore the hot classicals ("HC") are less variable than the cold classicals ("CC"), with the resonant population being consistent with either and our constraints on the scattering population being weaker.Detached TNOs are also decisively less variable than resonant or scattering TNOs in the same H r range.
This statistic compares the maximum posterior likelihood of H 2 to H 1 .In some conditions (which we do not actually meet), this statistic would have a χ 2 distribution with ν = 2 degrees of freedom, if the two populations were indeed drawn from the same model.Then the probability of ∆χ 2 > 4.6 would be 10% and ∆χ 2 > 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.By this criterion, the "Big" sample is indeed decisively distinct from all the higher-H r subpopulations, except for being indicative for the detached TNOs.
A more Bayesian statistic is the (log) evidence ratio, defined as Adopting the scale of Jeffreys (1961), values of R > 2.30, 3.45, and 4.61 are labelled as "strong," "very strong," and "decisive" in favor of H 2 .The first and last of these correspond to ≈ 10% and ≈ 1% chances of H 1 being correct, if we assign equal prior probability to H 1 and H 2 .
For our problem, we require a normalized prior p( Ā, s).We choose a prior which is uniform over the range 0.06 < Ā < 0.19, 0.7 < s < 1.9.The evidence ratio is indeed decisively in favor of the Big size distribution being distinct from all the other subsets, except the detached population, for which R = 0.54 is indecisive about H 1 vs H 2 .
Also conclusive from the ∆χ 2 test and "very strong" at R = 4.33 in the evidence ratio is the hypothesis that the cold classicals at i free < 5 • are distinct, namely more variable, than the hot classicals (i free > 5 • ).The CC's 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 labelled subset. 3The agreement is good, with deviations that are larger for the subsets with fewer members.
Next we examine the size-controlled samples of the other dynamically excited populations, the resonant, scattering, and detached TNOs.Neither the resonant and scattered TNO populations can be confidently distinguished from the HC's or CC's in terms of their q(A) parameters.Indeed the evidence ratio R = −3.53 for the HC-resonant pair indicates very strong preference for hypothesis of a shared distribution.By the ∆χ 2 criterion, there is indicative evidence that the resonants are distinct from the CCs, and that the scattering TNOs are distinct from the HC's.We can summarize this by saying that the variabilities of the resonant and scattered populations are fully consistent with being somewhere between the HC's and CC's.
The last dynamically excited population, the detached TNOs, lies between the HC and Big loci, and is decisively distinct from CC's in variability by all criteria.The evidence ratio indicates that .The best-fit β-distributions q(A) of the light-curve amplitude for the marked subsets are plotted in the lower panels as solid lines, and the corresponding cumulative distributions in the upper panels.The dotted lines are the posterior distributions for A of the observed sample when using these best-fit functions as priors.The agreement between the data and the β distributions are good, with larger deviations seen (as expected) in the samples with fewer members.The top axis converts our sinusoidal flux variation amplitude A into the peak-to-peak magnitude change of the light curve (under the sinusoidal assumption), for easier comparison to the literature.The upper axis only applies to the upper panels.
the detached TNOs could be from the same distribution as either the Big or HC, while the ∆χ 2 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% CL ellipses of two populations do not overlap in in the upper panel of Figure 7, then they are decisively distinct from each other by one or both of the more formal criteria.

Discussion
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 rigour 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: • 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.
• 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.
• Estimates of the variability yield a full probability distribution for A, rather than arbitrary decisions about and treatment of non-detections.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.
• 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 various if they have detected variation, leaving out non-detections.
• The 3-5 year temporal spread of the DES photometry is far longer than any TNO rotation period, insuring that light-curve phases are randomly sampled and our variability estimates are independent of the period.
• 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 traditional measurements of light curves, but there is still a great deal of information available from their colors and variability.

Variability distributions
The level of variability in our "Big" sample (3.3 < H r < 6.0, plus Eris) is decisively lower than that 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, find a correlation of light-curve amplitude vs H V at 3σ significance in a sample comprised of 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). 4 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 is also in general agreement with the targeted sample of primarily large objects by Kecskeméthy et al. (2023), who show 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 CC's than in non-CC's, have been noted in other samples.Benecchi & Sheppard (2013) report differences at ≈ 95% CL between the variabilities of classical TNOs and those of resonant and/or scattered objects, and detect a correlation between ∆m and i at similar confidence across the full TNO sample.The decisive difference between HC and CC populations' varibilities seen in our data reinforce the suspicions based on dynamical and color differences that the CC's 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 physical two populations, so our HC and CC samples are a bit mixed, which means we are potentially underestimating their distinction.
Thirouin & Sheppard (2019) have compared their own measures of CC variability against other populations' data in the literature, reporting that 36% of CC's 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 the dynamical classes, at fixed H range. Alexandersen et al. (2019) also report 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, 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 (∆χ 2 = 18.2, R = 6.5).The detached and resonant populations meet the frequentist criterion for decisive distinction (∆χ 2 = 12.0, 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 scattered 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.

Binary discoveries
We discover one new binary, in addition to the previously known 2014 LQ 28 , and characterize their colors and sizes.Our ability to identify these binaries are a proof-of-concept of the combination of scene modeling photometry 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 28 and 2013 RJ 124 are cold Classicals, adding to the now large sample of characterized binary CCs (Noll et al. 2020).Notably, 2014 LQ 28 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 124 is either half as dense or half as reflective as 2014 LQ 28 .Assuming the same albedo for both objects, the volume ratios between both objects is 2.8, while the mass ratio is ≈ 1.38.On the other hand, if we assume that 2014 LQ 28 's albedo is twice that of 2013 RJ 124 , 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 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.

Data release
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 kind 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 at [a site to be opened upon acceptance of this paper by the journal].Included in the data release are: • A table of all the valid individual photometric flux measurements of all the TNOs.
• A table of all photometric measurements of the resolved binaries.
• For the two TNO binaries we also include the MCMC chains with their derived mutual orbits.
• 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.
• For each TNO, a table giving the output MCMC chains produced from its photometry.
• 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 light-curve amplitudes A for each object, and derived from the MCMC chains.
• 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.
of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Figure 1 .
Figure 1.Examples of successful scene-modelling measurements for a few objects.The columns, from left to right for each object are: 1) a postage stamp of the data centered on the detection; 2) the best-fit model in Equation2; 3) the difference between the data and the background portion of the model (i.e.without subtracting the model of the TNO itself) and 4) the residuals of data minus the full model.The first row shows an example measurement of our brightest object, the magnitude m r ≈ 19.5 Eris.

Figure 2 .
Figure2.Same columns as in Figure1, but the four examples show failures of the scene-modelling procedure.The first row shows a case where the TNO was near an extended background galaxy.This example is remarkable because, even though there is no usable photometry for the TNO, the model reproduces the spiral structure of the galaxy really well, leading to a near perfect subtraction.The second row shows a similar case of a bright background source near the TNO, but in this case the model failed to reproduce this background.The third row corresponds to a bright source near the edge of the stamp, with its scattered light precluding any photometric measurement nearby.The final row is a case where the least squares procedure of Equation 3 led to an incorrect solution (namely, the model predicts a large negative flux for the central point source).

Figure 3 .
Figure3.Single (upper row of each object) and binary (bottom row) point source scene-modelling photometry applied to the same detections of two high-confidence binary TNOs identified in the DES data.In the first case (∆χ 2 = 130.6), a distinct quadrupole pattern is present in the residual image, a strong indication that this source is a binary.The second object (∆χ 2 = 12.8), while less visually apparent, is statistically significant across several images.

Figure 4 .
Figure 4. ∆χ 2 between the single source and binary models vs S/N of the fainter source for a randomly chosen subset of 10% of the ≈ 25, 000 photometric measurements (in gray) in the DES data.The two binary objects are shown with different colors, and the acceptance region is shaded in yellow.The other cases of objects where the binary fit yields a statistically significant improvement in the χ 2 arise when there is a poorly subtracted background object in the single-source fit, or when the binary fitting corrects small errors in the assumed pixel position of the primary object caused by atmospheric turbulence.
Figure 5. g − r (blue), r − i (red), i − z (green) and g − i (purple) colors and corresponding uncertainties measured from sequential exposures of each component of the 2014 LQ 28 system.The dashed gray diagonal line corresponds to the regime where the colors are equal in both members.The colors scatter around this diagonal line, and are all within 3σ of the nominal identity, implying that each member has the same colors (and therefore, similar surface composition).

Figure 6 .
Figure 6.The results of the modeling of the light curve amplitude and determination of mean flux and color are summarized for (612620) 2003 SQ 317 (top), a source with well-characterized variability; and 2016 SD 106 (bottom), a source consistent with no variability.For each source, the upper panel plots the measured fluxes over the six years of the survey, in particular the flux variation f /f b relative to the mean flux in each band.The bottom panels for each object shows histograms of the MCMC samples drawing from the inferred posterior probabilities of the absolute magnitude H r (left), the light curve amplitude (LCA) A (middle), and the colors (right) for the object.The horizontal band in the relative-flux figure and the dashed line in the LCA figure for (612620) 2003 SQ 317 show the light curve peak-to-peak amplitude measured byLacerda et al. (2014).Our method underestimates the LCA, likely because the true light curve is not sinusoidal.The vertical bands in the 2016 SD 106 color histograms show the 68% limits fromChen et al. (2022), which are in good agreement with our inferred colors.

Figure 7 .
Figure 7.For each of six disjoint subsamples of TNOs described in Table2, the top plot shows their 68% and 95% posterior credible regions of the ( Ā, s) distribution parameters as used in Equation (17).The maximum posterior likelihood of each is marked with an "X."The bottom panel projects these posterior distributions onto the single parameter Ā that gives the mean of the inferred distribution of noiseless light curve amplitudes for that population.It is immediately clear that larger TNOs ("Big") with H r < 6 are less variable than the smaller members of all the dynamical populations, with the possible exception of detached TNOs.Furthermore the hot classicals ("HC") are less variable than the cold classicals ("CC"), with the resonant population being consistent with either and our constraints on the scattering population being weaker.Detached TNOs are also decisively less variable than resonant or scattering TNOs in the same H r range.
Figure8.The best-fit β-distributions q(A) of the light-curve amplitude for the marked subsets are plotted in the lower panels as solid lines, and the corresponding cumulative distributions in the upper panels.The dotted lines are the posterior distributions for A of the observed sample when using these best-fit functions as priors.The agreement between the data and the β distributions are good, with larger deviations seen (as expected) in the samples with fewer members.The top axis converts our sinusoidal flux variation amplitude A into the peak-to-peak magnitude change of the light curve (under the sinusoidal assumption), for easier comparison to the literature.The upper axis only applies to the upper panels.

Table 1 .
Mutual orbits of the trans-Neptunian binaries

Table 2 .
TNO sample definitions 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 classification are taken from Bernardinelli et al. (