# Photometry of outer Solar System objects from the Dark Energy Survey I: photometric methods, light curve distributions and trans-Neptunian binaries

Pedro H Bernardinelli, Gary M Bernstein, Nicholas Jindal, T.M.C Abbott, M Aguena, F Andrade-Oliveira, J Annis, D Bacon, E Bertin, D Brooks, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Pedro H Bernardinelli, Gary M Bernstein, Nicholas Jindal, T.M.C Abbott, M Aguena, et al. Photometry of outer Solar System objects from the Dark Energy Survey I: photometric methods, light curve distributions and trans-Neptunian binaries. Astrophys.J.Suppl., 2023, 269 (1), pp.18. 10.3847/1538-4365/acf6bf. hal-04080097

HAL Id: hal-04080097

https://hal.science/hal-04080097

Submitted on 10 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

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

Pedro H. Bernardinelli<sup>1,43</sup>, Gary M. Bernstein<sup>2</sup>, Nicholas Jindal<sup>2,3</sup>, T. M. C. Abbott<sup>4</sup>, M. Aguena<sup>5</sup>, O. Alves<sup>6</sup>, F. Andrade-Oliveira<sup>6</sup>, J. Annis<sup>7</sup>, D. Bacon<sup>8</sup>, E. Bertin<sup>9,10</sup>, D. Brooks<sup>11</sup>, D. L. Burke<sup>12,13</sup>, A. Carnero Rosell<sup>5,14,15</sup>, M. Carrasco Kind<sup>16,17</sup>, J. Carretero<sup>18</sup>, L. N. da Costa<sup>5</sup>, M. E. S. Pereira<sup>19</sup>, T. M. Davis<sup>20</sup>, S. Desai<sup>21</sup>, H. T. Diehl<sup>7</sup>, P. Doel<sup>11</sup>, S. Everett<sup>22</sup>, I. Ferrero<sup>23</sup>, D. Friedel<sup>16</sup>, J. Frieman<sup>7,24</sup>, J. García-Bellido<sup>25</sup>, G. Giannini<sup>18</sup>, D. Gruen<sup>26</sup>, K. Herner<sup>7</sup>, S. R. Hinton<sup>20</sup>, D. L. Hollowood<sup>27</sup>, K. Honscheid<sup>3,28</sup>, D. J. James<sup>29</sup>, K. Kuehn<sup>30,31</sup>, J. L. Marshall<sup>32</sup>, J. Mena-Fernández<sup>33</sup>, F. Menanteau<sup>16,17</sup>, R. Miquel<sup>18,34</sup>, R. L. C. Ogando<sup>35</sup>, A. Palmese<sup>36</sup>, A. Pieres<sup>5,35</sup>, A. A. Plazas Malagón<sup>37</sup>, M. Raveri<sup>38</sup>, E. Sanchez<sup>33</sup>, I. Sevilla-Noarbe<sup>33</sup>, M. Smith<sup>39</sup>, E. Suchyta<sup>40</sup>, M. E. C. Swanson<sup>41</sup>, G. Tarle<sup>6</sup>, C. To<sup>28</sup>, A. R. Walker<sup>4</sup>, P. Wiseman<sup>39</sup>, and Y. Zhang<sup>4,42</sup>

(The DES Collaboration)

<figure><img src="image_4.png" /></figure>

<footnote n="1"><sup>1</sup> DiRAC Institute, Department of Astronomy, University of Washington, 3910 15th Avenue NE, Seattle, WA 98195, USA; phbern@uw.edu</footnote>

<footnote n="2"><sup>2</sup> Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA; garyb@upenn.edu</footnote>

<footnote n="3"><sup>3</sup> Department of Physics, Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="4"><sup>4</sup> Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile</footnote>

<footnote n="5"><sup>5</sup> Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="6"><sup>6</sup> Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="7"><sup>7</sup> Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="8"><sup>8</sup> Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK</footnote>

<footnote n="9"><sup>9</sup> CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="10"><sup>10</sup> Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="11"><sup>11</sup> Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK</footnote>

<footnote n="12"><sup>12</sup> Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="13"><sup>13</sup> SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="14"><sup>14</sup> Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="15"><sup>15</sup> Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="16"><sup>16</sup> Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA</footnote>

<footnote n="17"><sup>17</sup> Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="18"><sup>18</sup> Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona) Spain</footnote>

<footnote n="19"><sup>19</sup> Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="20"><sup>20</sup> School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="21"><sup>21</sup> Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="22"><sup>22</sup> Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="23"><sup>23</sup> Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="24"><sup>24</sup> Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="25"><sup>25</sup> Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="26"><sup>26</sup> University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="27"><sup>27</sup> Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="28"><sup>28</sup> Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="29"><sup>29</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="30"><sup>30</sup> Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="31"><sup>31</sup> Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA</footnote>

<footnote n="32"><sup>32</sup> George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote n="33"><sup>33</sup> Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="34"><sup>34</sup> Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="35"><sup>35</sup> Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="36"><sup>36</sup> Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15312, USA</footnote>

<footnote n="37"><sup>37</sup> Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="38"><sup>38</sup> Department of Physics, University of Genova and INFN, Via Dodecaneso 33, I-16146, Genova, Italy</footnote>

<footnote n="39"><sup>39</sup> School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK</footnote>

<footnote n="40"><sup>40</sup> Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA</footnote>

<footnote n="41"><sup>41</sup> Independent Researcher</footnote>

<footnote n="42"><sup>42</sup> Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote>Received 2023 April 10; revised 2023 August 28; accepted 2023 August 29; published 2023 October 31</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further</footnote>

<footnote>distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.</footnote>

<footnote n="43"><sup>43</sup> DiRAC Postdoctoral Fellow.</footnote>
Abstract

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.

Unified Astronomy Thesaurus concepts: Kuiper belt (893); Trans-Neptunian objects (1705); Photometry (1234); Asteroid satellites (2207)

1. 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 hundred 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, for a recent review, Gladman & Volk 2021). 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; Sánchez 2006; The Dark Energy Survey Collaboration 2016) received an allocation of 575 nights on the 4 m Blanco telescope at Cerro Tololo, using the Dark Energy Camera (DECam; Flaugher et al. 2015) to cover 5000 deg² of the sky in the *grizY* photometric system from 2013–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, 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 (LCA).

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, and ensures that the vast majority of objects have high-quality orbital arcs spanning multiple years (Bernardinelli et al. 2022). Each object has been observed many times: typically between 6 and 10 times in each of the *grizY* bands, over the 6 yr of data collection. This means that each object is imaged typically 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 for estimation of colors in the griz bands (*Y* is generally too low in signal-to-noise ratio, S/N) 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 nonsimultaneous 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 optimal fluxes from the pixel data in the images, and this process 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. Section 4 describes the extraction of colors and LCAs from the flux time series. In Section 5 we examine the distributions of LCAs in different TNO subpopulations as an indicator of different physical states. Scientific analyses of the color catalog will appear in future publications.

Accompanying this publication is a data release with the ≈30,000 photometric measurements, absolute magnitudes, colors, and LCAs of our objects, as well as some of the software required for the analysis we present here. The data release is available on Zenodo at doi:10.5281/zenodo.8231238, and discussed in Section 6.4.
## 2. 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 that the orbit predicts should contain the object, but there is no detection in the catalog (“nondetec-tions”)—similar to the subthreshold significance 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 observation using the scene modeling photometry (SMP) technique commonly used in supernovae type Ia cosmology (Brout et al. 2019). This technique enables the optimal determination of fluxes and uncertainties of the target in the presence of background sources. Similarly to difference imaging, this technique makes use of all images taken by the survey in the same band for a given region of the sky, being well-suited to a survey that serendipitously observes moving targets like DES. Unlike difference imaging, however, SMP does not require a template image to be created, which makes it very computationally efficient: as will be seen below, our SMP methodology requires only the solution of a linear system of equations to obtain each TNO flux measurement. 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. (2021), and comparisons between DES and Gaia show an exquisite calibration with 3 mmag rms 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 an $m \times m$ square grid, with spacing of roughly the $1\sigma$ 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'' \times 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 $\mu$, $M_{ij}^{\mu}$, for this mosaic of sources is modeled as the convolution of each $P_{uv}$ with exposure's PSF derived for its location $(u, v)$; plus some constant background level $b_{\mu}$:

$$ M_{ij}^{\mu} = \sum_{u,v} PSF_{\mu}(i-u, j-v)P_{uv} + b_{\mu}. \quad (1) $$

In the single exposure $\nu$ where the TNO is present, we adopt the same model with an additional point-source term:

$$ M_{ij}^{\nu} = \sum_{u,v} PSF_{\nu}(i-u, j-v)P_{uv} + b_{\nu} + f_{\text{TNO}} PSF_{\nu}(i, j). \quad (2) $$

That is, $f_{\text{TNO}}$ represents the integrated flux of the TNO at this epoch and band. For the maximal rate of motion for TNOs of

~5" hr$^{-1}$, an object moves 0.125 in our 90 s exposure times, almost a factor of 8 below the median 0.98 PSF FWHM for r-band exposures (Sevilla-Noarbe et al. 2021), so these sources are not trailed significantly. This means that we can safely ignore corrections such as pill apertures (e.g., Fraser et al. 2016). 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 \times 30$ pixels, corresponding to a sky area of $7'' \times 7''$). These choices mean that the model does not span the full postage stamp, to avoid fitting for potential contributions from sources on the edges or whose center lies outside the stamp without its full data.

We fit the model using a least-squares minimization, comparing the model to the measured pixel values $Im^{\mu}$:

$$ \chi^2 = \sum_{i,j} \frac{(k_{\nu} Im_{ij}^{\nu} - M_{ij}^{\nu})^2}{\sigma_{\nu,ij}^2} + \sum_{\mu,ij} \frac{(k_{\mu} Im_{i,j}^{\mu} - M_{ij}^{\mu})^2}{\sigma_{\mu,ij}^2}. \quad (3) $$

The (constant) $k_{\mu,\nu}$ terms correct for the different zero-points in each exposure, bringing them to a common flux scale (Burke et al. 2017), and $\sigma_{\mu,ij}^2$ 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_{\mu}, f_{\text{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. Note that if a background source has variable flux (e.g., a variable star), its corresponding $P_{uv}$ will represent the mean measured flux, and such a case might lead to a poor subtraction of that source.

The initial estimates of the $\sigma_{\nu,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 $\nu$ using the first least-squares solution for $M_{ij}^{\nu}$:

$$ \sigma_{\nu,ij}^2 \rightarrow \sigma_{\nu,ij}^2 + \frac{\max(M_{ij}^{\nu}, 0)}{g_{\nu,ij}}, \quad (4) $$

where $g_{\nu,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,44 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-modeling 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 \times 400$ pixel stamps in order to detect C/2014 UN271's

<footnote n="44">44 All DES PSF and photometric color dependence is parameterized by $g-i$ of a stellar source.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Examples of successful scene-modeling 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 Equation (2); (3) the difference between the data and the background portion of the model (i.e., without subtracting the model of the TNO itself); (4) the difference between the data and the point-source portion of the model (i.e., without subtracting the background model); and (5) the residuals of data minus the full model. The first row shows an example measurement of our brightest object, the magnitude $m_r \approx 19.5$ Eris.</figcaption></figure>

extended coma. A generalization to binary sources is presented in the next section.

### 3. Identifying Binaries with Scene-modeling Photometry

#### 3.1. Methodology

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

$$ M_{ij}^{\nu} = \sum_{u,v} \mathbf{PSF}_{\nu}(i-u, j-v)\mathbf{P}_{uv} + \mathbf{b}_{\nu} \\ + f_1 \mathbf{PSF}_{\nu}(i + \Delta x_1, j + \Delta y_1) \\ + f_2 \mathbf{PSF}_{\nu}(i + \Delta x_2, j + \Delta y_2), \quad (5) $$

where the indices 1 and 2 correspond to the primary and secondary sources, respectively. For a fixed set of {$\Delta x_1, \Delta y_1, \Delta x_2, \Delta y_2$}, the linearity of the least-squares procedure in Equation (3) is preserved.<fnref n="45" /> The minimization of the $\chi^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 $\chi^2$ between the two, $\Delta\chi^2 = \chi^2_{\text{single}} - \chi^2_{\text{binary}}$. As we have an additional five parameters (the flux of the secondary source and the position shifts), we expect $\Delta\chi^2$ to also follow a $\chi^2$ distribution with 5 degrees of freedom. However, we note that there are other reasons (other than a binary object) for the $\chi^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 photometry), or another transient in the same image (such as an asteroid). We visually inspect all cases where $\Delta\chi^2 \ge 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 ≈25,000 of our detections in the griz bands leads to two objects where several images led to both an improvement in the $\chi^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 LQ28 had been previously identified as a binary (Thirouin & Sheppard 2019), 2013 RJ124 is a new binary discovery. The objects (612620) 2003 SQ317, 2014 QL441, and 2016 TT94 also show two images each with $\Delta\chi^2 > 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

<footnote n="45">⁴⁵ Note that this is also the case for a single source with a free centroid, and so the same methodology can be used to adjust the location of the PSF. This is particularly useful for bright sources, where ≈10 mas shifts in the center can account for uncertainties due to atmospheric turbulence.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Same columns as in Figure 1, but the four examples show failures of the scene-modeling procedure. The first row shows a case where the TNO (with a small expected S/N) was near an extended background galaxy with a much more significant flux contribution. 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 rare 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), possibly due to excess noise fluctuations in the background and the low S/N of the target source.</figcaption></figure>

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.

## 3.2. Mutual Orbit Determination and Derived Properties

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 $\Delta\chi^2 \ge 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 $\chi^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 $\phi \in [0^\circ, 180^\circ]$, 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, \Omega_m, \omega_m, 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 $\mu$, we use a Markov Chain Monte Carlo (MCMC) to sample this seven-dimensional parameter space

with the likelihood

$$ \prod_{\mu} \exp \left[ - \frac{\mathcal{L}_{\infty} (r_{\mu}, \hat{r}(t_{\mu}), \phi, \hat{\phi}_{\mu}) \sum_{\mu}^{-1} (r_{\mu}, \hat{r}(t_{\mu}), \phi, \hat{\phi}_{\mu})^T}{2} \right] \quad (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 \le \Omega_m, \omega_m, M_m < 2\pi, 0 \le i_m \le \pi, 0 \le e_m \le 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, \Omega_m, \omega_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 well as the ratio between the mutual semimajor axis and the system’s Hill radius (see, e.g., Parker et al. 2011).

In addition to astrometric data, we can also determine flux ratios. We define the magnitude difference $\delta m \equiv |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
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Single (upper row of each object) and binary (bottom row) point-source scene-modeling photometry applied to the same detections of two high-confidence binary TNOs identified in the DES data. In the first case ($\Delta\chi^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 ($\Delta\chi^2 = 12.8$), while less visually apparent, is statistically significant across several images.</figcaption></figure>

colors without any ambiguity. These results are included in the data release (see Section 6.4).

2014 LQ28 has several observation pairs, presented in Figure 5. All color pairs are less than $3\sigma$ 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.

## 4. Extracting Colors and Light-curve Amplitudes

After the scene-modeling photometry, each TNO $i \in \{1, 2, ..., N_{\text{TNO}}\}$ has a series of measured fluxes $f_{ij}$ for $j \in \{1, ..., N_{\text{obs},i}\}$ observations in bands $b_{ij}$, with nearly Gaussian uncertainties $\sigma_{ij}$. We first adjust all fluxes for their heliocentric and geocentric distances to represent fluxes that would be observed with both distances at $d_{\text{ref}} = 30$ au. In the derivation below, we assume that this adjustment has been made for all fluxes. What are the best estimates of the mean fluxes $\bar{f}_{ib}$ for source $i$ in band $b$? From these, the best estimated absolute magnitudes $H_{ib}$ and colors $(b - b')_i \equiv H_{ib} - H_{ib'}$ can be determined, as well as 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 for 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 $\phi$, and get phase-averaged fluxes $\bar{f}_{ib}$ for each band $b$. Unfortunately the DES observing cadence does not admit either method. In each band, the source is typically observed eight times over a 5 yr span, far too sparse to determine a period, never mind construct a light curve. Furthermore only occasionally are two observations of a given TNO 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
<figure><img src="image_8.png" /><figcaption>**Figure 4.** $\Delta\chi^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 $\chi^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.</figcaption></figure>

in which the flux at observation j is assumed to equal

$$ \hat{f}_j = \bar{f}_{b_j} [1 + Ah(\phi_j)]. \quad (7) $$

Here we have dropped the index i of the TNO, for brevity, since the process is independent for each TNO. A phase function $h(\phi)$ is a function of the variability phase $\phi \in [0, 1)$ with $\langle h \rangle = 0$ when averaged over phase, and $\max_{\phi} h - \min_{\phi} h = 2$. The parameter A then gives the semiamplitude of the fractional variation of true flux, and $\bar{f}_b$ is the time-averaged mean flux in band b. Note that we define the mean and the light-curve semiamplitude in flux space rather than in magnitude space. The peak-to-peak magnitude variation would be $\Delta m \equiv 2.5 \log_{10}[(1 + A)/(1 - A)]$. We assume the convention $A \ge 0$, and only $A < 1$ is physically possible. We make the following critical assumptions:

1. All bands have the same phase function h and amplitude A, i.e., the variation is achromatic.

2. The variability phase $\phi_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 hr of each other, in which case we assume $\phi_j = \phi_k$ but the value is unknown.

3. We will assume that the phase function is $h(\phi) = \sin(2\pi\phi)$. This is likely inaccurate for non-ellipsoidal geometries, but our goal is to obtain A values that are indicative of variability level of each object. To admit comparison of populations' variabilities, we do not 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 $\sigma_j$ as

$$ p(\{f_j\} | \{\bar{f}_b\}, A, \{\phi_j\}) \propto \prod_j \exp\left\{-\frac{1}{2}\left(\frac{f_j - \bar{f}_{b_j}[1 + Ah(\phi_j)]}{\sigma_j}\right)^2\right\}. \quad (8) $$

To obtain $p(\{f_j\} | \{\bar{f}_b\}, A)$ we need to marginalize Equation (8) over all possible phases using the second assumption. We divide all of 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 1 hr of each other.<fnref n="46" /> With this convention, the marginalization over light-curve phases becomes

$$ p(\{f_j\} | \{\bar{f}_b\}, A) \propto \prod_s \int_0^1 d\phi_s \prod_{j \in s} \exp\left\{-\frac{1}{2}\left(\frac{f_j - \bar{f}_{b_j}[1 + Ah(\phi_s)]}{\sigma_j}\right)^2\right\}. \quad (9) $$

Numerically, the integral over $\phi$ can be executed as a sum over ≈20 equally spaced samples of $0 \le \phi < 1$. Using Bayes' theorem and assuming uniform priors for $\bar{f}_b$ and for $0 \le A < 1$, we can take

$$ p(\{\bar{f}_b\}, A | \{f_j\}) \propto p(\{f_j\} | \{\bar{f}_b\}, A). \quad (10) $$

We created a straightforward Metropolis-Hastings MCMC to sample values of $\bar{f}$ and A from this posterior probability distribution. We can create samples from the posterior distribution of a TNO color $b - b'$ by calculating

<footnote n="46">⁴⁶ In the DES wide survey, fields are not visited multiple times per night except in short time intervals.</footnote>
<table><caption>**Table 1**
Mutual Orbits of the Trans-Neptunian Binaries</caption><thead><tr><th>Object</th><th>a<sub>m</sub><br>(10<sup>6</sup> m)</th><th>e<sub>m</sub></th><th>i<sub>m</sub><br>(deg)</th><th>Ω<sub>m</sub><br>(deg)</th><th>ω<sub>m</sub><br>(deg)</th><th>M<sub>m</sub><br>(deg)</th><th>P<sub>m</sub><br>(days)</th><th>Epoch<br>(MJD)</th><th>Mass<br>(10<sup>18</sup> kg)</th><th>Hill Radius<br>(10<sup>7</sup> m)</th></tr></thead><tbody><tr><td>2014 LQ<sub>28</sub></td><td>32.9<sup>+7.1</sup><sub>-3.5</sub></td><td>0.64<sup>+0.14</sup><sub>-0.16</sub></td><td>43<sup>+2</sup><sub>-2</sub></td><td>5<sup>+5</sup><sub>-3</sub></td><td>256<sup>+4</sup><sub>-6</sub></td><td>268<sup>+18</sup><sub>-20</sub></td><td>1470<sup>+390</sup><sub>-230</sub></td><td>57,615.258</td><td>1.4<sup>+0.7</sup><sub>-0.5</sub></td><td>52<sup>+7</sup><sub>-7</sub></td></tr><tr><td>2013 RJ<sub>124</sub></td><td>34.9<sup>+6.6</sup><sub>-4.0</sub></td><td>0.14<sup>+0.13</sup><sub>-0.09</sub></td><td>90<sup>+14</sup><sub>-15</sub></td><td>103<sup>+49</sup><sub>-40</sub></td><td>152<sup>+162</sup><sub>-43</sub></td><td>38<sup>+121</sup><sub>-22</sub></td><td>1850<sup>+290</sup><sub>-250</sub></td><td>56,545.355</td><td>1.0<sup>+0.3</sup><sub>-0.2</sub></td><td>49<sup>+5</sup><sub>-4</sub></td></tr></tbody></table>

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

−2.5 log$_{10} \bar{f}_b / \bar{f}_{b'}$ 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 $\sigma_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-modeling residuals are excluded from the analysis. After the MCMC runs, any individual photometric data points that lie >3$\sigma$ outside of the span of the model light curve at the median A value are clipped, and the MCMC is rerun. We also visually inspect diagnostic plots of the MCMC (as in Figure 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 “stress test” 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) determined the period and phased light curve of Eris from the combination of DES photometry in the g and r bands with data from three other observing campaigns at different observatories, and derive a sinusoidal light curve with an amplitude of $A = 0.015 \pm 0.001$. The analysis of Bernstein et al. (2023) also showed that the i- and z-band photometry from DES do not match this sinusoidal light curve, and so we can speculate that the effect of spurious sources of fluctuation—such as calibration errors, illumination phase variation, photometric systematic errors, underestimation of measurement errors—is to induce perhaps an extra 0.01 mag of perceived variability semiamplitude in this case.

We also compare the measured A to the $\Delta m$ estimated by the survey of Thirouin & Sheppard (2019) for the three objects measured by both surveys. 2001 QS$_{322}$ and 2003 QY$_{111}$ were measured to have $\Delta m \ge 0.3$ and 0.2, respectively, by Thirouin & Sheppard (2019). The DES estimated $\Delta m$ for these objects have their 68% confidence limits constrained to 0.26–0.57, and 0.039–0.288, respectively, showing general agreement between the two surveys. 2014 OM$_{394}$ was determined to have $\Delta m \approx 0.1$ in 5.5 hr of observations by Thirouin & Sheppard (2019), while our 68% confidence limits are 0.41–0.78. This suggests that the light curve of this object varies over much

longer timescales than the 5.5 hr span of Thirouin & Sheppard (2019): follow-up observations are required to solve this discrepancy.

Another stress test on the determination of the LCA is comparison of our inferred $A \approx 0.2$ for (612620) 2003 SQ$_{317}$ (implying $\Delta m = 0.44$ mag) with the $\Delta m \approx 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 $\Delta m = 0.85$ mag. This discrepancy is at least partially attributable to the fact that 2003 SQ$_{317}$ has a very nonsinusoidal light curve (Snodgrass et al. 2010; Lacerda et al. 2014), 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 data set and is worthy of further investigation. In any case, the discrepancy between our A inference and the higher observed $\Delta m$ does not invalidate our methodology. We are not using the sinusoidal-assumption A as a definitive measure of peak-to-peak 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.

# 5. Constraining Population-variance Models

## 5.1. Probability Calculation

To make physical inferences about some selected subpopulation of TNOs, we will assume that their LCAs are drawn independently from some distribution $q(A|\theta_A)$, where $\theta_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
<figure><img src="image_9.png" /><figcaption>**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 LQ28 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).</figcaption></figure>

posterior indicates a clear detection of variability. We can now form a posterior likelihood for the distribution of the parameters $\theta_A$ of the LCA distribution for this class as

$$p(\theta_A|D) \propto p(D|\theta_A)p(\theta_A) \quad (11)$$

$$= p(\theta_A) \prod_i p(\{f_{ij}\}|\theta_A) \quad (12)$$

$$= p(\theta_A) \prod_i \int dA_i p(\{f_{ij}\}|A_i) q(A_i|\theta_A) \quad (13)$$

$$\propto p(\theta_A) \prod_i \left[ \frac{1}{N_{\text{samp},i}} \sum_{k=1}^{N_{\text{samp},i}} q(A_{ik}|\theta_A) \right]. \quad (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|\theta_A)$ to the MCMC outputs, e.g.,

$$p(b - b' | \{f_{ik}\}, \theta_A)$$

$$\propto \sum_{k=1}^{N_{\text{samp},i}} q(A_{ik}|\theta_A) \delta(b - b' + 2.5 \log_{10} \bar{f}_{ibk} / \bar{f}_{ib'k}). \quad (15)$$

## 5.2. 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 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 \le A < 1$ is the $\beta$ distribution, whose probability distribution function depends on two shape parameters $\alpha'$ and $\beta'$ and is usually written as

$$q(A|\alpha', \beta') \propto A^{\alpha'-1}(1-A)^{\beta'-1}. \quad (16)$$

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

$$q(A|\theta_A) = q(A|\bar{A}, s) \propto A^{\bar{A}10^{s-1}}(1-A)^{(1-\bar{A})10^{s-1}}. \quad (17)$$

The $\beta$ 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,
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The results of the modeling of the LCA and determination of mean flux and color are summarized for (612620) 2003 SQ317 (top), a source with well-characterized variability; and 2016 SD106 (bottom), a source consistent with no variability. For each source, the upper panel plots the measured fluxes over the 6 yr 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 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 SQ317 show the light curve peak-to-peak amplitude measured by Lacerda et al. (2014). Our method underestimates the LCA, likely because the true light curve is not sinusoidal. The vertical bands in the 2016 SD106 color histograms show the 68% limits from Chen et al. (2022), which are in good agreement with our inferred colors.</figcaption></figure>

$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
<table><caption>**Table 2**
TNO Sample Definitions</caption><thead><tr><th>Subset</th><th>N</th><th>Size</th><th>Median (H<sub>r</sub>)</th><th>Dynamical Class</th><th>Inclination</th><th>A&#773; 68% CL</th></tr></thead><tbody><tr><td>Cold classical (CC)</td><td>95</td><td>6 &lt; H<sub>r</sub> &lt; 8.2</td><td>7.18</td><td>Classical</td><td>i<sub>free</sub> &lt; 5°</td><td>0.149–0.165</td></tr><tr><td>Hot classical (HC)</td><td>261</td><td>6 &lt; H<sub>r</sub> &lt; 8.2</td><td>7.38</td><td>Classical</td><td>i<sub>free</sub> &gt; 5°</td><td>0.124–0.128</td></tr><tr><td>Resonant</td><td>170</td><td>6 &lt; H<sub>r</sub> &lt; 8.2</td><td>7.56</td><td>Resonant</td><td>...</td><td>0.127–0.136</td></tr><tr><td>Detached</td><td>133</td><td>6 &lt; H<sub>r</sub> &lt; 8.2</td><td>7.39</td><td>Detached</td><td>...</td><td>0.114–0.115</td></tr><tr><td>Scattering</td><td>34</td><td>6 &lt; H<sub>r</sub> &lt; 8.2</td><td>7.57</td><td>Scattering</td><td>...</td><td>0.148–0.176</td></tr><tr><td>Big</td><td>36</td><td>H<sub>r</sub> &lt; 6</td><td>5.43</td><td>...</td><td>...</td><td>0.081–0.095</td></tr></tbody></table>

**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_{\text{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.

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 $\langle H_r \rangle = 7.4 \pm 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 $\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 $H_1$ is that they are from a single $\beta$ distribution. Hypothesis $H_2$ is that they are from distinct distributions with two pairs of $\beta$-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

$$ \Delta\chi^2 \equiv 2 \log \frac{[\max_{\bar{A},s} p(\bar{A}, s|D_A)][\max_{\bar{A},s} p(\bar{A}, s|D_B)]}{\max_{\bar{A},s} p(\bar{A}, s|D_A, D_B)}. \quad (18) $$

This statistic compares the maximum posterior likelihood of $H_2$ to $H_1$. Under the assumption that the parameters are Gaussian distributed, this statistic would have a $\chi^2$ distribution with $\nu = 2$ degrees of freedom, if the two populations were indeed drawn from the same model. Then the probability of $\Delta\chi^2 > 4.6$ would be 10% and $\Delta\chi^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. 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{aligned} \mathcal{R} &= \log \frac{p(D|H_2)}{p(D|H_1)} \\ &= \log \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{aligned} \quad (19) $$

Adopting the scale of Jeffreys (1961), values of $\mathcal{R} > 2.30$, 3.45, and 4.61 are labeled as “strong,” “very strong,” and “decisive” in favor of $H_2$. The first and last of these correspond to $\approx 10\%$ and $\approx 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(\bar{A}, s)$. We choose a prior that is uniform over the range $0.06 < \bar{A} < 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 $H_1$ vs $H_2$.

Also conclusive from the $\Delta\chi^2$ test and “very strong” at $\mathcal{R} = 4.33$ in the evidence ratio is the hypothesis that the CCs at $i_{\text{free}} < 5^\circ$ are distinct, namely more variable, than the HCs ($i_{\text{free}} > 5^\circ$). 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 $\beta$ distributions are adequate descriptions of the measured size distributions. In the left panel of Figure 8, we plot in solid lines the $q(A)$ $\beta$-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)/\int dA p(D_i|A)q(A)$ over all members of the labeled subset.<fnref n="47" /> The 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 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 $\Delta\chi^2$ 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.

<footnote n="47">⁴⁷ Note that this posterior distribution would equal the observed distribution if every TNO's value of A were well constrained.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** For each of six disjoint subsamples of TNOs described in Table 2, the top plot shows their 68% and 95% posterior credible regions of the ($\bar{A}$, 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 $\bar{A}$ that gives the mean of the inferred distribution of noiseless LCAs for that population. It is immediately clear that larger TNOs (Big) with $H_r < 6$ are less variable than the smaller members of all of the dynamical populations, with the possible exception of detached TNOs. Furthermore, the hot classicals (HCs) 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.</figcaption></figure>

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 $\Delta\chi^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% 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.

### 5.3. 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
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The best-fit β-distributions $q(A)$ of the LCA for the marked subsets are plotted in the lower panels as solid lines, and the corresponding cumulative distributions are shown 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 is 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.</figcaption></figure>

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.

# 6. Conclusions

## 6.1. Methodologies

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.

## 6.2. *Variability Distributions*

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 $\Delta m$ > 0.3 mag among $H \lesssim 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\sigma$ significance in a sample comprising 32 objects at $H_v < 6$ that they measure, combined with 96 objects at $0 < H_v < 12$ with $\Delta 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 $\sigma_{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).$^{48}$ 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 \le 6$ are also in general agreement with the targeted sample of primarily large objects by Kecskeméthy et al. (2023), who showed that $\approx 55\%$ of their sample has $\Delta m \le 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 $\approx 95\%$ CL between the variabilities of classical TNOs and those of resonant and/or scattering objects, and detected a correlation between $\Delta 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 $\Delta 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 $\Delta m \le 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 ($\Delta\chi^2 = 18.2, R = 6.5$). The detached and resonant populations meet the frequentist criterion for decisive distinction ($\Delta\chi^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 $\approx 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.

## 6.3. *Binary Discoveries*

We discover one new binary, in addition to the previously known 2014 LQ28, 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 \approx 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 LQ28 and 2013 RJ124 are cold Classicals, adding to the now large sample of characterized binary CCs (Noll et al.

<footnote>$^{48}$ We focus on the $\sigma_{mag}$ statistics from Alexandersen et al. (2019) because it is less susceptible to measurement noise and outliers, such as from background sources, than their $\Delta_{mag}$ statistic of extremal measurements. Thirouin & Sheppard (2019) suggested that some of the outlying OSSOS photometry needs to be clipped, and the methodology also has no allowance for removal of measurement noise. Showalter et al. (2021) found the OSSOS distributions discrepant from other data, but the other data come from ill-defined samples.</footnote>
2020). Notably, 2014 LQ28 is the second-largest CC in the DES sample ($H_r \approx 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 RJ124 is either half as dense or half as reflective as 2014 LQ28. Assuming the same albedo for both objects, the volume ratio between both objects is 2.8, while the mass ratio is $\approx 1.38$. On the other hand, if we assume that 2014 LQ28's albedo is twice that of 2013 RJ124, 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.

## 6.4. 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 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 *Jupiter* 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 *Jupiter* 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 *Jupiter* notebook with some examples of its usability, in a different Zenodo repository at doi:10.5281/zenodo.8231232 (Bernardinelli et al. 2023).

## Acknowledgments

We thank the anonymous referee for providing a detailed review of our manuscript, and the AAS data editor for providing assistance with our data release. P.H.B. acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. Work by G.M.B., P.H.B., and N.J. was supported by National Science Foundation grants AST-2009210 and AST-2205808.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical-Infrared Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and
306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project No. 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 nonexclusive, 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.

*Software:* This work made use of the following public codes: NUMPY (Oliphant 2015), SCIPY (Jones 2001), ASTROPY (Astropy Collaboration 2013, 2018), MATPLOTLIB (Hunter 2007), IPYTHON (Pérez & Granger 2007), WCSFIT and PIXMAPPY (Bernstein et al. 2017), PIFF (Jarvis et al. 2020), EMCEE (Foreman-Mackey et al. 2013).

## ORCID iDs

Pedro H. Bernardinelli https://orcid.org/0000-0003-0743-9422

Gary M. Bernstein https://orcid.org/0000-0002-8613-8259

Nicholas Jindal https://orcid.org/0000-0002-2092-3545

T. M. C. Abbott https://orcid.org/0000-0003-1587-3931

M. Aguena https://orcid.org/0000-0001-5679-6747

O. Alves https://orcid.org/0000-0002-7394-9466

J. Annis https://orcid.org/0000-0002-0609-3987

D. Bacon https://orcid.org/0000-0002-2562-8537

E. Bertin https://orcid.org/0000-0002-3602-3664

D. Brooks https://orcid.org/0000-0002-8458-5047

D. L. Burke https://orcid.org/0000-0003-1866-1950

A. Carnero Rosell https://orcid.org/0000-0003-3044-5150

M. Carrasco Kind https://orcid.org/0000-0002-4802-3194

J. Carretero https://orcid.org/0000-0002-3130-0204

L. N. da Costa https://orcid.org/0000-0002-7731-277X

M. E. S. Pereira https://orcid.org/0000-0002-7131-7684

T. M. Davis https://orcid.org/0000-0002-4213-8783

S. Desai https://orcid.org/0000-0002-0466-3288

H. T. Diehl https://orcid.org/0000-0002-8357-7467

S. Everett https://orcid.org/0000-0002-3745-2882

I. Ferrero https://orcid.org/0000-0002-1295-1132

D. Friedel https://orcid.org/0000-0002-3632-7668

J. García-Bellido https://orcid.org/0000-0002-9370-8360

S. R. Hinton https://orcid.org/0000-0003-2071-9349

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

D. J. James https://orcid.org/0000-0001-5160-4486

K. Kuehn https://orcid.org/0000-0003-0120-0808

J. Mena-Fernández https://orcid.org/0000-0001-9497-7266

F. Menanteau https://orcid.org/0000-0002-1372-2534

R. Miquel https://orcid.org/0000-0002-6610-4836

R. L. C. Ogando https://orcid.org/0000-0003-2120-1154

A. Pieres https://orcid.org/0000-0001-9186-6042

A. A. Plazas Malagón https://orcid.org/0000-0002-2598-0514

E. Sanchez https://orcid.org/0000-0002-9646-8198

M. Smith https://orcid.org/0000-0002-3321-1432

E. Suchyta https://orcid.org/0000-0002-7047-9358

M. E. C. Swanson https://orcid.org/0000-0002-1488-8552

G. Tarle https://orcid.org/0000-0003-1704-0781

C. To https://orcid.org/0000-0001-7836-2261

A. R. Walker https://orcid.org/0000-0002-7123-8943

Y. Zhang https://orcid.org/0000-0001-5969-4631

## References

Abbott, T. M. C., Adamów, M., Aguena, M., et al. 2021, ApJS, 255, 20

Alexandersen, M., Benecchi, S. D., Chen, Y.-T., et al. 2019, ApJS, 244, 19

Astropy Collaboration 2013, A&A, 558, A33

Astropy Collaboration 2018, AJ, 156, 123

Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18

Benecchi, S. D., Noll, K. S., Grundy, W. M., et al. 2009, Icar, 200, 292

Benecchi, S. D., & Sheppard, S. S. 2013, AJ, 145, 124

Bernardinelli, P. H., Bernstein, G. M., & Jindal, N., 2023 Photometry of outer solar system objects from the Dark Energy Survey I: photometric methods, light curve distributions and trans-Neptunian binaries—scene modeling photometry, v1.0.0, Zenodo:10.5281/zenodo.8231232

Bernardinelli, P. H., Bernstein, G. M., Montet, B. T., et al. 2021, ApJL, 921, L37

Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2020, ApJS, 247, 32

Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2022, ApJS, 258, 41

Bernstein, G. M., Abbott, T. M., Armstrong, R., et al. 2018, PASP, 130, 054501

Bernstein, G. M., Armstrong, R., Plazas, A. A., et al. 2017, PASP, 129, 074503

Bernstein, G. M., Holler, B. J., Navarro-Escamilla, R., et al. 2023, PSJ, 4, 116

Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364

Brout, D., Sako, M., Scolnic, D., et al. 2019, ApJ, 874, 106

Brown, M. E., Barkume, K. M., Ragozzine, D., & Schaller, E. L. 2007, Natur, 446, 294

Burke, D. L., Rykoff, E. S., Allam, S., et al. 2017, AJ, 155, 41

Chen, Y.-T., Eduardo, M. R., Muñoz-Gutiérrez, M. A., et al. 2022, ApJL, 937, L22

Doressoundiram, A., Boehnhardt, H., Tegler, S. C., & Trujillo, C. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 91

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fraser, W., Alexandersen, M., Schwamb, M. E., et al. 2016, AJ, 151, 158

Fraser, W. C., & Brown, M. E. 2012, ApJ, 749, 33

Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100

Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, in Nomenclature in the Outer Solar System, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 43

Gladman, B., & Volk, K. 2021, ARA&A, 59, 203

Grundy, W., Noll, K., Nimmo, F., et al. 2011, Icar, 213, 678

Huang, Y., Gladman, B., Volk, K., et al. 2022, ApJS, 259, 54

Hunter, J. D. 2007, CSE, 9, 90

Jarvis, M., Bernstein, G. M., Amon, A., et al. 2020, MNRAS, 501, 1282

Jeffreys, H. 1961, Theory of Probability (3rd ed.; Oxford: Oxford Univ. Press)

Jones, E., Oliphant, T., Peterson, P., et al., 2001 SciPy: Open source scientific tools for Python, http://www.scipy.org/

Kavelaars, J. J., Petit, J.-M., Gladman, B., et al. 2021, ApJL, 920, L28

Kecskeméthy, V., Kiss, C., Szakáts, R., et al. 2023, ApJS, 264, 18

Khain, T., Becker, J. C., Lin, H. W., et al. 2020, AJ, 159, 133

Lacerda, P., & Luu, J. 2006, AJ, 131, 2314

Lacerda, P., McNeill, A., & Peixinho, N. 2014, MNRAS, 437, 3824

Nesvorný, D., & Vokrouhlický, D. 2019, Icar, 331, 49

Nesvorny, D. 2018, ARA&A, 56, 137

Noll, K. 2007, HST Proposal, 11113

Noll, K., Grundy, W. M., Nesvorný, D., & Thirouin, A. 2020, in The Trans-Neptunian Solar System, ed. D. Prialnik, M. A. Barucci, & L. Young (Amsterdam: Elsevier), 201

Oliphant, T. E. 2015, Guide to NumPy (2nd ed.; USA: CreateSpace Independent Publishing Platform)

Parker, A. H., Kavelaars, J. J., Petit, J.-M., et al. 2011, ApJ, 743, 1

Pérez, F., & Granger, B. E. 2007, CSE, 9, 21

Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131

Sánchez, E. 2006, in AIP Conf. Proc. 878, 2nd Int. Conf. on The Dark Side of the Universe, ed. C. Muñoz & G. Yepes (Melville, NY: AIP), 213

Schemel, M., & Brown, M. E. 2021, PSJ, 2, 40

Schwamb, M. E., Fraser, W. C., Bannister, M. T., et al. 2019, ApJS, 243, 12

Sevilla-Noarbe, I., Bechtol, K., Kind, M. C., et al. 2021, ApJS, 254, 24

Showalter, M. R., Benecchi, S. D., Buie, M. W., et al. 2021, Icar, 356, 114098
Snodgrass, C., Carry, B., Dumas, C., & Hainaut, O. 2010, A&A, 511, A72

Stephens, D. C., & Noll, K. S. 2006, AJ, 131, 1142

Strauss, R., Trilling, D. E., Bernardinelli, P. H., et al. 2023, arXiv:2309.04034

The Dark Energy Survey Collaboration 2016, MNRAS, 460, 1270

The Dark Energy Survey Collaboration 2022, PhRvD, 105, 023520

Thirouin, A., & Sheppard, S. S. 2019, AJ, 157, 228

Trilling, D., Gerdes, D. W., Juric, M., et al. 2023, arXiv:2309.03417

Vilenius, E., Kiss, C., Mommert, M., et al. 2012, A&A, 541, A94