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

# Milky Way Satellite Census. I. The Observational Selection Function for Milky Way Satellites in DES Y3 and Pan-STARRS DR1

A. Drlica-Wagner¹,²,³, K. Bechtol⁴,⁵, S. Mau²,³, M. McNanna⁴, E. O. Nadler⁶,⁷, A. B. Pace⁸,⁹, T. S. Li¹,²,¹⁰,¹¹,<fnref n="42" />, A. Pieres¹²,¹³, E. Rozo¹⁴, J. D. Simon¹⁰, A. R. Walker¹⁵, R. H. Wechsler⁶,⁷,¹⁶, T. M. C. Abbott¹⁵, S. Allam¹, J. Annis¹, E. Bertin¹⁷,¹⁸, D. Brooks¹⁹, D. L. Burke⁷,¹⁶, A. Carnero Rosell¹²,²⁰, M. Carrasco Kind²¹,²², J. Carretero²³, M. Costanzi²⁴,²⁵, L. N. da Costa¹²,¹³, J. De Vicente²⁰, S. Desai²⁶, H. T. Diehl¹, P. Doel¹⁹, T. F. Eifler²⁷,²⁸, S. Everett²⁹, B. Flaugher¹, J. Frieman¹,², J. García-Bellido³⁰, E. Gaztanaga³¹,³², D. Gruen⁶,⁷,¹⁶, R. A. Gruendl²¹,²², J. Gschwend¹²,¹³, G. Gutierrez¹, K. Honscheid³³,³⁴, D. J. James³⁵, E. Krause²⁷, K. Kuehn³⁶,³⁷, N. Kuropatkin¹, O. Lahav¹⁹, M. A. G. Maia¹²,¹³, J. L. Marshall⁹, P. Melchior¹¹, F. Menanteau²¹,²², R. Miquel²³,³⁸, A. Palmese¹,², A. A. Plazas¹¹, E. Sanchez²⁰, V. Scarpine¹, M. Schubnell³⁹, S. Serrano³¹,³², I. Sevilla-Noarbe²⁰, M. Smith⁴⁰, E. Suchyta⁴¹, and G. Tarle³⁹

(DES Collaboration)

¹ Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA; kadrlica@fnal.gov

² Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

³ Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

⁴ Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue Madison, WI 53706-1390, USA; bechtol@wisc.edu

⁵ LSST, 933 North Cherry Avenue, Tucson, AZ 85721, USA

⁶ Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

⁷ Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA

⁸ Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15312, USA

⁹ 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

¹⁰ Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA

¹¹ Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

¹² Laboratório Interinstitucional de e-Astronomia—LineA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

¹³ Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

¹⁴ Department of Physics, University of Arizona, Tucson, AZ 85721, USA

¹⁵ Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

¹⁶ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

¹⁷ CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

¹⁸ Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

¹⁹ Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

²⁰ Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

²¹ National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA

²² Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

²³ Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain

²⁴ INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

²⁵ Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

²⁶ Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

²⁷ Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

²⁸ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

²⁹ Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

³⁰ Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

³¹ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

³² Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

³³ Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

³⁴ Department of Physics, The Ohio State University, Columbus, OH 43210, USA

³⁵ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

³⁶ Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

³⁷ Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

³⁸ Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

³⁹ Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁴⁰ School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

⁴¹ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

Received 2019 December 4; revised 2020 March 1; accepted 2020 March 2; published 2020 April 15

## Abstract

We report the results of a systematic search for ultra-faint Milky Way satellite galaxies using data from the Dark Energy Survey (DES) and Pan-STARRS1 (PS1). Together, DES and PS1 provide multi-band photometry in optical/near-infrared wavelengths over ~80% of the sky. Our search for satellite galaxies targets ~25,000 deg² of the high-Galactic-latitude sky reaching a 10σ point-source depth of ≳22.5 mag in the g and r bands. While satellite galaxy searches have been performed independently on DES and PS1 before, this is the first time that a self-

<footnote n="42">⁴² NHFP Einstein Fellow.</footnote>
consistent search is performed across both data sets. We do not detect any new high-significance satellite galaxy candidates, recovering the majority of satellites previously detected in surveys of comparable depth. We characterize the sensitivity of our search using a large set of simulated satellites injected into the survey data. We use these simulations to derive both analytic and machine-learning models that accurately predict the detectability of Milky Way satellites as a function of their distance, size, luminosity, and location on the sky. To demonstrate the utility of this observational selection function, we calculate the luminosity function of Milky Way satellite galaxies, assuming that the known population of satellite galaxies is representative of the underlying distribution. We provide access to our observational selection function to facilitate comparisons with cosmological models of galaxy formation and evolution.

*Unified Astronomy Thesaurus concepts: Dwarf spheroidal galaxies (420); Local Group (929); Dark matter (353); Milky Way Galaxy (1054)*

*Supporting material: machine-readable tables*

# 1. Introduction

Faint dwarf galaxies dominate the universe by number, yet a precise census of these objects remains challenging, due to the limited sensitivity of observational surveys. Dwarf galaxies with stellar mass $\lesssim 10^6 M_\odot$ have only been identified within the Local Volume (distances of a few Mpc), either in isolation or as satellites of larger galaxies (e.g., Martin et al. 2013; Müller et al. 2015; Carlin et al. 2016; Smercina et al. 2018; Crnojević et al. 2019). At even lower masses, the census of ultra-faint satellites is incomplete, even within the Milky Way halo. Despite significant observational challenges, the demographics of ultra-faint dwarf galaxies offer a unique window into feedback processes in galaxy formation (e.g., Mashchenko et al. 2008; Wheeler et al. 2015, 2019; Agertz et al. 2020; Munshi et al. 2019), reionization and the first stars (e.g., Bullock et al. 2000; Shapiro et al. 2004; Weisz et al. 2014a, 2014b; Boylan-Kolchin et al. 2015; Ishiyama et al. 2016; Weisz & Boylan-Kolchin 2017; Tollerud & Peek 2018; Graus et al. 2019; Katz et al. 2019), and the nature of dark matter (e.g., Bergström et al. 1998; Spekkens et al. 2013; Malyshev et al. 2014; Ackermann et al. 2015; Geringer-Sameth et al. 2015; Bullock & Boylan-Kolchin 2017; Clesse & García-Bellido 2018; Nadler et al. 2019a).

The lowest-luminosity satellite galaxies are detected in optical imaging surveys as arcminute-scale statistical overdensities of individually resolved stars. Beginning with the Sloan Digital Sky Survey (SDSS), digitized wide-area multi-band optical imaging surveys—combined with automated search algorithms—have greatly increased the known population of Milky Way satellites (Willman et al. 2005a, 2005b; Belokurov et al. 2006, 2007, 2008, 2009, 2010; Grillmair 2006, 2009; Sakamoto & Hasegawa 2006; Zucker et al. 2006a, 2006b; Irwin et al. 2007; Walsh et al. 2007; Kim et al. 2015a). More recently, searches using data from the Dark Energy Survey (DES; Bechtol et al. 2015; Drlica-Wagner et al. 2015; Koposov et al. 2015; Kim & Jerjen 2015b; Luque et al. 2016), other DECam surveys (e.g., SMASH, MagLiteS, and DELVE; Martin et al. 2015; Drlica-Wagner et al. 2016; Koposov et al. 2018; Torrealba et al. 2018; Mau et al. 2020), ATLAS (Torrealba et al. 2016a, 2016b), Pan-STARRS1 (PS1; Laevens et al. 2015a, 2015b), and Gaia (Torrealba et al. 2019b) have further increased the sample of confirmed and candidate satellites to more than 50 (Figure 1).

Both observational and theoretical arguments suggest that the current census of Milky Way satellite galaxies is incomplete. From an observational standpoint, this incompleteness is demonstrated by the continued discovery of fainter, more distant, and lower surface brightness systems. For example, the first ~700 deg$^2$ of deep imaging with the Hyper Suprime-Cam Strategic Survey

Program (HSC SSP) has revealed three new satellites at sufficiently low luminosities and large heliocentric distances that escaped detection by earlier overlapping surveys (Homma et al. 2016, 2018, 2019). Moreover, several recently discovered Milky Way companions (e.g., Crater II, Virgo I, Aquarius II, Cetus III, Antlia II, and Boötes IV) are lower surface brightness than most ultra-faint dwarfs discovered in the SDSS era, implying that the current generation of surveys and search techniques are sensitive to systems that were previously undetectable. Searches using compact spatial kernels and a wider variety of stellar population ages and metallicities have revealed diverse Milky Way substructures (Torrealba et al. 2019a), and precise proper motion information for billions of nearby stars provided by Gaia has enlarged the sample of extremely low-surface-brightness satellites (Torrealba et al. 2019b).

Theoretical predictions for the smallest galaxies have advanced hand in hand with observations. Since galaxy formation is a nonlinear process, numerical simulations have long been used to predict the population statistics of these objects. Early simulations that resolved dark matter substructure within Milky Way-mass halos predicted far more surviving dark matter subhalos than the number of observed satellites (Klypin et al. 1999; Moore et al. 1999). This mismatch, dubbed the “missing satellites problem,” simply reflects the fact that mapping subhalos in dark-matter-only simulations to observed satellites is nontrivial. In particular, reionization and stellar feedback drastically suppress dwarf galaxy formation in low-mass halos (e.g., Bullock et al. 2000; Somerville 2002; Brown et al. 2014), and tidal interactions with the Galactic disk are expected to disrupt a significant number of systems (e.g., Garrison-Kimmel et al. 2017; Nadler et al. 2018; Kelley et al. 2019). Semi-empirical models that account for these effects— along with realistic satellite detection criteria—find that the observed satellite population is consistent with cold, collisionless dark matter (e.g., Jethwa et al. 2018; Kim et al. 2018; Newton et al. 2018; Bose et al. 2019; Nadler et al. 2019a, 2019b). Likewise, hydrodynamic simulations that self-consistently model galaxy formation in a cosmological context produce luminosity functions and radial distributions of satellites that are broadly consistent with observations of the Milky Way system (e.g., Wetzel et al. 2016; Garrison-Kimmel et al. 2019; Samuel et al. 2020). In concert, extremely high-resolution simulations of isolated ultra-faint systems suggest that low-mass dwarfs may be abundant (Wheeler et al. 2019).

Historically, the primary means of comparing Milky Way satellite observations to simulations has been through the total satellite luminosity function (i.e., the total number of satellites within the virial radius of the Milky Way halo as a function of satellite luminosity). Typically, an observational selection
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Equal-area Mollweide projection of the density of stars (log scale) observed with *r* < 22 mag over the DES Y3A2 and PS1 DR1 footprints (red and cyan borders, respectively). Filled markers indicate kinematically confirmed Milky Way satellite galaxies, and unfilled markers indicate satellite galaxy candidates that have not yet been kinematically confirmed. We mark classical Milky Way satellites (black circles), and satellites discovered by SDSS (blue squares), PS1 (cyan diamonds), DES (red upright triangles), other DECam surveys (purple inverted triangles), HSC SSP (green pluses), VLT ATLAS (magenta pentagons), and Gaia (brown crosses).</figcaption></figure>

function is built to predict the detectability of a satellite as a function of heliocentric distance, size, and luminosity. This type of analysis was pioneered by Koposov et al. (2008) and Walsh et al. (2009), who used simulations to characterize the satellite detection efficiency in SDSS, analyzing ~8000 deg² from SDSS DR5 and ~9500 deg² from SDSS DR6, respectively. The total luminosity function was derived by correcting the observed satellite population for observational selection effects, and the result was compared to cosmological predictions. Recently, several studies have begun to utilize more advanced model inference techniques that require a simple yet comprehensive mechanism to predict the detectability of a satellite (Jethwa et al. 2018; Newton et al. 2018; Nadler et al. 2019b). However, these studies have been limited by the lack of rigorous estimates for the selection functions of modern surveys.

In this paper, we present a systematic search for Milky Way satellites and a detailed quantitative measurement of the observational selection function for modern surveys. In particular, we performed an updated search for Milky Way satellites by applying two independent search algorithms to ~5000 deg² of data from DES Collaboration et al. (2018) and ~30,000 deg² of data from PS1 (Chambers et al. 2016). After quality cuts, our analysis covers approximately three times the sky area analyzed by Koposov et al. (2008) and Walsh et al. (2009). The DES data were collected during the first three years of survey operations and cover much of the southern Galactic cap (DES Y3A2; DES Collaboration et al. 2018; Shipp et al. 2018). When compared to previous DES satellite searches (i.e., Drlica-Wagner et al. 2015), DES Y3A2 has ~50% more exposure time, more homogeneous coverage, more accurate photometric calibration, and more efficient star-galaxy classification (e.g., Burke et al. 2018; DES Collaboration et al. 2018). To extend the coverage of our analysis to the Northern Hemisphere, we also apply our search algorithms to publicly available data from the first data release of PS1 (PS1 DR1; Chambers et al. 2016). Note that in most regions of the sky at high-Galactic latitude, the number density of background galaxies exceeds that of foreground Milky Way stars at magnitudes *r* ≥ 22. Accordingly, this analysis represents a systematic search over ~75% of the high-Galactic-latitude sky, reaching depths at which the stellar sample is limited primarily

by star-galaxy confusion rather than object detection (e.g., Fadely et al. 2012).

We quantify the observational selection function of our search to facilitate direct comparisons between the observed luminosity function and predictions from simulations. We simulate the resolved stellar populations of 10⁵ (10⁶) satellites and inject simulated stars into the DES Y3A2 (PS1 DR1) data at the catalog level. These simulations span a range of absolute magnitudes, heliocentric distances, physical sizes, ellipticities, position angles, ages, and metallicities. We run our search algorithms on each simulated satellite and find that the detectability of a satellite can be well described by its absolute magnitude, heliocentric distance, physical size, and local stellar density. We derive both analytic and machine-learning models that predict the detectability of a satellite as a function of these parameters.

The observational selection function derived in this paper can be used to test models that predict the abundance and properties of Milky Way satellites. As an illustrative example, we use our observational selection function to derive the total luminosity function of Milky Way satellites solely based on the properties of the observed population. In a companion paper (Nadler et al. 2019c, hereafter Paper II), we use high-resolution numerical simulations (including a model for the effects of baryons) to build a more rigorous model of the observed satellite population and to constrain models of galaxy formation. In deriving the observational selection function, we have intentionally set a high threshold for detection in order to provide a clean interpretation of the resulting satellite populations. The investigation of lower significance candidates is left to future work.

This paper is structured as follows. In Section 2, we provide a high-level overview of our simulation and analysis pipeline. The subsequent sections provide more detail on the survey data sets (Section 3), our catalog-level simulations (Section 4), and the satellite search algorithms (Section 5). In Section 6, we present the results of our search on the DES and PS1 data. The resulting observational selection functions derived from simulations are presented in Section 7, and our simple luminosity function inference is presented in Section 8. We conclude in Section 9. Our models for the observational selection functions of DES and PS1 are publicly available online.<fnref n="43" />

<footnote n="43">43 https://github.com/des-science/mw-sats</footnote>
## 2. Analysis Overview

In this section, we summarize the key components of the simulation and data analysis pipeline used to derive the observational selection function for Milky Way satellites. We applied two distinct algorithms to search for satellite galaxies in photometric catalog data from DES and PS1. To evaluate the sensitivity of our search, we embedded simulated satellite galaxies into these data and attempted to recover them with the same search algorithms. By self-consistently analyzing the data and simulations, we accurately characterized both the population of observed satellites and the population of satellites that remain undetected due to the limited sensitivity of our observations.

To generate realistic satellite galaxy simulations, we empirically modeled the survey coverage and photometric response of DES Y3A2 and PS1 DR1. We characterized the coverage, depth, completeness, and photometric measurement uncertainties as a function of sky location for each survey. We then simulated stellar catalogs for satellites covering a large range of physical properties, including sky location, luminosity, heliocentric distance, physical size, ellipticity, age, and metallicity. For each satellite, we generated a Poisson realization of the observable stellar distribution, simulating the position, flux, and photometric uncertainty of each star. These simulated satellites were injected one-by-one into the survey data sets, and two search algorithms were run at the location of each injected system. The analysis of these simulations produces a multi-dimensional vector containing the detection significance of each simulated satellite as a function of its intrinsic properties (e.g., luminosity, distance, and physical size) and global survey properties (e.g., survey depth, coverage, and local foreground stellar density). We refer to the mapping between satellite properties and satellite detectability as the observational selection function.

In parallel, we performed an untargeted search of DES Y3A2 and PS1 DR1 without any embedded simulations. This search produced a set of stellar overdensity “seeds” of varying significance. The most significant seeds are associated with physical systems reported in the literature, while less significant seeds can be attributed to statistical fluctuations, artificial variations in the stellar density due to survey systematics, and sub-threshold physical systems. We characterized the distribution of detection significances for the collection of seeds and defined a conservative detection threshold that recovers a large fraction of systems that were discovered in surveys of comparable depth. This a posteriori definition of a detection threshold was required to deal with systematic artifacts that contaminated the population of seeds and made it impossible to choose a statistical threshold a priori. However, by self-consistently applying the same detection threshold to the population of simulated satellites, we can determine the selection efficiency for any detection threshold and satellite properties. In this paper, we are primarily concerned with the demographics of the satellite galaxy population rather than detecting new, low-significance candidates, and therefore we set our significance threshold to yield a pure sample of “detected” satellites.

The detected population of satellites and the observational selection function can be combined to derive the Milky Way satellite galaxy luminosity function. Analyzing the results of the simulations directly can be cumbersome and computationally intensive, while representing the observational selection function with a simple analytic relationships discards some information. Therefore, to simplify the application of the selection function while retaining detailed information, we trained a gradient-

boosted decision tree classifier that takes as input characteristics of a satellite (e.g., size, luminosity, distance, and local stellar density) and outputs a probability that the satellite would be detected. When applying this classifier, we combine the satellite properties with the global geometric characteristics of each survey in the form of HEALPix maps of survey coverage.

## 3. Data Set

Data from DES Y3A2 and PS1 DR1 cover ~5000 deg² and ~30,000 deg² of the celestial sphere, respectively (Figure 2). The deep, multi-band, optical/near-infrared imaging of these surveys provides the photometric, astrometric, and morphological measurements necessary to separate stellar overdensities in the Milky Way halo from Milky Way field stars and unresolved background galaxies. In this section, we describe the selection of high-quality stellar samples for each of these surveys, the characterization of the survey geometry, and determination of survey response as a function of location on the sky. Additional technical details on our selections are provided in Appendix A.

### 3.1. DES Y3A2

DES is a broadband optical/near-infrared imaging survey of the southern Galactic cap using the Dark Energy Camera (DECam; Flaugher et al. 2015) mounted at the prime focus of the 4 m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO). Here, we analyze data from the first three years of DES operations (Diehl et al. 2016). The DES Y3A2 imaging data serves as the basis for the first DES public data release (DES DR1; DES Collaboration et al. 2018) and consists of ~45,000 wide-area survey exposures. Details of the DES image reduction and catalog generation can be found in Morganson et al. (2018), while more details on the internal DES Y3A2 data set can be found in I. Sevilla-Noarbe et al. (2020, in preparation).

**Photometry.** The internal DES Y3A2 object catalogs augment DES DR1 with additional multi-band, multi-epoch, forced photometry, which provides significantly improved photometric and morphological measurements of faint objects. These catalogs were generated in two steps. First, individual sources were detected in *riz* coadded images using Source-Extractor (Bertin & Arnouts 1996), with a detection threshold of S/N ~ 10 (Morganson et al. 2018). This coadd object catalog was then used as input to the *ngmix* multi-band, multi-epoch fitting routine, which performs a simultaneous fit of source parameters across the set of individual *griz* images for each object (Sheldon 2014; Drlica-Wagner et al. 2018). *ngmix* is run in two configurations: (1) fits are performed on single objects while masking nearby neighbors, referred to as the “single object fit” (SOF), and (2) fits are performed iteratively on groups of objects, referred to as the “multi-object fit” (MOF). The treatment of the PSF on an image-by-image basis substantially improves point-source photometry and star-galaxy separation (Sevilla-Noarbe et al. 2018). The relative top-of-the-atmosphere photometric accuracy across the DES footprint is estimated to be <7 mmag from a comparison to Gaia (Burke et al. 2018). Furthermore, DES Y3A2 includes SED-dependent chromatic corrections for each object based on an initial evaluation of the stellar spectral type (Li et al. 2016; I. Sevilla-Noarbe et al. 2020, in preparation).

**Coverage.** The observational coverage of DES Y3A2 was assembled in a vectorized format using *mangle* (Hamilton & Tegmark 2004; Swanson et al. 2008). The *mangle* representation
<figure><img src="image_3.png" /><figcaption>**Figure 2.** McBryde-Thomas flat polar quartic projection of the geometric masks applied to the DES and PS1 satellite searches. Colored regions are masked from our analysis because they lie outside the DES or PS1 survey footprint (dark blue), have interstellar extinction $E(B-V) > 0.2$ mag in the Schlegel et al. (1998) dust map (blue), or are associated with previously known stellar overdensities, galaxy clusters, or bright stars (green). The total unmasked area is 4844 deg$^2$ for DES, 21,123 deg$^2$ for PS1, and 24,343 deg$^2$ for the combined search.</figcaption></figure>

accounts for missing coverage at the boundary of the survey footprint, as well as gaps associated with saturated stars, bleed trails, and other instrument signatures. These vectorized maps were then converted into a subsampled *nside* = 4096 HEALPix map with a resolution of ~0.74 arcmin$^2$ per pixel (Drlica-Wagner et al. 2018). We restrict the survey footprint to HEALPix pixels where the *griz* sky coverage fraction is greater than 0.5, resulting in a total effective solid angle before masking of 4945 deg$^2$.

**Depth.** The depth of DES Y3A2 was estimated for each band in each *nside* = 4096 HEALPix pixel. This involved combining the *mangle* maps with additional survey characteristics following the procedure developed in Rykoff et al. (2015), as described in Section 7.1 of Drlica-Wagner et al. (2018). In brief, we trained a random forest classifier that combined survey characteristics, such as coverage, seeing, and sky brightness, to estimate the 10σ limiting magnitude in each pixel. The 10σ SOF CM\_MAG depth for DES Y3A2 data in regions with $E(B-V) < 0.2$ is $g = 23.9$ and $r = 23.7$. When simulating and analyzing satellites in DES Y3A2, we incorporated the (small) variations in depth over the footprint.

**Reddening.** We followed the procedure described in Section 4.2 of DES Collaboration et al. (2018) to correct for interstellar extinction. We applied an additive correction to each measured magnitude of $A_b = R_b \times E(B-V)$, where $E(B-V)$ comes

from Schlegel et al. (1998, SFD). We use the $R_b$ coefficients from DES Collaboration et al. (2018)—specifically $R_g = 3.186$ and $R_r = 2.140$. These fiducial coefficients are derived using the Fitzpatrick (1999) reddening law with $R_V = 3.1$ and incorporate the renormalization of the SFD reddening map ($N = 0.78$) suggested by Schlafly et al. (2010). Hereafter, all DES Y3A2 magnitudes are extinction corrected.

**Star–Galaxy Separation.** The star–galaxy separation efficiency of previous DES data sets was estimated to be >95% complete for stars at $i < 22$ (Drlica-Wagner et al. 2018). Sevilla-Noarbe et al. (2018) showed that the efficiency of star–galaxy separation could be considerably improved by using multi-epoch morphological measurements. We assessed the efficiency of star–galaxy classification in DES Y3A2 by comparing to overlapping deep data from the HSC SSP (Aihara et al. 2018). Accounting for both object detection and classification, we found that the DES Y3A2 SOF classifier$^{44}$ achieves >90% completeness for stars down to a magnitude limit of $r \sim 23.5$ mag. The stellar completeness of ground-based surveys such as DES and PS1 is largely dominated by

<footnote n="44">44 We use $0 \le \text{EXTENDED\_CLASS\_MASH\_SOF} \le 2$ to define the stellar sample, balancing stellar completeness and galaxy contamination. Details on the performance of star–galaxy classification for DES Y3A2 will be presented in I. Sevilla-Noarbe et al. (2020, in preparation).</footnote>
the efficiency of star–galaxy classification rather than the point-source detection limit.

## 3.2. PS1 DR1

Our PS1 data set consists of data from the first public data release of the PS1 3π Survey (Chambers et al. 2016). PS1 DR1 is assembled from images taken by the PS1 Gigapixel Camera #1 (Tonry et al. 2008) mounted on the 1.8 m PS1 telescope at Haleakala Observatories on the island of Maui, Hawaii. PS1 DR1 covers the northern and equatorial sky with $\delta_{2000} > -30$ deg in five optical bands, $g_{P1}$, $r_{P1}$, $i_{P1}$, $z_{P1}$, $y_{P1}$ (Tonry et al. 2012). The DES and PS1 g, r, and i bands are similar enough that we drop the “P1” subscript, though we analyzed each survey in its respective filter system.

**Photometry.** We selected PS1 DR1 objects from For-cedMeanObject-View with qualityFlag & 16 > 0, nDetections > 0, and nStackDetections > 1. We removed duplicate objects from the catalogs by selecting only objects that are the primary detection, detectPrimary = 1. We made several additional quality cuts based on the InfoFlag and InfoFlag2 variables to remove objects where the photometric fit failed, objects that were likely to be defects, objects with too few points measured to derive an elliptical contour, and objects where all model fits failed. The full set of selection criteria can be found in Appendix A. These criteria were validated by comparing against catalogs derived from the HSC SXDS ultra-deep field. We found that our cuts retain the majority of PS1-HSC matched objects while significantly reducing the incidence of spurious objects in the PS1 data. We converted the measured fluxes from PS1 DR1 into magnitudes by applying a stack zero-point of 8.9 (i.e., mag = -2.5 log₁₀(flux) + 8.9). When PS1 DR1 reports nega-tive flux values, we set the corresponding magnitude and magnitude uncertainty to a sentinel value of -999. Since we are using (g - r) colors for our search, we selected only objects with measured PSF magnitudes in both the g and r bands.

**Coverage.** We approximated the coverage of PS1 empiri-cally using the full PS1 DR1 catalog prior to any star–galaxy separation or photometric cuts. We define the PS1 footprint as the set of nside = 2048 (~2.95 arcmin²) HEALPix pixels that contain any PS1 object.<fnref n="45" /> This coverage map is not strictly accurate, since some HEALPix pixels that contain objects are not fully covered by the survey while other HEALPix pixels are covered by the survey yet contain no objects. However, at the level of accuracy necessary for our search algorithms, this coverage map is sufficient to avoid significantly biasing estimates of the local stellar density. We converted from nside = 2048 to nside = 4096 (i.e., setting the value of each nested subpixel to the value of its parent) for use by the analysis algorithms. The total area of the PS1 DR1 footprint that we consider before masking is 29,343 deg².

**Depth.** We estimated the photometric depth of PS1 DR1 by interpolating the median magnitude uncertainty as a function of magnitude for a set of low-reddening, high-Galactic-latitude regions. We determined the magnitude at which the median magnitude uncertainty is 0.1085, corresponding to the 10σ detection limit. We found typical 10σ magnitude limits for PS1 DR1 of g = 22.5 and r = 22.4. These depth estimates agree

with those of Chambers et al. (2016), who estimate that the PS1 DR1 catalog retains 98% completeness at g, r, i ~ 22.5, with a spatial variation of ±0.25 mag (see figure 17 of Chambers et al. 2016). We assumed a constant magnitude limit for PS1 DR1; however, variable interstellar extinction introduces a spatially dependent intrinsic magnitude limit for stars (i.e., we are less sensitive in regions of high extinction). We have found that stars fainter than the 10σ magnitude limit can contribute significantly to the detectability of faint satellites. For this reason, we used a magnitude limit of g, r = 23 when performing the likelihood-based search (Section 5.2). In the spatial matched-filter analysis (Section 5.1), which uses a binary selection in color–magnitude space rather than a likelihood-based weighting of photometric uncertainties, we applied a signal-to-noise threshold that limits the completeness of faint stars but better controls the number of spurious seeds returned by the algorithm.

**Reddening.** We corrected the PS1 DR1 measured magni-tudes for interstellar extinction following the procedure described in Schlafly & Finkbeiner (2011). E(B – V) values were calculated from each object by performing a bilinear interpolation to the maps of Schlegel et al. (1998). We then applied $R_b$ coefficients from table 6 of Schlafly & Finkbeiner (2011), assuming the reddening law of Fitzpatrick (1999) with $R_V$ = 3.1. For PS1 DR1, this corresponds to $R_g$ = 3.172 and $R_r$ = 2.271 for the g and r bands, respectively. We note that these $R_b$ values include a renormalization factor, as suggested by Schlafly et al. (2010). Hereafter, all PS1 DR1 magnitudes are extinction corrected unless explicitly stated otherwise.

**Star–Galaxy Separation.** We performed star–galaxy separa-tion by comparing the measured i-band PSF and adaptive aperture magnitudes, iFPSFMag and iFKronMag, respec-tively. The choice of i-band was motivated by the superior PSF and depth in this band. Our primary cut required that the measured PSF and aperture magnitudes agree, iFPSFMag – iFKronMag < 0.05. However, we found that the PS1 PSF fit often fails in dense stellar regions. To retain sensitivity in these regions, we also included objects where iFPSFMag = -999 or iFPSFMag – iFKronMag > 4. By comparing to HSC SSP (Aihara et al. 2018), we find that our PS1 DR1 stellar sample is >90% complete down to a magnitude of r ~ 21.7.

## 4. Satellite Simulations

We simulated Milky Way satellite galaxies with a wide range of properties to accurately quantify the detection efficiency of our algorithms. We randomly sampled values of stellar luminosity, heliocentric distance, physical size, ellipti-city, and position angle from the ranges described in Table 1. Stellar photometry was simulated based on stellar isochrones from Bressan et al. (2012), selected from a range of ages and metallicities characteristic of observed satellites (Table 1). Photometric error models were derived for each survey and were used to assign a photometric uncertainty to each star and randomize the measured photometry relative to the determi-nistic photometry provided from the isochrone. The simulated satellites were randomly assigned spatial locations in a region that slightly overcovered each of the survey footprints. The population of simulated satellites was not intended to mimic any realistic satellite population; rather, it was intended to cover the range of parameter space where variations in detection efficiency occur.

<footnote n="45">⁴⁵ The choice of nside = 2048 was determined empirically, as the highest resolution that did not contain a significant number of pixels that lacked objects but were clearly covered by PS1 (as determined by visual inspection of the coadd images).</footnote>
<table><caption>**Table 1**
Parameter Ranges of Simulated Satellites</caption><thead><tr><th>Parameter</th><th>Range</th><th>Unit</th><th>Sampling</th></tr></thead><tbody><tr><td>Stellar Mass</td><td>[10, 10<sup>6</sup>]</td><td>M<sub>☉</sub></td><td>log</td></tr><tr><td>Heliocentric Distance</td><td>[5, 10<sup>3</sup>]</td><td>kpc</td><td>log</td></tr><tr><td>2D Half-light Radius</td><td>[1, 2 × 10<sup>3</sup>]</td><td>pc</td><td>log</td></tr><tr><td>Ellipticity</td><td>[0.1, 0.8]</td><td>...</td><td>linear</td></tr><tr><td>Position Angle</td><td>[0, 180]</td><td>deg</td><td>linear</td></tr><tr><td>Age</td><td>{10, 12, 13.5}</td><td>Gyr</td><td>choice</td></tr><tr><td>Metallicity</td><td>{0.0001, 0.0002}</td><td>...</td><td>choice</td></tr></tbody></table>

Note. Simulated satellite properties are drawn from uniform distributions in log space with the parameter ranges listed above. Isochrones were generated using the models of Bressan et al. (2012) with an IMF from Chabrier (2001).

Satellites were simulated at the catalog level as collections of individually resolved stars. To generate realistic catalogs, we began with a probabilistic model for the spatial and flux distributions of stars in each satellite. We sampled the spatial distribution of stars according to a Plummer profile (Plummer 1911), which has been found to be a good description of known Milky Way satellite galaxies (e.g., Simon 2019).<fnref n="46" /> We use $a_h$ to indicate the elliptical semimajor axis containing half the light (arcmin) and $r_h = a_h\sqrt{1 - e}$ to represent the azimuthally averaged half-light radius (arcmin), where $e$ is the ellipticity. $a_{1/2}$ and $r_{1/2}$ represent the equivalent quantities as projected physical lengths (pc) at the heliocentric distance, $D$.

The initial masses of satellite member stars were drawn from a Chabrier (2001) initial mass function (IMF), which has been found to be a reasonable description of known satellite galaxies (Simon 2019). Initial stellar masses were used to assign current absolute magnitudes from a Bressan et al. (2012) isochrone. When sampling from the IMF, the lower mass bound was set to the hydrogen-burning limit of $0.08 M_{☉}$ and the upper bound was set by the star with the largest initial mass in the evolved isochrone (white dwarfs are ignored). Using the Bressan et al. (2012) isochrones, we transformed from initial stellar mass to current absolute magnitude in the g and r bands for each survey, and then to apparent magnitudes using the distance modulus of the simulated satellite. We applied interstellar extinction to the apparent magnitudes of each simulated star using the same reddening coefficients described in Section 3.

We estimated the photometric uncertainty on the simulated stellar magnitudes based on the depth of the survey at the location of each star according to the formula

$$ \sigma_m = 0.01 + 10^{f(\Delta m)}. \qquad (1) $$

The function, $f(\Delta m)$, maps the difference between the apparent magnitude of a star and the $10\sigma$ survey magnitude limit at the location of the star, $\Delta m = m_{lim} - m$, to the median magnitude uncertainty. We derive $f(\Delta m)$ for each survey by calculating the median magnitude uncertainty as a function of magnitude and magnitude limit. In the middle panel of Figure 3, we plot our photometric uncertainty model as a function of r-band magnitude, given the characteristic depth of DES ($r_{lim} = 23.7$) and PS1 ($r_{lim} = 22.4$).

To assess the sensitivity of our search algorithms, we inserted simulated stellar catalogs for each satellite into the real data and ran our satellite search algorithms at the location of each injected satellite. We simulated stellar catalogs for $10^5$ ($10^6$) satellites in the DES (PS1) footprint.<fnref n="47" /> To make economical use of compute time and simulated data volume, satellites with high surface brightness $\mu < 23.5$ mag arcsec$^{-2}$ and $>10^3$ detected stars brighter than $g = 22$ were not fully simulated and were instead assumed to be detected if they reside within the geometric survey coverage masks (see Appendix B for details). We record the detection significance of each simulated satellite, along with metadata about the survey characteristics at the injected location.

When analyzing the simulated satellites, we use the same configuration that was used to search the real data. However, to save on computational time, we fixed the spatial location and distance modulus of our analysis to the value of the search grid that best matched the location and distance of the simulated satellite. This yields a conservative estimate of the detection significance, since we are ignoring the possibility that background fluctuations could slightly enhance the detection significance at other locations or distances. To assess the impact of this choice, we freed the distance modulus for a small set of simulated satellites and found that the detection probability increased by, at most, a few percentage points for satellites close to the detection threshold.

Our catalog-level insertion procedure does not account for effects of blending in regions of high object density that might affect the detection and/or photometric measurements of member stars. However, the constraints that we placed on the number of bright member stars and surface brightness typically limit our simulated satellite population to surface densities below a few stars per square arcminute (Appendix B). Based on studies of the performance of the DESDM pipeline in crowded regions, blending will not substantially decrease the detectability of satellite galaxies with these surface densities (Wang et al. 2019). In addition, diffuse light from unresolved stars is a subdominant component of the flux for resolved systems at these distances. These assumption are violated for bright nearby globular clusters and classical dwarf galaxies, but we assert that searches are complete for such systems in our survey area.

## 5. Search Algorithms

Milky Way satellites are detected as arcminute-scale overdensities of old, metal-poor stars located in the outer halo of the Milky Way. The brightest satellites were predominantly discovered in visual searches of photographic plates (Shapley 1938a, 1938b; Harrington & Wilson 1950; Wilson 1955; Cannon et al. 1977; Irwin et al. 1990; Ibata et al. 1994). The advent of large digital sky surveys enabled the discovery of fainter systems using statistical matched-filter techniques (Willman et al. 2005a, 2005b; Belokurov et al. 2006, 2007, 2008, 2009, 2010; Grillmair 2006, 2009; Sakamoto & Hasegawa 2006; Zucker et al. 2006a, 2006b; Irwin et al. 2007; Walsh et al. 2007). Matched-filter searches have been broadly applied to the current generation of large surveys to detect larger, fainter, and more distant systems (Kim et al. 2015a, 2015b; Kim & Jerjen 2015b; Koposov et al. 2015, 2018; Laevens et al. 2015a, 2015b;

<footnote n="46">⁴⁶ Assuming a different spatial profile (e.g., an exponential profile) has little impact on detectability; however, see Moskowitz & Walker (2019) for a detailed analysis of alternative spatial kernels.</footnote>

<footnote n="47">⁴⁷ On average, the DES simulations generate many more member stars per satellite due to the deeper DES imaging (Figure 3). This makes it computationally challenging to simulate more DES satellites.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (Left) Stellar completeness models for the DES Y3A2 and PS1 DR1 data sets. The effects of both object detection and star-galaxy classification are included. The dotted and solid curves for PS1 indicate stellar completeness with and without the S/N selection criteria, respectively. (Center) Photometric uncertainty models for both data sets. Dashed lines correspond to statistical uncertainty alone, while the solid lines impose a minimum photometric uncertainty of 0.01 mag. (Right) Color-magnitude distribution of stars in a simulated satellite with heliocentric distance *D* = 28.8 kpc and stellar mass *M*<sub>*</sub> = 880 M<sub>☉</sub> as realized in the DES (red) and PS1 (blue) surveys. Open circles for PS1 indicate stars excluded by the S/N selection criteria. DES field stars in a region of 0.1 deg radius are shown for comparison (gray).</figcaption></figure>

Martin et al. 2015; Homma et al. 2016, 2018, 2019; Torrealba
et al. 2016a, 2016b, 2018, 2019b; Luque et al. 2017; Mau et al.
2020). In addition, maximum-likelihood-based algorithms have
been developed to simultaneously combine morphological and
photometric information to increase sensitivity (Bechtol et al.
2015; Drlica-Wagner et al. 2015, 2016).

Our search employed two automated algorithms to detect low-
surface-brightness, arcminute-scale stellar overdensities. The
first algorithm uses a conventional matched-filter approach,
while the second uses a more complex maximum-likelihood
framework. Both search algorithms were optimized to detect old,
metal-poor stellar populations using their distinct locus in color-
magnitude space. Our search focused on conventional ultra-faint
galaxies and has slightly reduced efficiency for very large stellar
systems (e.g., Torrealba et al. 2016b, 2019b; Pieres et al. 2017),
or especially young and/or metal-rich systems (e.g., Torrealba
et al. 2019a). Our search was also optimized for high Galactic
latitude, where the foreground stellar density does not vary
significantly over degree scales. Importantly, our two search
methods employ different strategies to evaluate the local stellar
density and to filter candidate member stars of Milky Way
satellites according to their spatial and color measurements.

5.1. *Spatial Matched-filter Search*

The first search algorithm, simple, is inspired by the
matched-filter methods of Koposov et al. (2008) and Walsh
et al. (2009), and uses a simple isochrone filter to enhance the
contrast of halo substructures at a given distance relative to the
foreground field of Milky Way stars. The specific implementa-
tion builds upon the technique described by Bechtol et al.
(2015) and Drlica-Wagner et al. (2015).<fnref n="48" /> When analyzing
DES Y3A2, we required that objects be detected in both *g* and *r*
bands and be brighter than *g* = 24.5 mag. When analyzing PS1
DR1, we adopted a signal-to-noise threshold S/N > 10 in the *r*
band. A matched-filter search for spatial overdensities of old,
metal-poor stars was performed, scanning in distance modulus
from 16 ≤ *m* − *M* ≤ 24 mag (16 ≤ *m* − *M* ≤ 22 mag) for
DES Y3A2 (PS1 DR1) in steps of 0.5 mag. These searches

correspond to heliocentric distances of 16 kpc ≪ *D* ≪ 620 kpc and 16 kpc ≪ *D* ≪ 251 kpc, respectively.<fnref n="49" /> At each distance modulus, we selected stars with *g*- and *r*-band magnitudes consistent with the synthetic isochrone of Bressan et al. (2012) with metallicity *Z* = 0.0001 and age *τ* = 12 Gyr. We required that the color difference between each star and the template isochrone be Δ(*g* − *r*) < √0.1² + *σ*²<sub>*g*</sub> + *σ*²<sub>*r*</sub>, where *σ*<sub>*g*</sub> and *σ*<sub>*r*</sub> are the statistical uncertainties on the *g*- and *r*-band magnitudes, respectively.

The survey footprint was partitioned into HEALPix pixels
of nside = 32 (~3.4 deg²) for individual analysis. For each
nside = 32 pixel and distance modulus step, we applied the
isochrone filter described previously and created a map of the
filtered stellar density field, including the central pixel of
interest along with the eight surrounding HEALPix pixels. The
eight surrounding pixels were used to more accurately estimate
the average stellar density in the central pixel of interest. The
filtered stellar density field in the central pixel was smoothed by
a Gaussian kernel (σ = 2'), and we identified local density
peaks by iteratively raising a density threshold until there are
fewer than 10 disconnected regions above the threshold value.
In practice, only the most prominent of these stellar over-
densities passed our minimal statistical significance thresholds.

At the central location of each density peak, we determined
the angular size of a surrounding aperture that maximizes the
significance of the density peak with respect to the distribution
of field stars. Specifically, we iterate through circular apertures
with radii from 1' to 18', and for each radius, we compute the
Poisson significance for the observed stellar counts within the
aperture given the local field density. The local field density is
estimated from an annulus between 18' and 30' surrounding the
peak. When calculating the stellar density, we account for the
coverage of the survey, which is mapped at square arcminute
scales, as described in Sections 3.1 and 3.2. After consolidating
spatially coincident peaks at different distance moduli, all
peaks with Poisson significance SIG > 5.5σ are considered
seeds for subsequent analysis. simple has a high significance

<footnote n="48">⁴⁸ https://github.com/DarkEnergySurvey/simple</footnote>

<footnote n="49">⁴⁹ Our search is less sensitive to systems at larger distances where the apparent magnitude of horizontal branch stars is fainter than the detection limit of the surveys.</footnote>
ceiling at SIG = 37.5σ, corresponding to the numerical limit of the inverse survival function of the normal distribution implemented in scipy.

## 5.2. Likelihood-based Search

The second search algorithm employs a likelihood-based approach implemented with the ugali framework (Bechtol et al. 2015; Drlica-Wagner et al. 2015).<fnref n="50" /> A likelihood function is constructed from the product of Poisson probabilities to detect individual stars based upon their spatial positions, measured fluxes, photometric uncertainties, and local imaging depth, given a model that includes a putative dwarf galaxy and empirical estimation of the local stellar field population (see Appendix C for more details). When calculating the likelihood, we account for both missing survey area and local depth variations mapped on square-arcminute scales (Sections 3.1 and 3.2). We assumed a radially symmetric Plummer profile, scanning over half-light radii, $r_h = \{1.2, 4.2, 9.0\}$, and a spectral model composed of four Bressan et al. (2012) isochrones of $\tau = \{10\,\text{Gyr}, 12\,\text{Gyr}\}$ and $Z = \{0.0001, 0.0002\}$, each weighted by a Chabrier (2001) IMF. This spatial-spectral template was rastered over a spatial grid of HEALPix pixels (nside = 4096; spatial resolution of ~0.7') and range of distance moduli from $16 < m - M < 23$ (heliocentric distances of 16 kpc $< D < 400$ kpc) in steps of 0.5 mag. At each coordinate, we evaluated the likelihood ratio between models with and without a candidate satellite galaxy to generate a three-dimensional map of detection significance. We define a test statistic, $TS = -2\Delta\log(\mathcal{L})$, as our criterion for detection. In the asymptotic limit, the TS will follow a $\chi^2$-distribution with $n$ degrees of freedom (Wilks 1938; Chernoff 1954). In our case, the grid scan maximizes over a grid of satellite sky location, distance, richness, and size, yielding $n \sim 5$, and our threshold of $\sqrt{TS} > 6$ corresponds to a statistical significance of ~4.9σ. Isolated peaks in the TS map were extracted as seeds for further characterization.

# 6. Search Results

The DES and PS1 searches each returned several thousand “seeds” (i.e., locations where the local significance exceeds a minimum threshold). The observed distribution of detection significance falls steeply for each algorithm (Figure 4), with the majority of high-significance seeds coinciding with real resolved stellar systems, regions of spatially varying Galactic extinction or stellar density, or survey artifacts caused by bright stars and incomplete coverage. To robustly infer the population of Milky Way satellites, we define additional criteria to produce a high-purity sample of “detected” satellite candidates. These criteria were self-consistently applied to both survey data and simulations.

## 6.1. Detection Criteria

We developed a set of detection criteria intended to refine the set of raw seeds to a list of high-quality satellite galaxy candidates. These criteria were partially motivated by the observed distribution of seeds and the recovery of previously known satellites. We sought to design a set of selection criteria for which our recovery of known satellites was relatively complete, but restrictive enough that only a small number of additional objects would pass our criteria and could be visually inspected.

Our detection criteria can be broadly categorized into four different types: (i) a set of geometric criteria intended to mask known stellar systems and problematic regions of the survey (Figure 2); (ii) a detection significance threshold; (iii) a spatial match between seeds from the simple and ugali searches; and (iv) a visual inspection and masking of residual survey artifacts. These criteria mimic the criteria applied for satellite discovery in previous studies (i.e., Bechtol et al. 2015; Drlica-Wagner et al. 2015) and are described in more detail below.

(i) *Geometric Criteria*. We applied a set of geometric masks to exclude regions of the survey footprint where systematic features result in the detection of spurious seeds. We began by restricting the seeds from both surveys to regions of low interstellar extinction, $E(B-V) < 0.2$ (Schlegel et al. 1998). While our algorithms incorporate the effects of reddening, regions of high reddening generally trace regions of high Milky Way stellar density. In these regions, the reddening and the stellar field density can vary over relatively small spatial scales, and we find the incidence of false-positive seeds is disproportionately high. This mask removes ~7400 deg² from the PS1 footprint and is negligible for DES (Figure 2).

Our empirical model of the PS1 DR1 footprint is inaccurate near the survey boundaries, resulting in the detection of spurious seeds due to mis-estimation of the stellar density. To remove these spurious seeds, we applied a declination selection of $\delta > -25$ deg. We also removed regions of the PS1 footprint where very high stellar density led to memory overflow issues during application of the likelihood search algorithm. These regions largely overlap with the reddening mask and only remove ~20 deg² of additional area.

We also masked regions around astronomical objects that are known to produce spurious seeds. These masks can generally be separated into regions around bright stars (Hoffleit & Jaschek 1991), Milky Way globular clusters (Harris 1996, 2010 edition), open clusters (WEBDA)<fnref n="51" />, and nearby galaxies that are resolved into individual stars (Nilson 1973; Webbink 1985; Corwin 2004; Bica et al. 2008; Kharchenko et al. 2013). We also mask regions around overdensities in two narrow stellar streams, ATLAS (Koposov et al. 2014; Shipp et al. 2018) and Phoenix (Balbinot et al. 2016). For extended objects, we masked regions consistent with the half-light radii of these objects, with a minimum masked radius of 0.05 deg. For bright stars and objects where size information is unavailable, we masked a circular region with a 0.1 deg radius.

(ii) *Significance Threshold*. We require a fiducial significance threshold of $SIG > 6$ for the simple search algorithm and $\sqrt{TS} > 6$ for the ugali search algorithm. These significance thresholds were chosen such that the observed number of unassociated seeds increases rapidly if the threshold is reduced (Figure 4). In addition, most of the seeds above these thresholds can be readily classified as either genuine stellar systems or obvious survey artifacts (Appendix D), whereas seeds below these thresholds are more often ambiguous.

(iii) *Detection by Both Algorithms*. We required that seeds be detected above the stated significance thresholds by both the simple and ugali search algorithms. To apply this criteria, we matched seeds between the two searches. We defined two seeds to be matched if their centroids were within 12.0 (0.2 deg). We reiterate that our objective is to derive a statistically rigorous observational selection function that yields

<footnote n="50">50 https://github.com/DarkEnergySurvey/ugali</footnote>

<footnote n="51">51 https://webda.physics.muni.cz</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution of detection significance for seeds returned by the *simple* and *ugali* search algorithms in our DES and PS1 satellite searches. Histograms show seeds that survive successive application of selection criteria, including the geometric masks shown in Figure 2. The solid vertical line denotes our fiducial detection significance threshold for each survey. Note that the *simple* algorithm returns a maximum significance of 37.5.</figcaption></figure>

a high-purity list of candidates. Requiring detection by both algorithms greatly increases the purity of our candidate list at a moderate cost to detection efficiency (Figure 4). We find that achieving similar purity with a single algorithm would require a higher detection threshold and would result in lower overall efficiency. Requiring detection by both algorithms does result in the rejection of several known satellites in the PS1 footprint that were only detected by *ugali* (i.e., Aquarius II, Columba I, Leo IV, Leo V, Pisces II, Ursa Major I) or *simple* (Boötes II). We discuss these specific cases in Section 6.2.

(iv) *Residual Survey Artifacts*. We identified several seeds in the PS1 DR1 footprint that passed the previous selection criteria but were clearly associated with imaging or processing artifacts (e.g., stray and scattered light around bright stars, abrupt changes in survey depth, inaccurate survey coverage map, PSF fitting failures). These regions were visually identified and masked with a circular mask of radius of ~0.3–0.5 deg. An additional 12 seeds in the PS1 DR1 passed our fiducial selection criteria and were visually inspected. These seeds showed a poorly defined stellar sequence in color–magnitude space and/or poorly defined spatial morphology. These regions likely result from a combination of less obvious survey artifacts (i.e., mis-estimation of the sky background, excess sensor noise, or low level scattered light) and contamination from background galaxy clusters. We mask regions of ~0.3 deg around each of these seeds. More details on the identification of these regions can be found in Appendix D.

The total search area after masking is 4844 deg² in DES, 21,123 deg² in PS1, and 24,343 deg² together (DES and PS1 overlap in a region of ~1600 deg²). Figure 4 shows the effect

of the successive selections and demonstrates that our two independent search algorithms yield reasonably consistent results, particularly for objects that pass the selection criteria described above. The consistency in the number of high-significance objects returned by *both* search algorithms lends confidence to our final list of satellite systems.

## 6.2. Recovery of Real Satellites

We compare the results of our search to the population of confirmed and candidate dwarf galaxies (Table 2). To assemble our catalog of satellites, we augmented McConnachie (2012) with other recently discovered ultra-faint satellites collected in Simon (2019). Our structural parameters are taken primarily from Muñoz et al. (2018), incorporating improved measurements from deeper data when available (as noted in the table). The kinematic classification of ultra-faint satellite systems is notoriously difficult, due to their faintness and small intrinsic velocity dispersions. We thus assemble our catalog from larger satellite systems ($r_{1/2} > 100$ pc) or smaller satellites with low average surface brightness ($10 \text{ pc} \le r_{1/2} \le 100 \text{ pc}$ and $\mu > 25 \text{ mag arcsec}^{-2}$). Classification for systems with $10 \text{ pc} < r_{1/2} < 20 \text{ pc}$ is particularly challenging since velocity dispersions are rarely resolved and classification arguments are often based on non-uniform chemical and structural measurements. We divide systems in this table into kinematically classified dwarf galaxies (class 4); probable dwarf galaxies based on structural or metallicity measurements (class 3); probable star clusters based on structural, age, or metallicity arguments (class 2); and unconfirmed systems that were discovered in shallower DES data but were not detected in
<table><caption>**Table 2**

Confirmed and Candidate Milky Way Satellites</caption><thead><tr><th>(1)<br>Name</th><th>(2)<br>Survey</th><th>(3)<br>Classification</th><th>(4)<br>R.A.<br>(deg)</th><th>(5)<br>Decl.<br>(deg)</th><th>(6)<br>m - M</th><th>(7)<br>a<sub>h</sub><br>(')</th><th>(8)<br>ε</th><th>(9)<br>Distance<br>(kpc)</th><th>(10)<br>r<sub>1/2</sub><br>(pc)</th><th>(11)<br>M<sub>v</sub><br>(mag)</th><th>(12)<br>References</th></tr></thead><tbody><tr><td>Antlia II</td><td></td><td>4</td><td>143.8868</td><td>-36.7673</td><td>20.6</td><td>76.2</td><td>0.38</td><td>132</td><td>2301</td><td>-9.03</td><td>(1)</td></tr><tr><td>Aquarius II</td><td>PS1</td><td>4</td><td>338.4813</td><td>-9.3274</td><td>20.2</td><td>5.1</td><td>0.39</td><td>108</td><td>125</td><td>-4.4</td><td>(2)</td></tr><tr><td>Boötes I</td><td>PS1</td><td>4</td><td>210.0200</td><td>14.5135</td><td>19.1</td><td>9.97</td><td>0.30</td><td>66</td><td>160</td><td>-6.02</td><td>(3)</td></tr><tr><td>Boötes II</td><td>PS1</td><td>4</td><td>209.5141</td><td>12.8553</td><td>18.1</td><td>3.17</td><td>0.25</td><td>42</td><td>33</td><td>-2.94</td><td>(3)</td></tr><tr><td>Boötes III<sup>a</sup></td><td>PS1</td><td>4</td><td>209.3</td><td>26.8</td><td>18.4</td><td>30.0</td><td>0.5</td><td>47</td><td>289</td><td>-5.75</td><td>(4)</td></tr><tr><td>Boötes IV</td><td>PS1</td><td>3</td><td>233.689</td><td>43.726</td><td>21.6</td><td>7.6</td><td>0.64</td><td>209</td><td>277</td><td>-4.53</td><td>(5)</td></tr><tr><td>Canes Venatici I</td><td>PS1</td><td>4</td><td>202.0091</td><td>33.5521</td><td>21.7</td><td>7.12</td><td>0.44</td><td>218</td><td>338</td><td>-8.80</td><td>(3)</td></tr><tr><td>Canes Venatici II</td><td>PS1</td><td>4</td><td>194.2927</td><td>34.3226</td><td>21.0</td><td>1.52</td><td>0.40</td><td>160</td><td>55</td><td>-5.17</td><td>(3)</td></tr><tr><td>Carina</td><td></td><td>4</td><td>100.4065</td><td>-50.9593</td><td>20.1</td><td>10.1</td><td>0.36</td><td>105</td><td>248</td><td>-9.43</td><td>(3)</td></tr><tr><td>Carina II</td><td></td><td>4</td><td>114.1066</td><td>-57.9991</td><td>17.8</td><td>8.69</td><td>0.34</td><td>36</td><td>77</td><td>-4.5</td><td>(6)</td></tr><tr><td>Carina III</td><td></td><td>4</td><td>114.6298</td><td>-57.8997</td><td>17.2</td><td>3.75</td><td>0.55</td><td>28</td><td>20</td><td>-2.4</td><td>(6)</td></tr><tr><td>Centaurus I</td><td></td><td>3</td><td>189.585</td><td>-40.902</td><td>20.3</td><td>2.9</td><td>0.4</td><td>116</td><td>76</td><td>-5.55</td><td>(7)</td></tr><tr><td>Cetus II</td><td>PS1, DES</td><td>3</td><td>19.47</td><td>-17.42</td><td>17.4</td><td>1.9</td><td>&lt;0.4</td><td>30</td><td>17</td><td>0.0</td><td>(8)</td></tr><tr><td>Cetus III</td><td>PS1, DES</td><td>3</td><td>31.331</td><td>-4.270</td><td>22.0</td><td>1.23</td><td>0.76</td><td>251</td><td>44</td><td>-2.5</td><td>(9)</td></tr><tr><td>Columba I</td><td>PS1, DES</td><td>3</td><td>82.86</td><td>-28.01</td><td>21.3</td><td>2.2</td><td>0.3</td><td>183</td><td>98</td><td>-4.2</td><td>(10)</td></tr><tr><td>Coma Berenices</td><td>PS1</td><td>4</td><td>186.7454</td><td>23.9069</td><td>18.2</td><td>5.64</td><td>0.37</td><td>44</td><td>57</td><td>-4.38</td><td>(3)</td></tr><tr><td>Crater II</td><td>PS1</td><td>4</td><td>177.310</td><td>-18.413</td><td>20.4</td><td>31.2</td><td>&lt;0.1</td><td>117</td><td>1066</td><td>-8.2</td><td>(11)</td></tr><tr><td>DES J0225 + 0304</td><td>DES</td><td>1</td><td>36.4267</td><td>3.0695</td><td>16.9</td><td>2.68</td><td>0.61</td><td>24</td><td>12</td><td>-1.1</td><td>(12)</td></tr><tr><td>Draco</td><td>PS1</td><td>4</td><td>260.0684</td><td>57.9185</td><td>19.4</td><td>9.67</td><td>0.29</td><td>76</td><td>180</td><td>-8.71</td><td>(3)</td></tr><tr><td>Draco II</td><td>PS1</td><td>3</td><td>238.174</td><td>64.579</td><td>16.7</td><td>3.0</td><td>0.23</td><td>22</td><td>17</td><td>-0.8</td><td>(13)</td></tr><tr><td>Eridanus II</td><td>DES</td><td>4</td><td>56.0925</td><td>-43.5329</td><td>22.9</td><td>1.77</td><td>0.35</td><td>380</td><td>158</td><td>-7.21</td><td>(3)</td></tr><tr><td>Fornax</td><td>DES</td><td>4</td><td>39.9583</td><td>-34.4997</td><td>20.8</td><td>19.6</td><td>0.29</td><td>147</td><td>707</td><td>-13.46</td><td>(3), (14)</td></tr><tr><td>Grus I</td><td>DES</td><td>3</td><td>344.1797</td><td>-50.18</td><td>20.4</td><td>0.81</td><td>0.45</td><td>120</td><td>21</td><td>-3.47</td><td>(3)</td></tr><tr><td>Grus II</td><td>DES</td><td>3</td><td>331.02</td><td>-46.44</td><td>18.6</td><td>6.0</td><td>&lt;0.2</td><td>53</td><td>92</td><td>-3.9</td><td>(8)</td></tr><tr><td>Hercules</td><td>PS1</td><td>4</td><td>247.7722</td><td>12.7852</td><td>20.6</td><td>5.63</td><td>0.69</td><td>132</td><td>120</td><td>-5.83</td><td>(3)</td></tr><tr><td>Horologium I</td><td>DES</td><td>4</td><td>43.8813</td><td>-54.116</td><td>19.5</td><td>1.59</td><td>0.27</td><td>79</td><td>31</td><td>-3.55</td><td>(3)</td></tr><tr><td>Horologium II</td><td>DES</td><td>3</td><td>49.1077</td><td>-50.0486</td><td>19.5</td><td>2.09</td><td>0.52</td><td>78</td><td>33</td><td>-2.6</td><td>(15)</td></tr><tr><td>Hydra II</td><td></td><td>3</td><td>185.4251</td><td>-31.9860</td><td>20.9</td><td>1.52</td><td>0.24</td><td>151</td><td>58</td><td>-4.60</td><td>(3)</td></tr><tr><td>Hydrus I</td><td></td><td>4</td><td>37.389</td><td>-79.3089</td><td>17.2</td><td>7.42</td><td>0.21</td><td>28</td><td>53</td><td>-4.71</td><td>(16)</td></tr><tr><td>Indus II</td><td>DES</td><td>1</td><td>309.72</td><td>-46.16</td><td>21.7</td><td>2.9</td><td>&lt;0.4</td><td>214</td><td>180</td><td>-4.3</td><td>(8)</td></tr><tr><td>Kim 2</td><td>DES</td><td>2</td><td>317.2020</td><td>-51.1671</td><td>20.0</td><td>0.48</td><td>0.32</td><td>100</td><td>12</td><td>-3.32</td><td>(3)</td></tr><tr><td>Laevens I</td><td>PS1</td><td>2</td><td>174.0668</td><td>-10.8772</td><td>20.8</td><td>0.51</td><td>0.17</td><td>145</td><td>20</td><td>-4.80</td><td>(3)</td></tr><tr><td>LMC</td><td></td><td>4</td><td>80.8938</td><td>-69.7561</td><td>18.5</td><td>323.0</td><td>0.15</td><td>50</td><td>4735</td><td>-18.12</td><td>(17), (18)</td></tr><tr><td>Leo I</td><td>PS1</td><td>4</td><td>152.1146</td><td>12.3059</td><td>22.0</td><td>3.65</td><td>0.3</td><td>254</td><td>226</td><td>-11.78</td><td>(3)</td></tr><tr><td>Leo II</td><td>PS1</td><td>4</td><td>168.3627</td><td>22.1529</td><td>21.8</td><td>2.52</td><td>0.07</td><td>233</td><td>165</td><td>-9.74</td><td>(3)</td></tr><tr><td>Leo IV</td><td>PS1</td><td>4</td><td>173.2405</td><td>-0.5453</td><td>20.9</td><td>2.54</td><td>0.17</td><td>154</td><td>104</td><td>-4.99</td><td>(3)</td></tr><tr><td>Leo V</td><td>PS1</td><td>4</td><td>172.7857</td><td>2.2194</td><td>21.3</td><td>1.00</td><td>0.43</td><td>178</td><td>39</td><td>-4.40</td><td>(3)</td></tr><tr><td>Pegasus III</td><td>PS1</td><td>4</td><td>336.102</td><td>5.405</td><td>21.7</td><td>0.85</td><td>0.38</td><td>215</td><td>42</td><td>-3.4</td><td>(19)</td></tr><tr><td>Phoenix II</td><td>DES</td><td>4</td><td>354.996</td><td>-54.4115</td><td>19.6</td><td>1.49</td><td>0.67</td><td>83</td><td>21</td><td>-3.30</td><td>(3)</td></tr><tr><td>Pictor I</td><td>DES</td><td>3</td><td>70.949</td><td>-50.2854</td><td>20.3</td><td>0.88</td><td>0.63</td><td>114</td><td>18</td><td>-3.45</td><td>(3)</td></tr><tr><td>Pictor II</td><td></td><td>3</td><td>101.180</td><td>-59.897</td><td>18.3</td><td>3.8</td><td>0.13</td><td>46</td><td>47</td><td>-3.2</td><td>(20)</td></tr><tr><td>Pisces II</td><td>PS1</td><td>4</td><td>344.6345</td><td>5.9526</td><td>21.3</td><td>1.12</td><td>0.34</td><td>182</td><td>48</td><td>-4.22</td><td>(3)</td></tr><tr><td>Reticulum II</td><td>DES</td><td>4</td><td>53.9203</td><td>-54.0513</td><td>17.4</td><td>5.52</td><td>0.58</td><td>30</td><td>31</td><td>-3.88</td><td>(3)</td></tr><tr><td>Reticulum III</td><td>DES</td><td>3</td><td>56.36</td><td>-60.45</td><td>19.8</td><td>2.4</td><td>&lt;0.4</td><td>92</td><td>64</td><td>-3.3</td><td>(8)</td></tr><tr><td>Sagittarius</td><td></td><td>4</td><td>283.8313</td><td>-30.5453</td><td>17.1</td><td>342.0</td><td>0.64</td><td>26</td><td>1565</td><td>-13.5</td><td>(21)</td></tr><tr><td>Sagittarius II</td><td>PS1</td><td>4</td><td>298.1647</td><td>-22.0651</td><td>19.2</td><td>1.6</td><td>&lt;0.1</td><td>69</td><td>32</td><td>-5.2</td><td>(22)</td></tr><tr><td>Sculptor</td><td>DES</td><td>4</td><td>15.0183</td><td>-33.7186</td><td>19.6</td><td>11.17</td><td>0.33</td><td>84</td><td>223</td><td>-10.82</td><td>(3), (23)</td></tr><tr><td>Segue</td><td>PS1</td><td>4</td><td>151.7504</td><td>16.0756</td><td>16.8</td><td>3.62</td><td>0.33</td><td>23</td><td>20</td><td>-1.30</td><td>(3)</td></tr><tr><td>Segue 2</td><td>PS1</td><td>3</td><td>34.8226</td><td>20.1624</td><td>17.7</td><td>3.76</td><td>0.22</td><td>35</td><td>34</td><td>-1.86</td><td>(3)</td></tr><tr><td>Sextans</td><td>PS1</td><td>4</td><td>153.2628</td><td>-1.6133</td><td>19.7</td><td>16.5</td><td>0.30</td><td>86</td><td>345</td><td>-8.72</td><td>(3)</td></tr><tr><td>SMC</td><td></td><td>4</td><td>13.1867</td><td>-72.8286</td><td>19.0</td><td>151.0</td><td>0.40</td><td>62</td><td>2728</td><td>-17.18</td><td>(17), (24)</td></tr><tr><td>Triangulum II</td><td>PS1</td><td>3</td><td>33.3252</td><td>36.1702</td><td>17.4</td><td>1.99</td><td>0.46</td><td>30</td><td>13</td><td>-1.60</td><td>(3)</td></tr><tr><td>Tucana II</td><td>DES</td><td>4</td><td>342.9796</td><td>-58.5689</td><td>18.8</td><td>12.59</td><td>0.39</td><td>58</td><td>165</td><td>-3.8</td><td>(25)</td></tr><tr><td>Tucana III</td><td>DES</td><td>3</td><td>359.15</td><td>-59.60</td><td>17.0</td><td>6.0</td><td>0.0</td><td>25</td><td>44</td><td>-2.4</td><td>(8)</td></tr><tr><td>Tucana IV</td><td>DES</td><td>4</td><td>0.73</td><td>-60.85</td><td>18.4</td><td>11.8</td><td>0.4</td><td>48</td><td>128</td><td>-3.5</td><td>(8)</td></tr><tr><td>Tucana V</td><td>DES</td><td>3</td><td>354.35</td><td>-63.27</td><td>18.7</td><td>1.8</td><td>0.7</td><td>55</td><td>16</td><td>-1.6</td><td>(8)</td></tr><tr><td>Ursa Major I</td><td>PS1</td><td>4</td><td>158.7706</td><td>51.9479</td><td>19.9</td><td>8.31</td><td>0.59</td><td>97</td><td>151</td><td>-5.12</td><td>(3)</td></tr><tr><td>Ursa Major II</td><td>PS1</td><td>4</td><td>132.8726</td><td>63.1335</td><td>17.5</td><td>13.8</td><td>0.56</td><td>32</td><td>85</td><td>-4.25</td><td>(3)</td></tr><tr><td>Ursa Minor</td><td>PS1</td><td>4</td><td>227.2420</td><td>67.2221</td><td>19.4</td><td>18.3</td><td>0.55</td><td>76</td><td>272</td><td>-9.03</td><td>(3)</td></tr><tr><td>Virgo I</td><td>PS1</td><td>3</td><td>180.038</td><td>-0.681</td><td>19.8</td><td>1.76</td><td>0.59</td><td>91</td><td>30</td><td>-0.33</td><td>(9)</td></tr><tr><td>Willman I</td><td>PS1</td><td>4</td><td>162.3436</td><td>51.0501</td><td>17.9</td><td>2.51</td><td>0.47</td><td>38</td><td>20</td><td>-2.53</td><td>(3)</td></tr></tbody></table>

Notes. Column descriptions: (1) satellite name; (2) survey(s) in which the system resides; (3) system classification (see below); (4), (5) published R.A. and decl.; (6)–(8) published distance modulus, observed semimajor axis of an ellipse containing half of the light, and ellipticity; (9), (10) derived heliocentric distance and azimuthally averaged physical half-light radius; (11) published absolute V-band magnitude; (12) literature reference. When two references are listed, the second was used for the distance measurement. Classifications are: (1) unconfirmed systems, (2) probable star clusters, (3) probable dwarfs, (4) kinematically confirmed dwarfs. Literature references are: (1) Torrealba et al. (2019b), (2) Torrealba et al. (2016a), (3) Muñoz et al. (2018), (4) Grillmair (2009), (5) Homma et al. (2019), (6) Torrealba et al. (2018), (7) Mau et al. (2020), (8) Drlica-Wagner et al. (2015), (9) Homma et al. (2018), (10) Carlin et al. (2017), (11) Torrealba et al. (2016b), (12) Luque et al. (2017), (13) Longeard et al. (2018), (14) Pietrzyński et al. (2009), (15) Kim & Jerjen (2015b), (16) Koposov et al. (2018), (17) Makarov et al. (2014), (18) Pietrzyński et al. (2013), (19) Kim et al. (2016), (20) Drlica-Wagner et al. (2016), (21) McConnachie (2012), (22) Mutlu-Pakdil et al. (2018), (23) Martínez-Vázquez et al. (2015), (24) Graczyk et al. (2014), (25) Koposov et al. (2015).

a Approximate half-light radius and ellipticity estimated from Grillmair (2009).

(This table is available in machine-readable form.)
<table><caption>**Table 3**
Recovery of Confirmed and Candidate Milky Way Satellite Galaxies in DES Y3A2</caption><thead><tr><th>(1)<br>Name</th><th>(2)<br>√TS</th><th>(3)<br>SIG</th><th>(4)<br>P<sub>det</sub></th><th>(5)<br>Distance<br>(kpc)</th><th>(6)<br>r<sub>1/2</sub><br>(pc)</th><th>(7)<br>M<sub>V</sub><br>(mag)</th><th>(8)<br>ρ<sub>*</sub><br>(arcmin<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Fornax</td><td>480.07</td><td>37.5</td><td>1.00</td><td>147</td><td>707</td><td>-13.46</td><td>16.52</td></tr><tr><td>Sculptor</td><td>415.08</td><td>37.5</td><td>1.00</td><td>84</td><td>223</td><td>-10.82</td><td>5.87</td></tr><tr><td>Reticulum II</td><td>54.56</td><td>37.5</td><td>1.00</td><td>30</td><td>31</td><td>-3.88</td><td>1.17</td></tr><tr><td>*Eridanus II</td><td>36.17</td><td>27.41</td><td>1.00</td><td>380</td><td>158</td><td>-7.21</td><td>1.10</td></tr><tr><td>Tucana II</td><td>23.85</td><td>12.88</td><td>0.91</td><td>58</td><td>165</td><td>-3.8</td><td>1.87</td></tr><tr><td>Grus II</td><td>21.82</td><td>13.29</td><td>0.97</td><td>53</td><td>92</td><td>-3.9</td><td>2.03</td></tr><tr><td>Horologium I</td><td>21.78</td><td>20.68</td><td>0.99</td><td>79</td><td>31</td><td>-3.55</td><td>1.24</td></tr><tr><td>Tucana III</td><td>18.60</td><td>12.72</td><td>0.91</td><td>25</td><td>44</td><td>-2.4</td><td>1.62</td></tr><tr><td>Tucana IV</td><td>16.03</td><td>10.88</td><td>0.91</td><td>48</td><td>128</td><td>-3.5</td><td>1.60</td></tr><tr><td>Phoenix II</td><td>13.36</td><td>12.18</td><td>0.98</td><td>83</td><td>21</td><td>-3.30</td><td>1.35</td></tr><tr><td>Horologium II</td><td>11.62</td><td>11.21</td><td>0.87</td><td>78</td><td>33</td><td>-2.6</td><td>1.07</td></tr><tr><td>Tucana V</td><td>11.27</td><td>9.97</td><td>0.89</td><td>55</td><td>16</td><td>-1.6</td><td>1.78</td></tr><tr><td>Pictor I</td><td>10.91</td><td>8.55</td><td>0.97</td><td>114</td><td>18</td><td>-3.45</td><td>1.50</td></tr><tr><td>Columba I</td><td>10.67</td><td>9.45</td><td>0.33</td><td>183</td><td>98</td><td>-4.2</td><td>2.09</td></tr><tr><td>Cetus II</td><td>10.47</td><td>7.18</td><td>0.62</td><td>30</td><td>17</td><td>0.0</td><td>0.79</td></tr><tr><td>Grus I</td><td>9.78</td><td>8.80</td><td>0.97</td><td>120</td><td>21</td><td>-3.47</td><td>1.57</td></tr><tr><td>*Kim 2</td><td>8.17</td><td>9.31</td><td>0.93</td><td>100</td><td>12</td><td>-3.32</td><td>3.13</td></tr><tr><td>Reticulum III</td><td>8.12</td><td>7.46</td><td>0.92</td><td>92</td><td>64</td><td>-3.3</td><td>1.46</td></tr><tr><td>*Cetus III</td><td>3.96</td><td>3.68</td><td>0.01</td><td>251</td><td>44</td><td>-2.5</td><td>0.89</td></tr><tr><td>*Indus II</td><td>3.86</td><td>3.65</td><td>0.01</td><td>214</td><td>180</td><td>-4.3</td><td>4.02</td></tr><tr><td>*DES J0225+0304</td><td>...</td><td>...</td><td>0.94</td><td>24</td><td>12</td><td>-1.1</td><td>0.98</td></tr></tbody></table>

Note. Column descriptions: (1) satellite name; (2) square-root of the test statistic from ugali search; (3) statistical significance value from simple search (maximum of 37.5); (4) detection probability from survey selection function; (5), (6) heliocentric distance and azimuthally averaged physical half-light radius, calculated from observed parameters listed in Table 2; (7) published absolute V-band magnitude; (8) local stellar density. Satellites denoted with asterisks are not included in the statistical sample used to derive the luminosity function.

(This table is available in machine-readable form.)

our search (class 1). Class 2 notably includes Crater I/Laevens 1 (Belokurov et al. 2014; Laevens et al. 2014) and Kim 2 (Kim et al. 2015b), which have been proposed to be star clusters due to structural, age, and metallicity arguments (Kirby et al. 2015; Weisz et al. 2016), but pass our selection on size and surface brightness. Our primary sample of probable and kinematically classified satellite galaxies are systems with class ≥3. Our search recovers the majority of previously discovered Milky Way satellite galaxies in the PS1 DR1 and DES Y3A2 footprints (Tables 3 and 4).

We recovered 18 out of 21 confirmed and candidate satellite galaxies in the DES footprint above our nominal threshold of $\sqrt{TS} > 6$ and $SIG > 6$. Two of the ultra-faint galaxy candidates initially detected in DES Y2Q1 data were marked as “lower-confidence” candidates due to their locations in regions of non-uniform survey coverage during the first two years of DES observations (Drlica-Wagner et al. 2015). With the deeper and more homogeneous imaging of the Y3A2 data set, Cetus II (DES J0117–1725) is detected with $\sqrt{TS} = 10.5$ and $SIG = 7.2$. This statistical significance is comparable to other confirmed candidates (e.g., Columba I and Pictor I). Meanwhile, the second low-confidence candidate, Indus II (DES J2038–4609), drops below our detection threshold.<fnref n="52" /> We recovered the Eridanus III system reported in Bechtol et al. (2015) and Koposov et al. (2015) with high significance ($\sqrt{TS} = 9.5$ and $SIG = 9.0$); however, due to its small physical size

$(r_{1/2} = 5 \text{ pc}, \text{Conn et al. 2018})$, we do not include it in our list of confirmed and candidate dwarf galaxies. Neither of the two objects reported in Luque et al. (2017) was detected as a significant seed in our automated search, and visual inspection of the DES Y3A2 data coincident with these candidates did not reveal any significant excesses. Cetus III, an ultra-faint galaxy candidate identified in early data from HSC SSP (Homma et al. 2018), is located within the DES footprint but falls below our detection threshold. This is expected, given the large distance (251 kpc) and low azimuthally averaged surface brightness (~29.3 mag arcsec$^{-2}$) of this candidate.

In comparison, we recovered 20 of the 32 confirmed and candidate satellite galaxies known to reside in the PS1 DR1 footprint. The lower recovery rate in PS1 is expected since many of the satellites in the PS1 footprint were discovered with significantly deeper data. Of the 12 satellites that fall below our detection threshold, 5 were discovered in deeper surveys and are not expected to be detected by PS1: Boötes IV, Cetus III, and Virgo I were discovered in HSC SSP; Columba I was discovered in DES; and Aquarius II was discovered in VST ATLAS. The seven remaining satellites were discovered using data from SDSS, but several of these objects were confirmed with deeper follow-up observations before publication. Leo V had deep follow-up imaging from the Isaac Newton Telescope and spectroscopy from the Hectochelle fiber spectrograph at the Multiple Mirror Telescope (Belokurov et al. 2008), while Pisces II had follow-up imaging from the MOSAIC camera at the 4 m Mayall Telescope (Belokurov et al. 2010). Pegasus III was announced after deep follow-up observations with DECam (Kim & Jerjen 2015a). Leo IV was discovered in data from SDSS without additional follow-

<footnote n="52">⁵² Deep follow-up imaging with Magellan/Megacam supports the hypothesis that Indus II (DES J2038–4609) is a chance alignment of stars (S. A. Cantu et al. 2020, in preparation).</footnote>
<table><caption>**Table 4**
Recovery of Confirmed and Candidate Milky Way Satellite Galaxies in PS1 DR1</caption><thead><tr><th>(1)<br>Name</th><th>(2)<br>√TS</th><th>(3)<br>SIG</th><th>(4)<br>P<sub>det</sub></th><th>(5)<br>Distance<br>(kpc)</th><th>(6)<br>r<sub>1/2</sub><br>(pc)</th><th>(7)<br>M<sub>V</sub><br>(mag)</th><th>(8)<br>ρ<sub>*</sub><br>(arcmin<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Leo I</td><td>157.63</td><td>37.5</td><td>1.00</td><td>254</td><td>226</td><td>-11.78</td><td>1.18</td></tr><tr><td>Leo II</td><td>104.05</td><td>37.5</td><td>1.00</td><td>233</td><td>165</td><td>-9.74</td><td>3.09</td></tr><tr><td>Draco</td><td>96.94</td><td>37.5</td><td>1.00</td><td>76</td><td>180</td><td>-8.71</td><td>3.02</td></tr><tr><td>Ursa Minor</td><td>83.14</td><td>37.5</td><td>1.00</td><td>76</td><td>272</td><td>-9.03</td><td>3.19</td></tr><tr><td>Sextans</td><td>58.62</td><td>24.63</td><td>1.00</td><td>86</td><td>345</td><td>-8.72</td><td>3.00</td></tr><tr><td>Canes Venatici I</td><td>36.00</td><td>25.33</td><td>1.00</td><td>218</td><td>338</td><td>-8.80</td><td>1.01</td></tr><tr><td>Boötes I</td><td>25.29</td><td>11.63</td><td>0.95</td><td>66</td><td>160</td><td>-6.02</td><td>1.29</td></tr><tr><td>Ursa Major II</td><td>18.66</td><td>8.86</td><td>0.94</td><td>32</td><td>85</td><td>-4.25</td><td>1.76</td></tr><tr><td>Coma Berenices</td><td>15.29</td><td>9.75</td><td>0.93</td><td>44</td><td>57</td><td>-4.38</td><td>1.07</td></tr><tr><td>Sagittarius II</td><td>15.19</td><td>11.66</td><td>0.55</td><td>69</td><td>32</td><td>-5.2</td><td>12.60</td></tr><tr><td>Willman 1</td><td>15.03</td><td>12.54</td><td>0.54</td><td>38</td><td>20</td><td>-2.53</td><td>0.95</td></tr><tr><td>Canes Venatici II</td><td>11.70</td><td>8.78</td><td>0.93</td><td>160</td><td>55</td><td>-5.17</td><td>0.98</td></tr><tr><td>Segue 1</td><td>10.79</td><td>8.55</td><td>0.48</td><td>23</td><td>20</td><td>-1.30</td><td>1.22</td></tr><tr><td>Segue 2</td><td>10.75</td><td>7.25</td><td>0.05</td><td>35</td><td>34</td><td>-1.86</td><td>1.55</td></tr><tr><td>Crater II</td><td>10.42</td><td>6.08</td><td>0.06</td><td>117</td><td>1066</td><td>-8.2</td><td>2.44</td></tr><tr><td>*Ursa Major I</td><td>10.18</td><td>5.99</td><td>0.24</td><td>97</td><td>151</td><td>-5.12</td><td>0.99</td></tr><tr><td>*Laevens 1</td><td>9.89</td><td>9.52</td><td>0.96</td><td>145</td><td>20</td><td>-4.8</td><td>1.67</td></tr><tr><td>Draco II</td><td>9.76</td><td>7.90</td><td>0.24</td><td>22</td><td>17</td><td>-0.8</td><td>1.63</td></tr><tr><td>Triangulum II</td><td>9.46</td><td>6.76</td><td>0.23</td><td>30</td><td>13</td><td>-1.60</td><td>3.04</td></tr><tr><td>Hercules</td><td>9.11</td><td>6.44</td><td>0.44</td><td>132</td><td>120</td><td>-5.83</td><td>3.91</td></tr><tr><td>*Leo IV</td><td>8.25</td><td>4.94</td><td>0.18</td><td>154</td><td>104</td><td>-4.99</td><td>1.48</td></tr><tr><td>Cetus II</td><td>7.42</td><td>6.14</td><td>0.02</td><td>30</td><td>17</td><td>0.0</td><td>1.08</td></tr><tr><td>*Aquarius II</td><td>7.27</td><td>5.07</td><td>0.01</td><td>108</td><td>125</td><td>-4.4</td><td>1.67</td></tr><tr><td>*Leo V</td><td>6.95</td><td>4.14</td><td>0.14</td><td>178</td><td>39</td><td>-4.4</td><td>1.31</td></tr><tr><td>*Pisces II</td><td>6.25</td><td>4.39</td><td>0.03</td><td>182</td><td>48</td><td>-4.22</td><td>1.58</td></tr><tr><td>*Columba I<sup>a</sup></td><td>6.15</td><td>5.34</td><td>0.00</td><td>183</td><td>98</td><td>-4.2</td><td>2.46</td></tr><tr><td>*Boötes II</td><td>5.95</td><td>6.46</td><td>0.42</td><td>42</td><td>33</td><td>-2.94</td><td>1.27</td></tr><tr><td>*Boötes IV</td><td>5.41</td><td>4.69</td><td>0.00</td><td>209</td><td>277</td><td>-4.53</td><td>1.56</td></tr><tr><td>*Pegasus III</td><td>4.80</td><td>3.68</td><td>0.00</td><td>215</td><td>42</td><td>-3.4</td><td>2.09</td></tr><tr><td>*Virgo I</td><td>4.06</td><td>4.08</td><td>0.00</td><td>91</td><td>30</td><td>-0.33</td><td>1.58</td></tr><tr><td>*Boötes III<sup>b</sup></td><td>4.00</td><td>4.34</td><td>0.24</td><td>47</td><td>289</td><td>-5.75</td><td>1.16</td></tr><tr><td>*Cetus III</td><td>...</td><td>4.55</td><td>0.00</td><td>251</td><td>44</td><td>-2.5</td><td>1.00</td></tr></tbody></table>

Notes. Column descriptions are the same as Table 3. Satellites denoted with asterisks are not included in the statistical sample used to derive the luminosity function.

(This table is available in machine-readable form.)

up (Belokurov et al. 2007) and is detected significantly above threshold by ugali ($\sqrt{TS} = 8.2$). However, the simple significance ($SIG = 4.9$) falls below our threshold. Similarly, Ursa Major I is detected significantly with ugali, but falls just below the threshold for simple. Boötes II falls just below our threshold for detection and has a comparably low detection probability reported by Koposov et al. (2008). Boötes III is a diffuse object (extending ~1.5 deg) with a complex morphology that was identified visually in filtered stellar density maps from SDSS DR5 (Grillmair 2009). The large size and complex morphology of Boötes III have made it challenging to detect with automated search algorithms (Koposov et al. 2008; Walsh et al. 2009).<fnref n="53" />

We found one candidate in the PS1 DR1 search that passed all selection criteria and is unassociated with our catalogs of known satellites (see Appendix D for more details). This candidate, located at $(\alpha_{2000}, \delta_{2000}) = (247.725, -0.971)$, is a compact, $r_{1/2} = 3.7$ pc, low-luminosity, $M_V = 0.6$, system residing at a heliocentric distance of $D = 15.6$ kpc. While

the physical nature of such a faint system is ambiguous without internal kinematics measurements, the small physical size of this candidate is consistent with other low-luminosity outer-halo star clusters that have been discovered in recent surveys (e.g., Mau et al. 2019; Torrealba et al. 2019a). Due to the small physical size of this candidate, we classify it as a likely star cluster and do not include it in the sample of confirmed and candidates satellite galaxies used for deriving the Milky Way satellite galaxy luminosity function.

The non-detection of several known satellites is not unexpected, given that we prioritize purity over completeness in our selection criteria. Discovery-driven searches often set a lower significance threshold, relying on visual inspection and targeted follow-up observations to reject false positives. This subjective selection function is difficult to characterize with an automated analysis relying on simulations. By setting more restrictive detection criteria, we can be confident that every satellite that passes our automated selection criteria would pass subjective selection. Such a requirement is critical to self-consistently interpret the recovery of simulated satellites and the derived observational selection function.

<footnote><sup>a</sup> Located in a masked region of the PS1 footprint ($\delta < -25$ deg).</footnote>

<footnote><sup>b</sup> Approximate half-light radius and ellipticity estimated from Grillmair (2009).</footnote>

<footnote n="53">⁵³ We also examined the regions around five SDSS candidates proposed by Liu et al. (2008), but we do not find any significant excesses at these locations.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Detection efficiency of searches for Milky Way satellites in DES Y3A2. Detection efficiency ranges from 0% (blue) to 100% (white) and is shown as a function of azimuthally averaged physical half-light radius and absolute V-band magnitude in different bins of heliocentric distance (logarithmically spaced from 8 kpc to 512 kpc). The physical parameters of known satellites located within the DES, and/or PS1 footprints are indicated in black. The black dashed line shows an analytic approximation to the 50% detectability contour, while the red dashed line shows the 50% detection efficiency contour for SDSS DR5 from Koposov et al. (2008). The DES search is significantly more sensitive than the SDSS DR5 search of Koposov et al. (2008).</figcaption></figure>

### 6.3. Recovery of Simulated Satellites

We illustrate the recovery of simulated satellites in the physically motivated parameter space of satellite absolute magnitude, $M_V$; heliocentric distance, $D$; and azimuthally averaged physical half-light radius, $r_{1/2}$. Figures 5 and 6 show the detectability of simulated satellites as a function of $M_V$ and $r_{1/2}$ in six slices of distance. The coloring of each bin corresponds to the fraction of satellites in that bin that pass the detection criteria presented in Section 6.1 (i.e., $\sqrt{TS} > 6$ and $SIG > 6$). The DES simulations generally contain ~30 simulated satellites per bin, while the PS1 simulations contain ~250 simulated satellites per bin. We overplot the physical properties of the known Milky Way satellites residing within the DES and PS1 footprints from Tables 3 and 4. In general, the physical parameters of the known satellites recovered by our search (Section 6.2) are consistent with the sensitivity envelope derived from simulated satellites (Figure 6).

We indicate the approximate surface brightness and absolute magnitude limits of SDSS as derived by Koposov et al. (2008) with a dashed red line. We find that our search on DES Y3A2 is significantly more sensitive than SDSS, while the sensitivity of our PS1 DR1 search is roughly comparable to the SDSS search in many regions of parameter space. When directly comparing PS1 and SDSS catalogs in overlapping fields, and applying quality and star–galaxy separation criteria to both surveys, the stellar efficiency curves of SDSS and PS1 are similar.<fnref n="54" /> The

sensitivity of our PS1 search is largely driven by the conservative detection thresholds we set and the requirement that candidates be detected by both search algorithms. If these restrictions are relaxed, we find that our search becomes significantly more sensitive, but with a corresponding increase in the number of false-positive seeds that need to be rejected using visual inspection. It may be possible to increase the sensitivity of the PS1 search if the incidence of spurious seeds can be reduced by subsequent PS1 data releases (i.e., PS1 DR2) or through a more precise estimate of the PS1 survey coverage. We self-consistently applied the same detection criteria to both the simulated satellites and the real systems when deriving our luminosity function.

## 7. Observational Selection Function

The observational selection function defines the detectability of a satellite as a function of its physical properties and location on the sky. Following the convention of Koposov et al. (2008) and Walsh et al. (2009), we present our observational selection function in terms of the physical parameters of heliocentric distance ($D$), absolute magnitude ($M_V$), and azimuthally averaged projected physical half-light radius ($r_{1/2}$). The detectability of a satellite with a given set of parameters can be predicted directly from our simulations through a nearest-neighbors search (e.g., Knuth 1973). However, nearest-neighbor searches can be unwieldy, and we offer two alternative parameterizations of satellite detectability based on an analytic approximation and a machine-learning classifier. When training the machine-learning classifier, we include the local stellar density as an additional feature for predicting satellite detectability.

<footnote n="54">⁵⁴ Care must be taken in comparing the quoted depths of PS1 and SDSS. SDSS conventionally quotes a 95% completeness limit for point sources of $g$, $r = 22.2$ before star–galaxy separation. This is not directly comparable to the 10$\sigma$ magnitude limit of $g$, $r \sim 22.5$ calculated for PS1 DR1 in Section 3. The completeness depth of each data set is dependent on the exact selection criteria applied to the catalogs.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Similar to Figure 5, but for the PS1 DR1 search. The 50% detectability of the PS1 DR1 search (dashed black line) is significantly less sensitive than DES and is somewhat less sensitive than the 50% detection efficiency of the SDSS DR5 search reported by Koposov et al. (2008) (red dashed line).</figcaption></figure>

<table><caption>**Table 5**
Parameterization of 50% Detectability Contour</caption><thead><tr><th rowspan="2">Distance (kpc)</th><th colspan="3">DES</th><th colspan="3">PS1</th></tr><tr><th>A<sub>0</sub></th><th>M<sub>V,0</sub> (mag)</th><th>log<sub>10</sub>(r<sub>1/2,0</sub>) (log<sub>10</sub>(pc))</th><th>A<sub>0</sub></th><th>M<sub>V,0</sub> (mag)</th><th>log<sub>10</sub>(r<sub>1/2,0</sub>) (log<sub>10</sub>(pc))</th></tr></thead><tbody><tr><td>11.3</td><td>21.5</td><td>7.8</td><td>3.8</td><td>22.8</td><td>7.1</td><td>4.0</td></tr><tr><td>22.6</td><td>24.1</td><td>8.3</td><td>4.2</td><td>19.0</td><td>5.0</td><td>4.1</td></tr><tr><td>45.2</td><td>17.2</td><td>5.2</td><td>4.3</td><td>14.1</td><td>1.8</td><td>4.2</td></tr><tr><td>90.5</td><td>8.6</td><td>1.2</td><td>4.1</td><td>11.0</td><td>-0.3</td><td>4.3</td></tr><tr><td>181</td><td>6.6</td><td>-1.1</td><td>4.1</td><td>7.5</td><td>-2.2</td><td>4.2</td></tr><tr><td>362</td><td>6.3</td><td>-2.3</td><td>4.3</td><td>6.8</td><td>-4.0</td><td>4.4</td></tr></tbody></table>

## 7.1. Analytic Approximation

We first present a simple analytic approximation for the contour defining the parameters of satellites with 50% detection probability, $P_{\text{det}}(D, M_V, r_{1/2}) = 0.5$. We find that at fixed distance, this $P_{\text{det}}$ contour can be well described by

$$ \log_{10}(r_{1/2}) = \frac{A_0(D)}{(M_V - M_{V,0}(D))} + \log_{10}(r_{1/2,0}(D)), \quad (2) $$

where $r_{1/2}$ is in units of pc, $M_V$ is in units of mag, and $D$ is in units of kpc. This equation contains three distance-dependent constants ($A_0, M_{V,0}, r_{1/2,0}$), which were fit to the $P_{\text{det},50}$ contour in each of our six slices of distance (Table 5). In particular, $M_{V,0}$ and $\log_{10}(r_{1/2,0})$ represent asymptotic limits in absolute magnitude and physical half-light radius as a function of satellite distance. The scale parameter, $A_0(D)$, describes the “radius of curvature” of the $P_{\text{det},50}$ contour at a given distance. These fits to $P_{\text{det},50}$ are overplotted as dashed black lines on Figures 5 and 6. The parameters describing $P_{\text{det},50}$ vary smoothly as a function of distance, and interpolating between them can provide a reasonable approximation for $P_{\text{det},50}$ for any distance within the range studied.

Koposov et al. (2008) provided a similar description of satellite detectability in SDSS DR5, in terms of limiting surface brightness, $\mu_{\text{lim}}$, and limiting absolute magnitude, $M_{V,\text{lim}}$ (red dashed lines in Figures 5 and 6). Koposov et al. show that this simple model is sufficient to describe the sensitivity of their search, as derived from a small sample of 8000 simulated galaxies. However, the functional form of this model is insufficient to fully capture the shape of the detectability contours of our search, which was derived from a much larger set of simulated galaxies. Fitting the model of Koposov et al. to our simulations systematically underestimates our sensitivity to faint, compact satellites. The difference between our model and that of Koposov et al. is not unexpected, since the exact shape of the detectability contour depends on the data set and search algorithm.

## 7.2. Machine-learning Model

Walsh et al. (2009) emphasized that $P_{\text{det},50}$ serves as a useful approximation for detectability, but that it does not fully capture the intermediate region between 100% detectability and 0% detectability. Most of the known satellites lie in the region of parameter space of intermediate detection probability, and accordingly, accurate treatment of the detection efficiency gradient in this region is an important component of the interpretation. It is expected that many of the Milky Way satellites lie just beyond the current detection threshold (e.g., Hargis et al. 2014; Jethwa et al. 2018; Newton et al. 2018), and thus any sensitivity beyond the $P_{\text{det},50}$ envelope provides valuable information on this population of faint, distant, and low-surface-brightness satellites.

To more fully encapsulate the results of our simulations, we trained a gradient-boosted decision tree classifier to predict the detectability of a satellite based on its physical properties. This represents a binary classification problem, where we seek to predict the relationship between a set of input features, $X$, and a set of labels, $Y$. We treated each simulated satellite as a training
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fraction of simulated satellites in the test set that pass our detection criteria for DES (left) and PS1 (right) versus the detection probability predicted by our ML classifier. The dashed black one-to-one line indicates perfect performance, and the colorbar indicates the number of simulated satellites in each bin of predicted detection probability. The number of simulated satellites in the training and test sets for each survey is indicated above each plot.</figcaption></figure>

instance, *i*, with a feature vector, $X_i$, composed of the physical properties of the satellites. We labeled satellites as “detected” ($Y_i = 1$) if they satisfied the detection criteria described in Section 6.1 and “undetected” ($Y_i = 0$) otherwise. The output of the machine-learning classifier is the probability that a satellite will be detected.<fnref n="55" /> For each survey, we classified satellites as detected/undetected, depending on whether they pass the detection criteria defined in Section 6.1.

We trained a gradient-boosted decision tree classifier using XGBoost (Chen & Guestrin 2016) and scikit-learn (Pedregosa et al. 2011) as follows:

1. Randomly split the simulated satellites into training and test sets that contain 90% and 10% of the simulated satellites, respectively.

2. Randomly split the training set from the previous step into *k* hold-out cross-validation subsets. We chose *k* = 3 for this analysis by performing a manual grid search over different numbers of cross-validation folds.

3. Train a XGBClassifier using GridSearchCV to select hyperparameters that yield the best test-set classification score. Hyperparameters include the learning rate, number of trees, and maximum tree depth (see Appendix E).

Our feature vector consisted of the absolute magnitude, the heliocentric distance, the azimuthally averaged projected half-light radius, and the stellar density at the location of each simulated satellite in the training set.

We evaluate the performance of our trained classifier using several metrics. To assess the robustness and accuracy of the model, we evaluate the fraction of detected and undetected objects in the test set that are classified correctly. We find that the DES classifier is 97% accurate for both classes, while the PS1 classifier is 94% accurate for detected objects and 97% accurate for undetected objects. The fact that our test-set classification is accurate indicates that the training and test sets (unsurprisingly) represent the same underlying distribution. Moreover, because both detected and undetected objects are classified accurately, we conclude that the algorithm is not systematically biased toward either class. Figure 7 illustrates

the true fraction of detected objects in the test set versus the detection probability predicted by the classifier. Even though the majority of objects are either always detected or never detected for both DES and PS1, our algorithm accurately predicts the detection probability of satellites in the intermediate regime. The region of intermediate detection probability can be attributed to stochasticity in the distribution of stellar fluxes and spatial positions in statistical realizations of a given satellite, as well as local variations in the field population and survey characteristics.

We also trained random forest (RF) classifiers on the same sample of simulated satellites, since RFs provide easily interpretable estimates of feature importance. We trained one RF using a minimal feature vector (absolute magnitude, heliocentric distance, and physical size) and another using a larger set of simulated satellite properties (absolute magnitude, heliocentric distance, physical size, surface brightness, $N(g < 22)$, ellipticity, sky position, and local stellar density). We found that the RF model trained on the full feature vector was slightly more accurate compared to the three-feature RF, although both were biased for high- and low-detection probability objects with respect to our nominal algorithm. The relatively small improvement from using the full feature vector gives us confidence that our nominal feature vector adequately captures much of the necessary information to predict detectability. We note that retraining our selection-function classifier using additional features is straightforward.

We publicly distribute the trained classifiers to encapsulate the sensitivity of DES and PS1 (see footnote 43). These classifiers can be used to predict the probability that a satellite would be detected as a function of its physical parameters and location on the sky. These detection probabilities can be used to predict the number of observed galaxies from a model of the full underlying Milky Way satellite galaxy population (e.g., Paper II).

## 8. Satellite Luminosity Function

The observational selection functions derived for DES and PS1 allow us to predict the detectability of a satellite, given its physical properties. As a simple demonstration, we apply these selection functions to estimate the luminosity function of Milky Way satellite galaxies assuming that satellites are distributed isotropically and that the physical properties of the observed

<footnote n="55">⁵⁵ We could similarly define a regression task, which would predict the detection significance directly, but we find that the output of a classification task is more easily interpreted when deriving population statistics.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Luminosity functions of Milky Way satellite galaxies within $R_{vir}$ = 300 kpc inferred from our analysis under the assumptions of an isotropic satellite distribution and a cored NFW radial distribution. (Left) Differential number of satellites within bins of absolute V-band magnitude. This calculation uses satellites detected with our selection threshold of $\sqrt{\text{TS}} > 6$ and SIG > 6. The observed satellites have been volume corrected using the observational selection functions for DES and PS1 derived here. Uncertainties are estimated from the Poisson error on the number of satellites in each bin. The power-law function $dN/dM_V = 10 \times 10^{0.1(M_V+5)}$ from Koposov et al. (2008) is shown as a dashed gray line. (Right) Cumulative number of Milky Way satellites. All known satellites are shown by the dashed gray line, while the satellites detected in our search of DES Y3A2 and PS1 DR1 are shown by the black dashed line. The solid black line shows detected satellites in DES and PS1 weighted by the observable volume correction for each satellite. Also shown are results from Newton et al. (2018) and Nadler et al. (2019b) that combine the observed satellite population with numerical simulations of dark matter and galaxy formation.</figcaption></figure>

satellites are representative of the total population. The results of this analysis may be compared to similar analyses based on data from SDSS (e.g., Koposov et al. 2008; Tollerud et al. 2008; Walsh et al. 2009); however, recent observations suggest that a more complex modeling framework is necessary (e.g., Jethwa et al. 2018; Newton et al. 2018; Nadler et al. 2019b). In Paper II, we perform a more rigorous analysis of the Milky Way satellite population based on cosmological simulations, which includes the effects of satellite disruption, galaxy formation efficiency, and the presence of the LMC.

Following a procedure similar to that of Koposov et al. (2008), we used the trained machine-learning classifier to derive the probability that each satellite will be detected as a function of Galactocentric position. For each confirmed and candidate dwarf galaxy (indexed by i) in the DES and PS1 footprints (class $\ge$3 in Table 2), we calculated the probability that a galaxy with the same parameters would be found within the survey volume,

$$C_{V,i} = \frac{\int \int_{R_{\text{min}}}^{R_{\text{max}}} P_{\text{det}}(D(\mathbf{r}), \rho_{*}(\mathbf{r}), M_{V,i}, r_{1/2,i}) n(r) r^2 dr \, d\Omega}{\int \int_{R_{\text{min}}}^{R_{\text{max}}} n(r) r^2 dr \, d\Omega}, \quad (3)$$

where $P_{\text{det}}$ is the probability that a satellite with absolute magnitude $M_{V,i}$ and physical half-light radius $r_{1/2,i}$ will be detected at a heliocentric distance D and local stellar density $\rho_*$, which are functions of the Galactocentric location of the satellite, $\mathbf{r}$.<fnref n="56" /> Following Koposov et al. (2008), we assumed that the Galactocentric radial distribution of satellites follows a cored NFW profile, $n(r) \propto r^{-2}(r + r_c)^{-1}$, where a core radius of $r_c = 20$ kpc was adopted to prevent divergence in the inner regions. Moreover, we adopted an inner cutoff of $R_{\text{min}} = 20$ kpc in our radial integration to crudely estimate the tidal effects of the Galactic disk, which generally disrupts subhalos that pass within this radius (e.g., Garrison-Kimmel et al. 2017). We integrate the satellite population out to $R_{\text{max}} = 300$ kpc, which is consistent with similar analyses in the literature (e.g., Hargis et al. 2014;

Jethwa et al. 2018; Newton et al. 2018; Nadler et al. 2019b) and is roughly comparable to the virial radius of the Milky Way (e.g., Garrison-Kimmel et al. 2014). The angular distribution of satellites is assumed to be isotropic in Galactocentric coordinates and is transformed into celestial equatorial coordinates to estimate the local stellar density.

To incorporate information from both surveys, we defined a weight for each satellite,

$$w_i = \frac{1}{C_{V,i,\text{DES}} + C_{V,i,\text{PS1}}}, \quad (4)$$

which combines the volumetric correction factors for the two surveys. To avoid double counting, we removed the area from PS1 that overlapped with DES.

The resulting volume-weighted Milky Way satellite luminosity function is shown in Figure 8. The left panel of Figure 8 shows the differential volume-weighted number of satellites in bins of absolute magnitude, $dN/dM_V$. Both PS1 and DES are largely complete out to the virial radius for satellites with $M_V \lesssim -7.5$ and $r_{1/2} \lesssim 200$ pc. Thus, for brighter satellites, the volume correction amounts to an area correction due to incomplete sky coverage. Uncertainties on the weighted number of satellites in each bin are calculated as $\sqrt{\sum_i w_i^2}$, which is equivalent to a weighted Poisson uncertainty.

For comparison, we also plot the differential luminosity function measured by Koposov et al. (2008) using data from SDSS DR5 (gray points with uncertainties). When compared to previous analyses of SDSS (Koposov et al. 2008; Walsh et al. 2009), our analysis covers ~3 times the sky area and includes ~3 times as many satellite galaxies. This allows us to extend the direct calculation of the differential Milky Way satellite luminosity function to both brighter and fainter systems. To guide the eye, we include the power-law model, $dN/dM_V = 10 \times 10^{0.1(M_V+5)}$, described by Koposov et al. (2008).

The right panel of Figure 8 shows the cumulative number of satellites brighter than a given absolute magnitude. We show both the total number of known satellites (gray dashed line) and the number of satellites detected in our search of DES and PS1 (black dashed line). We correct the detected satellite curve via the volumetric weights described previously to yield the volume-corrected cumulative satellite luminosity function

<footnote n="56">⁵⁶ We assume the Galactocentric distance of the Sun is 8.3 kpc (Gillessen et al. 2009).</footnote>
(black solid line). This estimate of the satellite luminosity function assumes that satellites are distributed isotropically, that the radial distribution of satellites goes as $r^{-3}$ at large radii, and that the physical sizes of the known satellites are representative of the population as a whole. This analysis predicts that the Milky Way contains ~270 satellite galaxies within 300 kpc with $M_V \lesssim 0$. The resulting cumulative number of satellites is slightly higher than predictions from some cosmological simulations (e.g., Newton et al. 2018; Nadler et al. 2019b) but lower than others (e.g., Kelley et al. 2019). The sharp upturn in $N(<M_V)$ around $M_V \sim -4$ comes from recent discoveries in the Southern Hemisphere, combined with increasing volumetric weights at these faint magnitudes. This upturn at faint magnitudes has been attributed to the existence of a population of faint satellites associated with the LMC (Jethwa et al. 2018; Newton et al. 2018) and/or reionization physics (Bose et al. 2018). We discuss the importance of modeling the LMC satellite system in more detail in Paper II.

The analysis described herein offers few explicit modeling choices (i.e., the functional form of $n(r)$); however, it contains several implicit assumptions about the satellite population. First, we have assumed that the physical properties of the undetected satellite population are consistent with the population of detected satellites. Such an assumption breaks down if a large population of extremely low-surface-brightness systems exists (e.g., satellites like Crater II and Antlia II). A population of large, luminous, and distant Milky Way satellites has been suggested to resolve the observed radial distribution discrepancy between the satellite populations of the Milky Way and Andromeda (e.g., Yniguez et al. 2014; Samuel et al. 2020). In addition, we assumed that the intrinsic properties of satellites (e.g., luminosity and half-light radius) are independent of heliocentric distance. Again, this assumption may break down due to the larger influence of baryonic physics on the properties of satellites close to the Milky Way (e.g., Nadler et al. 2018). Finally, we made the simplifying assumption that the satellite population is isotropic, while observations suggest that this is very unlikely to be true for the Milky Way (e.g., Pawlowski et al. 2012; Bechtol et al. 2015; Drlica-Wagner et al. 2015; Koposov et al. 2015) and Andromeda (e.g., Ibata et al. 2013). These assumptions strongly motivate the more rigorous analysis in Paper II, which simultaneously includes the effects of satellite disruption (enhanced by the presence of the Galactic disk), uncertainties in the efficiency of galaxy formation in low-mass halos, and the influence of the LMC.

## 9. Conclusions

We have conducted a search for low-luminosity Milky Way satellite galaxies over ~24,000 deg² (~75% of the high-Galactic-latitude sky) using data from DES Y3A2 and PS1 DR1. We recover most of the satellites previously discovered in the DES and PS1 footprints, though a number of low-luminosity satellites are beyond the sensitivity limit of our search. We set a detection threshold intended to minimize the number of false-positive candidates, and no new, high-confidence satellite galaxy candidates were discovered. The only significant new candidate is a compact, low-luminosity outer-halo star cluster (Appendix D).

We determined the sensitivity of our search by simulating the resolved stellar populations of Milky Way satellite galaxies, injecting them into the DES and PS1 catalogs, and running the same satellite detection algorithms that were applied to the real data. We quantified the observational sensitivity of our search in

terms of satellite properties (i.e., absolute magnitude, physical size, distance, and local stellar density) and provide both analytic and machine-learning models of the observational selection function. Finally, we demonstrated the application of our observational selection function to derive a data-driven luminosity function from the observed Milky Way satellite galaxy population. By encapsulating the observational selection function in a flexible and accurate machine-learning model, we facilitate more rigorous statistical-inference-based approaches to extract the properties of galaxy formation and dark matter physics.

Deep, multi-band optical imaging over the entire sky at depths fainter than $g \sim 23$ is now within reach. DES covers only one-sixth of the sky accessible to DECam, and the past several years have seen an active community campaign to complete contiguous DECam coverage of the entire southern sky. Programs like MagLiteS (Drlica-Wagner et al. 2016; Torrealba et al. 2018), BLISS (Mau et al. 2019), DECaLS (Dey et al. 2019), and DELVE (Mau et al. 2020)<fnref n="57" /> are actively collecting, processing, and mining these data for fainter and more distant satellites. Meanwhile, HSC SSP on the 8.2 m Subaru telescope will achieve $r \sim 26$ over ~1400 deg², thereby extending the search for the faintest galaxies to unprecedented distances (Homma et al. 2016, 2018, 2019). In the early 2020s, the Legacy Survey of Space and Time (LSST) conducted from the Vera C. Rubin Observatory will expand this depth of coverage to the entire southern sky. The power of these upcoming surveys, combined with advanced modeling techniques, promise to shed new light on the darkest galaxies.

This research was supported in part by the National Science Foundation (NSF) under Grant No. NSF PHY-1748958 through the Kavli Institute for Theoretical Physics program “The Small-Scale Structure of Cold(?) Dark Matter,” and grant no. NSF DGE-1656518 through the NSF Graduate Research Fellowship received by E.O.N. This research made use of computational resources at SLAC National Accelerator Laboratory, a U.S. Department of Energy Office; the authors thank the support of the SLAC computational team. This research has made use of the WEBDA database, operated at the Department of Theoretical Physics and Astrophysics of the Masaryk University.

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

<footnote n="57">⁵⁷ https://delve-survey.github.io</footnote>
Chicago, University College London, the DES-Brazil Con-
sortium, the University of Edinburgh, the Eidgenössische
Technische Hochschule (ETH) Zürich, Fermi National Accel-
erator 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 NSFs National Optical-Infrared
Astronomy Laboratory, 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.

This work is based in part on observations at Cerro Tololo
Inter-American Observatory, NSFs National Optical-Infra-
red Astronomy Laboratory, 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 Frame-
work Program (FP7/2007-2013), including ERC grant agree-
ments 240672, 291329, and 306478. We acknowledge support
from 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 U.S. Government retains—and the
publisher, by accepting the article for publication, acknowl-
edges that the U.S. government retains—a non-exclusive, paid-
up, irrevocable, worldwide license to publish or reproduce the
published form of this manuscript, or allow others to do so, for
U.S. government purposes.

# Appendix A
## Data Selection

A.1. DES Y3A2

The DES data were selected from the Y3A2 internal release of
the GOLD catalog (v2.0) accessed via a bulk download from the
easyaccess SQL command line interpreter (Carrasco Kind
et al. 2018). We removed objects with FLAGS_GOLD &
0b111100, which identifies SOF fit failures, objects with
SExtractor FLAGS > 3, objects with bad pixels in their isophotes
IMAFLAGS_ISO != 0, objects characterized as bright blue
color outliers (Section 6.2 of Drlica-Wagner et al. 2018), and
other extreme color outliers (Section 6.2 of Drlica-Wagner
et al. 2018). Magnitudes were corrected for a CCD-dependent
magnitude adjustment to the zero-points (DELTA MAG_V3), an
SED-dependent chromatic correction (DELTA MAG_CHROM),
and extinction using the extinction maps from Schlegel
et al. (1998; A_SED_SFD98). Star-galaxy classification used

the EXTENDED_CLASS_MASH_SOF variable, with stars selected
with 0 ≤ EXTENDED_CLASS_MASH_SOF ≤ 2. The construc-
tion of the extended classification variable is described in Section
2 of Shipp et al. (2018). Briefly, this classifier uses the SOF size
parameter CM_T (and associated error) to classify most objects.
When the SOF parameters are unavailable (due to a small fraction
of objects where the SOF fit fails), the weighted average of the
single-epoch measurements WAVG_SPREAD MODEL (and asso-
ciated error) are used for classification. For faint objects where
both SOF and WAVG values are unavailable, classification is
performed with the SPREAD MODEL (and associated uncertainty)
parameters derived from the coadded images.

A.2. PS1 DR1

The PS1 DR1 data were downloaded from the MAST
CasJobs server using queries similar to the example provided as
follows. Note that the query selects objects in an interval in
objid. The full bulk download is ~10,000 slices in objid,
designed to have approximately equal numbers of objects per
slice. Each slice is effectively a narrow interval in decl. In a
post-processing step, the objects were partitioned into a
HEALPix grid for rapid access during the search phase.

```sql
SELECT
  f.raMean, f.decMean,
  f.objID, f.uniquePspsoBID,
  f.objInfoFlag, f.qualityFlag,
  f.nStackDetections, f.nDetections, f.ng, f.nr,
  f.ni,
  f.gFPSFFlux, f.rFPSFFlux, f.iFPSFFlux,
  f.gFPSFFluxErr, f.rFPSFFluxErr, f.iFPSFFluxErr,
  f.gFKronFlux, f.rFKronFlux, f.iFKronFlux,
  f.gFKronFluxErr, f.rFKronFluxErr, f.
  iFKronFluxErr,
  f.gFlags, f.rFlags, f.iFlags,
  s.gInfoFlag, s.rInfoFlag, s.iInfoFlag,
  s.gInfoFlag2, s.rInfoFlag2, s.iInfoFlag2,
  s.gInfoFlag3, s.rInfoFlag3, s.iInfoFlag3,
  s.primaryDetection, s.bestDetection
INTO MyDB.ps1_dr1_objid_bin_09991_10000_v02
FROM ForcedMeanObjectView f
JOIN StackObjectThin s on f.objid = s.objid
WHERE f.objid BETWEEN 212390000000001521 AND
2125900000000007936
AND f.qualityFlag & 16 > 0;
```

We remove duplicate objects by requiring

```
sel=(data["primaryDetection"==1])
```

and apply additional selection criteria to filter the high-
quality catalog objects.

```
sel &= (data["nStackDetections" > 1]
  && (data["nDetections" > 0]
  && (data["gInfoFlag" >= 0]
  && (data["rInfoFlag" >= 0]
  && (data["iInfoFlag" >= 0]
  && ((data["gInfoFlag"] & (8 + 2048)) == 0)
  && ((data["rInfoFlag"] & (8 + 2048)) == 0)
  && ((data["iInfoFlag"] & (8 + 2048)) == 0)
  && ((data["gInfoFlag2"] & (8192 + 4194304)) == 0)
  && ((data["rInfoFlag2"] & (8192 + 4194304)) == 0)
  && ((data["iInfoFlag2"] & (8192 + 4194304)) == 0)
```
These criteria were determined by comparing object detections in PS1 to detections in deeper imaging data from DES and HSC, to identify the quantities that correlate with higher rates of unmatched PS1 objects.

Our star-galaxy selection was performed using the measured aperture and PSF magnitudes in i-band:

$$ 
\begin{aligned}
\texttt{sel} &= ((\texttt{data["iFPSFMag"]} - \texttt{data["IFKronMag"]}) < 0.05) \\
\texttt{sel} &= ((\texttt{data["iFPSFMag"]} == -999.) \\
\texttt{sel} &= ((\texttt{data["iFPSFMag"]} - \texttt{data["IFKronMag"]}) > 4.0)
\end{aligned}
$$

The i-band was chosen due to its superior PSF and depth.

For the simple analysis, we also require moderate signal-to-noise object detections (S/N > 10) to ensure a sample of high-confidence stars:

$$ \texttt{sel} = ((\texttt{data["rFPSFMagErr"]} < 0.1) $$

While this S/N selection limited the number of faint stars used in the search (Figure 3), it significantly reduced the number of spurious candidates returned by the simple algorithm. The S/N selection was not applied to the ugali search on PS1.

# Appendix B

## Validation of the Parameter Space for Simulated Satellites

The ability to detect satellites depends largely on the surface brightness and number of resolved member stars recovered by a survey (e.g., Walsh et al. 2009). These observed properties depend both on the physical characteristics of the satellite

(e.g., luminosity, size, distance) and the characteristics of the survey (e.g., completeness, coverage). For DES and PS1, we parameterize these two attributes using the average surface brightness within the half-light radius, $\mu$, and the number of satellite member stars brighter than $g = 22$, $N(g < 22)$. The $N(g < 22)$ parameter is challenging to compare directly to models of galaxy formation due to its dependence on the survey completeness, but is readily obtained from the numerical simulations of satellites, and is an effective means to gauge whether a given satellite could be detected.

To assess the sensitivity of our two search algorithms (simple and ugali) applied to our two data sets (DES Y3A2 and PS1 DR1), we show the recovered significance of satellites as a function of their surface brightness, $\mu$, and the number of bright resolved member stars, $N(g < 22)$, in Figure B1. The detectability of satellites is easily characterized in this parameter space; in particular, for DES Y3A2 (PS1 DR1), satellites brighter than $\mu \sim 30$ mag arcsec$^{-2}$ ($\mu \sim 28$ mag arcsec$^{-2}$) with more than $N(g < 22) \gtrsim 10^2$ resolved member stars are reliably detected, while lower surface brightness systems or those with fewer resolved stars are detected less efficiently. The cross-hatching indicates regions where satellites start to become large enough that sky projection effects come into consideration ($a_h > 5$ deg; DIFFICULTY = 1) or have so many stars that direct simulation becomes too computationally intensive ($N(g < 22) > 10^3$ and $\mu < 23.5$ mag arcsec$^{-2}$; DIFFICULTY = 2). We only used simulations with DIFFICULTY = 0 when training the machine-learning model, though the classifier accurately predicts the detectability of simulated satellites in regions of parameter space with DIFFICULTY != 0.
<figure><img src="image_10.png" /><figcaption>Figure B1. Detection significance of simulated satellites as a function of surface brightness, $\mu$, and the number of member stars brighter than $g = 22$, $N(g < 22)$. The top (bottom) row shows results for DES Y3A2 (PS1 DR1) with each point corresponding to a simulated satellite. The colorbar indicates the detection significance corresponding to the simple SIG value (left column) or the ugali $\sqrt{\text{TS}}$ (right column). Satellites with $a_h > 5$ deg are flagged with DIFFICULTY = 1 and are not used in our training. Satellites in the region labeled DIFFICULTY = 2 have both high surface brightness and many resolved stars. We assume satellites in this region are always detected and are therefore not simulated for computational efficiency.</figcaption></figure>

## Appendix C
### Maximum-likelihood Formalism

The extended maximum-likelihood formalism (e.g., Orear 1958; Barlow 1990) is widely used in physics and astronomy as a means to perform statistical searches for overdensities in a Poisson random field with unknown signal normalization (e.g., galaxy cluster searches as described in Kepner et al. 1999 or Rykoff et al. 2012). In this appendix, we describe the application of this formalism to the search for stellar overdensities, as implemented in the ugali software package.<fnref n="58" />

We assume that the angular and magnitude distribution of stars in a survey can be described by a Poisson realization of the underlying stellar density field. Dividing the survey into

bins, we can write the binned Poisson likelihood, $\mathcal{L}$, as

$$ \mathcal{L} = \prod_{i}^{\text{bins}} P_i = \prod_{i}^{\text{bins}} \frac{n_i k_i e^{-n_i}}{k_i!}, \qquad (C1) $$

where in bin $i$, $P_i$ is the Poisson probability to observe $k_i$ stars, given a model expectation of $n_i$ stars. This likelihood is valid for binning over any domain, but for our specific application to satellite galaxy searches, we commonly utilize two spatial dimensions (e.g., R.A. and decl.) and two magnitude dimensions (e.g., $g$ and $r$). For numerical simplicity, we generally work with the logarithm of the likelihood

$$ \log \mathcal{L} = \sum_{i}^{\text{bins}} \{-n_i + k_i \log(n_i) - \log(k_i!}\}. \qquad (C2) $$

The last term in this equation, $\log(k_i!)$, does not depend on any of the model parameters and can be safely discarded as an additive constant.

<footnote n="58">58 https://github.com/DarkEnergySurvey/ugali</footnote>
The number of model-predicted stars in a bin is the integral of the probability density function (PDF) of the model over that bin, $n_i = \int_i m(\lambda, \theta)dV$. In this equation, the integral is over the observable “volume” of the bin, while $\lambda$ and $\theta$ are parameters of our model. We explicitly formulate our model for the stellar counts in terms of a contribution from the satellite and a contribution from Milky Way field stars,

$$m(\lambda, \theta) = \lambda s(\theta) + b(\theta). \qquad (C3)$$

In the preceding equations, we have been explicit about the dependence of $m$ on the signal normalization parameter, $\lambda$, which we call the “richness.”<fnref n="59" />

In the limit of infinitesimally small bins, each bin contains either zero or one star ($k_i \in \{0, 1\}$), and the log-likelihood can be expressed as

$$\log \mathcal{L} = - \sum_{i}^{\text{empty}} n_i - \sum_{i}^{\text{filled}} n_i + \sum_{i}^{\text{filled}} \log(n_i). \qquad (C4)$$

By definition, the sum over empty and filled bins is equivalent to the sum over all bins. Thus, the first two terms can be written as the integral of the model over all observable space,

$$\sum_{i}^{\text{bins}} n_i = \int_{V} m(\lambda, \theta) dV \qquad (C5)$$

$$= \lambda N_s + N_b. \qquad (C6)$$

In the last line, we have used $N_s$ and $N_b$ to represent the total number of expected satellite stars and Milky Way field stars, respectively. The only filled bins will be located on the observed stars. Thus, the last term in Equation (C4) can be rewritten as a sum over stars, where we replace the model-predicted counts with the model PDF evaluated at the location of each star,

$$\sum_{i}^{\text{filled}} \log(n_i) = \sum_{j}^{\text{stars}} \log(\lambda s_j + b_j), \qquad (C7)$$

where $j$ indexes over stars, and we use $s_j$, $b_j$ as a shorthand for the components of the model PDF evaluated at the location of each star. Substituting back into Equation (C4) yields

$$\log \mathcal{L} = - \lambda N_s - N_b + \sum_{j}^{\text{stars}} \log(\lambda s_j + b_j), \qquad (C8)$$

which is conventionally referred to as the unbinned Poisson log-likelihood function.

When searching for stellar overdensities, we are interested in maximizing the likelihood with respect to the richness parameter, $\lambda$. Differentiating Equation (C8) with respect to $\lambda$ gives

$$\frac{d \log \mathcal{L}}{d \lambda} = -N_s + \sum_{j}^{\text{stars}} \frac{s_j}{\lambda s_j + b_j}. \qquad (C9)$$

The maximum-likelihood estimator for the richness, $\lambda = \hat{\lambda}$, will occur when the derivative of the likelihood is zero; thus,

$$\hat{\lambda} N_s = \sum_{j}^{\text{stars}} \frac{\hat{\lambda} s_j}{\hat{\lambda} s_j + b_j}. \qquad (C10)$$

The left-hand side of this equation represents the number of observable stars predicted by the signal model, while the right-hand side represents the sum of the probabilities that each star belongs to the signal distribution,

$$p_j \equiv \frac{\lambda s_j}{\lambda s_j + b_j}. \qquad (C11)$$

If we care only about the dependence of the likelihood on $\lambda$, we can rewrite Equation (C8) in terms of the observable fraction and the membership probabilities,<fnref n="60" />

$$\log \mathcal{L} = - \lambda N_s - \sum_{j}^{\text{stars}} \log(1 - p_j). \qquad (C12)$$

This formulation makes explicit the dependence of the likelihood on the membership probabilities, but discards terms that depend on the background model alone.

From Equation (C3), it is clear that there is a degeneracy between the normalization of the signal PDF, $s$, and the richness, $\lambda$. Some authors (e.g., Martin et al. 2008) choose to normalize the signal PDF to unity over the observed space, thus requiring an observational correction to predict the number of model stars below the detection threshold of a survey. In contrast, the ugali framework normalizes the signal PDF to unity over the entire observable domain, $\int_{\text{all}} s(\theta)dV = 1$, and defines the observable fraction of the signal, $f$, at a given location in a survey. Thus, the number of observable stars is

$$\lambda N_s = \lambda \int_{\text{all}} f s(\theta) dV = \lambda f. \qquad (C13)$$

This definition allows us to interpret $\lambda$ as the total number of stars in the satellite, rather than the number of observed stars.

## Appendix D: Remaining Candidates

After applying the first three selection criteria described in Section 6.1, we are left with a list of 28 candidates in the PS1 DR1 footprint. We visually inspected each of these candidates to determine whether they were viable new satellite galaxies. We find that most of the remaining candidates could be clearly identified as artifacts in the PS1 DR1 catalogs. The most obvious artifact manifested as rectangular regions of increased or decreased object density (often affecting the density of both stars and galaxies) with characteristic sizes of ~0.4 × 0.4 deg (Figure D1). These regions correspond to the skycells over which the PS1 DR1 stack images were created (Flewelling et al. 2016). We visually inspected these skycells and found that they contained bright stars that were causing issues in the automated image processing and catalog creation. The corresponding under/overdensities in catalog objects bias estimates of the field density and lead to spurious candidate detection. We masked regions of radius 0.5 deg around each of the seven skycells identified in this manner. In addition, we find eight cases where reflected and scattered light from bright stars leads to spurious overdensities of blue objects. These artifacts are less extreme than the skycell failures described previously, and we mask a circular region with a radius matched to the diagonal dimension of a skycell (~0.3 deg). We visually inspected the remaining 13 candidates and found that all but one of them lacked

<footnote n="59">⁵⁹ The term “richness” is inherited from the analysis of galaxy clusters (e.g., Abell 1958).</footnote>

<footnote n="60">⁶⁰ We note a typo in Bechtol et al. (2015) that omitted the “log” from the right-hand side of Equation (C12).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure D1.** Artifact resulting from the bright star (BD+18 240) biasing the processing and catalog creation for a stack image in a single PS1 DR1 skycell. Artifacts like this were identified visually in the PS1 DR1 candidate list, and regions around the affected skycells were masked. Similar figure to Figure D2, but with the stellar and galactic densities convolved with a Gaussian kernel of 0.88.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure D2.** DELVE 1 is a candidate faint outer-halo star cluster discovered in our search of PS1 DR1 data and independently in DELVE (Mau et al. 2020). The small physical size ($r_{1/2}$ = 3.7 pc) and low luminosity ($M_V$ = 0.6) make this system similar to other outer-halo star systems that have been discovered in recent surveys (e.g., Mau et al. 2019; Torrealba et al. 2019b). (Left) Stellar density convolved with a Gaussian kernel of 0.67. (Center) Galaxy density convolved with a Gaussian kernel of 0.67. (Right) Hess diagram corresponding to foreground stars within $3r_h$ of the centroid and background stars in a concentric annulus from $5r_h$ to $5.8r_h$. The best-fit isochrone is shown in black.</figcaption></figure>

a well-defined stellar sequence in color–magnitude space and/or a well-defined spatial morphology. These regions likely result from a combination of less obvious survey artifacts (i.e., mis-estimation of the sky background, excess sensor noise, or low level scattered light) and contamination from background galaxy clusters. We mask regions of ~0.3 deg around each of these seeds. These masks were applied to our candidate list and were included in the geometric masking cuts described in Section 6.1.

As discussed in Section 6, one candidate passed all of our detection criteria (Figure D2). Visual inspection of this candidate did not identify any obvious survey artifacts that would significantly alter the stellar density in this region. Similarly, we found no evidence of a correlated overdensity in the galaxy sample. This candidate appeared similar to other recently discovered compact outer-halo star clusters (e.g., Mau et al. 2019; Torrealba et al. 2019b), and we followed the procedure of Bechtol et al. (2015) and Drlica-Wagner et al. (2015) to derive the physical parameters of this candidate using the maximum- likelihood fitting formalism of ugali. We simultaneously fit the richness, centroid position, angular extension, ellipticity, position angle, and distance modulus of this system using an affine invariant Markov Chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013). Best-fit parameters and uncertainties were derived from the marginalized posterior distribution (sampled with $2.5 \times 10^5$ steps) and are reported in Table D1. During this sampling, the age and metallicity were

<table><caption>**Table D1**
Observed and Derived Properties of DELVE 1</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>α<sub>2000</sub></td><td>247.725<sup>+0.004</sup><sub>-0.004</sub></td><td>degree</td></tr><tr><td>δ<sub>2000</sub></td><td>-0.971<sup>+0.005</sup><sub>-0.004</sub></td><td>degree</td></tr><tr><td>a<sub>h</sub></td><td>0.8<sup>+0.5</sup><sub>-0.3</sub></td><td>arcmin</td></tr><tr><td>r<sub>h</sub></td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>arcmin</td></tr><tr><td>r<sub>1/2</sub></td><td>3.7<sup>+1</sup><sub>-1</sub></td><td>pc</td></tr><tr><td>ε</td><td>0.04</td><td>...</td></tr><tr><td>P.A.</td><td>25<sup>+71</sup><sub>-62</sub></td><td>degree</td></tr><tr><td>m - M<sup>a</sup></td><td>16.0<sup>+0.4</sup><sub>-0.2</sub> ± 0.1</td><td>...</td></tr><tr><td>D<sub>☉</sub></td><td>15.6<sup>+3.5</sup><sub>-1.1</sub></td><td>kpc</td></tr><tr><td>Σp<sub>t</sub></td><td>29<sup>+5</sup><sub>-5</sub></td><td>...</td></tr><tr><td>TS</td><td>69</td><td>...</td></tr><tr><td>M<sub>V</sub><sup>b</sup></td><td>0.6<sup>+0.7</sup><sub>-1.9</sub></td><td>mag</td></tr><tr><td>M<sub>*</sub></td><td>86<sup>+40</sup><sub>-31</sub></td><td>M<sub>☉</sub></td></tr><tr><td>μ</td><td>26.4</td><td>mag arcsec<sup>-2</sup></td></tr></tbody></table>

**Notes.** Uncertainties were derived from the highest density interval containing the peak and 68% of the marginalized posterior distribution.

<sup>a</sup> We assume a systematic uncertainty of ±0.1 associated with isochrone modeling.

<sup>b</sup> The uncertainty in $M_V$ was calculated following Martin et al. (2008) and does not include uncertainty in the distance.
fixed at $\tau = 10$ Gyr and $Z = 0.0001$, respectively. We conclude that this candidate is likely a compact ($r_{1/2} = 3.7$ pc), low-luminosity ($M_V = 0.6$) star cluster residing at a heliocentric distance of $D = 15.6$ kpc. We note that this system was independently discovered in deeper data from the DECam Local Volume Exploration survey (DELVE), where it was investigated in more detail and named DELVE 1 (Mau et al. 2020). We use the same name here to avoid confusion.

## Appendix E

### Machine-learning Classifier Parameters

We modeled the observational selection function with a gradient-boosted decision tree classifier, XGBClassifier, as implemented in Python in the `xgboost` package version 0.82 (Chen & Guestrin 2016).<sup>61</sup> We trained separate classifiers on the DES and PS1 simulations. Hyperparameters were selected using GridSearchCV from `scikit-learn` version 0.19.1 (Pedregosa et al. 2011). The hyperparameters scanned are `learning_rate` = {0.01, 0.05, 0.1}, `max_depth` = {6, 7, 8}, and `n_estimators` = {100, 250, 500}, with fixed `max_delta_step` = 1. The optimal training hyperparameters for the DES simulations were `learning_rate` = 0.01, `max_depth` = 8, and `n_estimators` = 500. The optimal training hyperparameters for the PS1 simulations were `learning_rate` = 0.05, `max_depth` = 7, and `n_estimators` = 250. The trained classifiers are available at https://github.com/des-science/mw-sats.

### ORCID iDs

A. Drlica-Wagner https://orcid.org/0000-0001-8251-933X
K. Bechtol https://orcid.org/0000-0001-8156-0429
S. Mau https://orcid.org/0000-0003-3519-4004
M. McNanna https://orcid.org/0000-0001-5435-7820
E. O. Nadler https://orcid.org/0000-0002-1182-3825
A. B. Pace https://orcid.org/0000-0002-6021-8760
T. S. Li https://orcid.org/0000-0002-9110-6163
A. Pieres https://orcid.org/0000-0001-9186-6042
E. Rozo https://orcid.org/0000-0002-1666-6275
A. R. Walker https://orcid.org/0000-0002-7123-8943
R. H. Wechsler https://orcid.org/0000-0003-2229-011X
S. Allam https://orcid.org/0000-0002-7069-7857
J. Annis https://orcid.org/0000-0002-0609-3987
D. Brooks https://orcid.org/0000-0002-8458-5047
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
J. De Vicente https://orcid.org/0000-0001-8318-6813
S. Desai https://orcid.org/0000-0002-0466-3288
H. T. Diehl https://orcid.org/0000-0002-8357-7467
T. F. Eifler https://orcid.org/0000-0002-1894-3301
B. Flaugher https://orcid.org/0000-0002-2367-5049
J. Frieman https://orcid.org/0000-0003-4079-3263
J. García-Bellido https://orcid.org/0000-0002-9370-8360
E. Gaztanaga https://orcid.org/0000-0001-9632-0815
D. Gruen https://orcid.org/0000-0003-3270-7644
R. A. Gruendl https://orcid.org/0000-0002-4588-6517
J. Gschwend https://orcid.org/0000-0003-3023-8362
G. Gutierrez https://orcid.org/0000-0003-0825-0517
K. Honscheid https://orcid.org/0000-0002-6550-2023

D. J. James https://orcid.org/0000-0001-5160-4486
E. Krause https://orcid.org/0000-0001-8356-2014
K. Kuehn https://orcid.org/0000-0003-0120-0808
N. Kuropatkin https://orcid.org/0000-0003-2511-0946
O. Lahav https://orcid.org/0000-0002-1134-9035
M. A. G. Maia https://orcid.org/0000-0001-9856-9307
J. L. Marshall https://orcid.org/0000-0003-0710-9474
P. Melchior https://orcid.org/0000-0002-8873-5065
F. Menanteau https://orcid.org/0000-0002-1372-2534
R. Miquel https://orcid.org/0000-0002-6610-4836
A. Palmese https://orcid.org/0000-0002-6011-0530
A. A. Plazas https://orcid.org/0000-0002-2598-0514
E. Sanchez https://orcid.org/0000-0002-9646-8198
M. Schubnell https://orcid.org/0000-0001-9504-2059
I. Sevilla-Noarbe https://orcid.org/0000-0002-1831-1953
M. Smith https://orcid.org/0000-0002-3321-1432
E. Suchyta https://orcid.org/0000-0002-7047-9358
G. Tarle https://orcid.org/0000-0003-1704-0781

### References

Abell, G. O. 1958, ApJS, 3, 211
Ackermann, M., Albert, A., Anderson, B., et al. 2015, PhRvL, 115, 231301
Agertz, O., Pontzen, A., Read, J. I., et al. 2020, MNRAS, 491, 1656
Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8
Balbinot, E., Yanny, B., Li, T. S., et al. 2016, ApJ, 820, 58
Barlow, R. 1990, NIMPA, 297, 496
Bechtol, K., Drlica-Wagner, A., Balbinot, E., et al. 2015, ApJ, 807, 50
Belokurov, V., Irwin, M. J., Koposov, S. E., et al. 2014, MNRAS, 441, 2124
Belokurov, V., Walker, M. G., Evans, N. W., et al. 2008, ApJL, 686, L83
Belokurov, V., Walker, M. G., Evans, N. W., et al. 2009, MNRAS, 397, 1748
Belokurov, V., Walker, M. G., Evans, N. W., et al. 2010, ApJL, 712, L103
Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJL, 647, L111
Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007, ApJ, 654, 897
Bergström, L., Ullio, P., & Buckley, J. H. 1998, APh, 9, 137
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
Bica, E., Bonatto, C., Dutra, C. M., & Santos, J. F. C. 2008, MNRAS, 389, 678
Bose, S., Deason, A. J., Belokurov, V., & Frenk, C. S. 2019, arXiv:1909.04039
Bose, S., Deason, A. J., & Frenk, C. S. 2018, ApJ, 863, 123
Boylan-Kolchin, M., Weisz, D. R., Johnson, B. D., et al. 2015, MNRAS, 453, 1503
Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
Brown, T. M., Tumlinson, J., Geha, M., et al. 2014, ApJ, 796, 91
Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343
Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517
Burke, D. L., Rykoff, E. S., Allam, S., et al. 2018, AJ, 155, 41
Cannon, R. D., Hawarden, T. G., & Tritton, S. B. 1977, MNRAS, 180, 81P
Carlin, J. L., Sand, D. J., Muñoz, R. R., et al. 2017, AJ, 154, 267
Carlin, J. L., Sand, D. J., Price, P., et al. 2016, ApJL, 828, L5
Carrasco Kind, M., Drlica-Wagner, A., Koziol, A., & Petravick, D. 2018, JOSS, 4, 1022
Chabrier, G. 2001, ApJ, 554, 1274
Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560
Chen, T., & Guestrin, C. 2016, in Proc. 22nd ACM SIGKDD Int. Conf. on Knowledge Discovery and Data Mining, KDD '16 (New York, NY: ACM), 785
Chernoff, H. 1954, Ann. Math. Statist., 25, 573
Clesse, S., & García-Bellido, J. 2018, PDU, 22, 137
Conn, B. C., Jerjen, H., Kim, D., & Schirmer, M. 2018, ApJ, 852, 68
Corwin, H. G. 2004, yCat, 7239, 0
Crnojević, D., Sand, D. J., Bennet, P., et al. 2019, ApJ, 872, 80
DES Collaboration, Abbott, T. M. C., Abdalla, F. B., et al. 2018, ApJS, 239, 18
Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
Diehl, H. T., Neilsen, E., Gruendl, R., et al. 2016, Proc. SPIE, 9910, 99101D
Drlica-Wagner, A., Bechtol, K., Allam, S., et al. 2016, ApJL, 833, L5
Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al. 2015, ApJ, 813, 109
Drlica-Wagner, A., Sevilla-Noarbe, I., Rykoff, E. S., et al. 2018, ApJS, 235, 33
Fadely, R., Hogg, D. W., & Willman, B. 2012, ApJ, 760, 15
Fitzpatrick, E. L. 1999, PASP, 111, 63
Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150
Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2016, arXiv:1612.05243

<footnote n="61">61 https://xgboost.readthedocs.io</footnote>
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Garrison-Kimmel, S., Boylan-Kolchin, M., Bullock, J. S., & Lee, K. 2014, MNRAS, 438, 2578

Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2019, MNRAS, 487, 1380

Garrison-Kimmel, S., Wetzel, A., Bullock, J. S., et al. 2017, MNRAS, 471, 1709

Geringer-Sameth, A., Koushiappas, S. M., & Walker, M. G. 2015, PhRvD, 91, 083535

Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075

Graczyk, D., Pietrzyński, G., Thompson, I. B., et al. 2014, ApJ, 780, 59

Graus, A. S., Bullock, J. S., Kelley, T., et al. 2019, MNRAS, 488, 4585

Grillmair, C. J. 2006, ApJL, 645, L37

Grillmair, C. J. 2009, ApJ, 693, 1118

Hamilton, A. J. S., & Tegmark, M. 2004, MNRAS, 349, 115

Hargis, J. R., Willman, B., & Peter, A. H. G. 2014, ApJL, 795, L13

Harrington, R. G., & Wilson, A. G. 1950, PASP, 62, 118

Harris, W. E. 1996, AJ, 112, 1487

Hoffleit, D., & Jaschek, C. 1991, The Bright Star Catalogue (New Haven, CT: Yale Univ. Observatory)

Homma, D., Chiba, M., Komiyama, Y., et al. 2019, PASJ, 71, 94

Homma, D., Chiba, M., Okamoto, S., et al. 2016, ApJ, 832, 21

Homma, D., Chiba, M., Okamoto, S., et al. 2018, PASJ, 70, S18

Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Natur, 370, 194

Ibata, R. A., Lewis, G. F., Conn, A. R., et al. 2013, Natur, 493, 62

Irwin, M. J., Belokurov, V., Evans, N. W., et al. 2017, ApJL, 656, L13

Irwin, M. J., Bunclark, P. S., Bridgeland, M. T., & McMahon, R. G. 1990, MNRAS, 244, 16P

Ishiyama, T., Sudo, K., Yokoi, S., et al. 2016, ApJ, 826, 9

Jethwa, P., Erkal, D., & Belokurov, V. 2018, MNRAS, 473, 2060

Katz, H., Ramsoy, M., Rosdahl, J., et al. 2019, arXiv:1905.11414

Kelley, T., Bullock, J. S., Garrison-Kimmel, S., et al. 2019, MNRAS, 487, 4409

Kepner, J., Fan, X., Bahcall, N., et al. 1999, ApJ, 517, 78

Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R.-D. 2013, A&A, 558, A53

Kim, D., & Jerjen, H. 2015a, ApJ, 799, 73

Kim, D., & Jerjen, H. 2015b, ApJL, 808, L39

Kim, D., Jerjen, H., Geha, M., et al. 2016, ApJ, 833, 16

Kim, D., Jerjen, H., Mackey, D., Da Costa, G. S., & Milone, A. P. 2015a, ApJL, 804, L44

Kim, D., Jerjen, H., Milone, A. P., Mackey, D., & Da Costa, G. S. 2015b, ApJ, 803, 63

Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, PhRvL, 121, 211302

Kirby, E. N., Simon, J. D., & Cohen, J. G. 2015, ApJ, 810, 56

Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82

Knuth, D. E. 1973, The Art of Computer Programming. Vol. 3: Sorting and Searching (Reading, MA: Addison-Wesley)

Koposov, S., Belokurov, V., Evans, N. W., et al. 2008, ApJ, 686, 279

Koposov, S. E., Belokurov, V., Torrealba, G., & Evans, N. W. 2015, ApJ, 805, 130

Koposov, S. E., Irwin, M., Belokurov, V., et al. 2014, MNRAS, 442, L85

Koposov, S. E., Walker, M. G., Belokurov, V., et al. 2018, MNRAS, 479, 5343

Laevens, B. P. M., Martin, N. F., Bernard, E. J., et al. 2015a, ApJ, 813, 44

Laevens, B. P. M., Martin, N. F., Ibata, R. A., et al. 2015b, ApJL, 802, L18

Laevens, B. P. M., Martin, N. F., Sesar, B., et al. 2014, ApJL, 786, L3

Li, T. S., DePoy, D. L., Marshall, J. L., et al. 2016, AJ, 151, 157

Liu, C., Hu, J., Newberg, H., & Zhao, Y. 2008, A&A, 477, 139

Longeard, N., Martin, N., Starkenburg, E., et al. 2017, MNRAS, 480, 2609

Luque, E., Pieres, A., Santiago, B., et al. 2017, MNRAS, 468, 97

Luque, E., Queiroz, A., Santiago, B., et al. 2016, MNRAS, 458, 603

Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13

Malyshev, D., Neronov, A., & Eckert, D. 2014, PhRvD, 90, 103506

Martin, N. F., de Jong, J. T. A., & Rix, H.-W. 2008, ApJ, 684, 1075

Martin, N. F., Ibata, R. A., McConnachie, A. W., et al. 2013, ApJ, 776, 80

Martin, N. F., Nidever, D. L., Besla, G., et al. 2015, ApJL, 804, L5

Martínez-Vázquez, C. E., Monelli, M., Bono, G., et al. 2015, MNRAS, 454, 1509

Mashchenko, S., Wadsley, J., & Couchman, H. M. P. 2008, Sci, 319, 174

Mau, S., Cerny, W., Pace, A. B., et al. 2020, ApJ, 890, 136

Mau, S., Drlica-Wagner, A., Bechtol, K., et al. 2019, ApJ, 875, 154

McConnachie, A. W. 2012, AJ, 144, 4

Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19

Morganson, E., Gruendl, R. A., Menanteau, F., et al. 2018, PASP, 130, 074501

Moskowitz, A. G., & Walker, M. G. 2019, arXiv:1910.10134

Müller, O., Jerjen, H., & Binggeli, B. 2015, A&A, 583, A79

Muñoz, R. R., Côté, P., Santana, F. A., et al. 2018, ApJ, 860, 66

Munshi, F., Brooks, A. M., Christensen, C., et al. 2019, ApJ, 874, 40

Mutlu-Pakdil, B., Sand, D. J., Carlin, J. L., et al. 2018, ApJ, 863, 25

Nadler, E. O., Gluscevic, V., Boddy, K. K., & Wechsler, R. H. 2019a, ApJL, 878, L32

Nadler, E. O., Mao, Y.-Y., Green, G. M., & Wechsler, R. H. 2019b, ApJ, 873, 34

Nadler, E. O., Mao, Y.-Y., Wechsler, R. H., Garrison-Kimmel, S., & Wetzel, A. 2018, ApJ, 859, 129

Nadler, E. O., Wechsler, R., Bechtol, K., et al. 2019c, ApJ, 893, 48

Newton, O., Cautun, M., Jenkins, A., Frenk, C. S., & Helly, J. C. 2018, MNRAS, 479, 2853

Nilson, P. 1973, Uppsala General Catalogue of Galaxies (Uppsala: Roy. Soc. Sci. Uppsala)

Orear, J. 1958, Notes on Statistics for Physicists UCRL-8017

Pawlowski, M. S., Pflamm-Altenburg, J., & Kroupa, P. 2012, MNRAS, 423, 1109

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825

Pieres, A., Santiago, B. X., Drlica-Wagner, A., et al. 2017, MNRAS, 468, 1349

Pietrzyński, G., Górski, M., Gieren, W., et al. 2009, AJ, 138, 459

Pietrzyński, G., Graczyk, D., Gieren, W., et al. 2013, Natur, 495, 76

Plummer, H. C. 1911, MNRAS, 71, 460

Rykoff, E. S., Koester, B. P., Rozo, E., et al. 2012, ApJ, 746, 178

Rykoff, E. S., Rozo, E., & Keisler, R. 2015, arXiv:1509.00870

Sakamoto, T., & Hasegawa, T. 2006, ApJL, 653, L29

Samuel, J., Wetzel, A., Tollerud, E., et al. 2020, MNRAS, 491, 1471

Schlaflay, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlaflay, E. F., Finkbeiner, D. P., Schlegel, D. J., et al. 2010, ApJ, 725, 1175

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Sevilla-Noarbe, I., Hoyle, B., Marchà, M. J., et al. 2018, MNRAS, 481, 5451

Shapiro, P. R., Iliev, I. T., & Raga, A. C. 2004, MNRAS, 348, 753

Shapley, H. 1938a, BHarO, 908, 1

Shapley, H. 1938b, Natur, 142, 715

Sheldon, E. S. 2014, MNRAS, 444, L25

Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114

Simon, J. D. 2019, ARA&A, 57, 375

Smercina, A., Bell, E. F., Price, P. A., et al. 2018, ApJ, 863, 152

Somerville, R. S. 2002, ApJL, 572, L23

Spekkens, K., Mason, B. S., Aguirre, J. E., & Nhan, B. 2013, ApJ, 773, 61

Swanson, M. E. C., Tegmark, M., Hamilton, A. J. S., & Hill, J. C. 2008, MNRAS, 387, 1391

Tollerud, E. J., Bullock, J. S., Strigari, L. E., & Willman, B. 2008, ApJ, 688, 277

Tollerud, E. J., & Peek, J. E. G. 2018, ApJ, 857, 45

Tonry, J. L., Burke, B. E., Isani, S., Onaka, P. M., & Cooper, M. J. 2008, Proc. SPIE, 7021, 702105

Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2018, MNRAS, 475, 5085

Torrealba, G., Belokurov, V., & Koposov, S. E. 2019a, MNRAS, 484, 2181

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2019b, MNRAS, 488, 2743

Torrealba, G., Koposov, S. E., Belokurov, V., et al. 2016a, MNRAS, 463, 712

Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016b, MNRAS, 459, 2370

Walsh, S. M., Jerjen, H., & Willman, B. 2007, ApJL, 662, L83

Walsh, S. M., Willman, B., & Jerjen, H. 2009, AJ, 137, 450

Wang, M. Y., de Boer, T., Pieres, A., et al. 2019, ApJ, 881, 118

Webbink, R. F. 1985, in IAU Symp. 113, Dynamics of Star Clusters, ed. J. Goodman & P. Hut (Dordrecht: Reidel), 541

Weisz, D. R., & Boylan-Kolchin, M. 2017, MNRAS, 469, L83

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014a, ApJ, 789, 147

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014b, ApJ, 789, 148

Weisz, D. R., Koposov, S. E., Dolphin, A. E., et al. 2016, ApJ, 822, 32

Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, 827, L23

Wheeler, C., Hopkins, P. F., Pace, A. B., et al. 2019, MNRAS, 490, 4447

Wheeler, C., Oñorbe, J., Bullock, J. S., et al. 2015, MNRAS, 453, 1305

Wilks, S. 1938, Ann.Math.Statist., 9, 60

Willman, B., Blanton, M. R., West, A. A., et al. 2005a, AJ, 129, 2692

Willman, B., Dalcanton, J. J., Martinez-Delgado, D., et al. 2005b, ApJL, 626, L85

Wilson, A. G. 1955, PASP, 67, 27

Yniguez, B., Garrison-Kimmel, S., Boylan-Kolchin, M., & Bullock, J. S. 2014, MNRAS, 439, 73

Zucker, D. B., Belokurov, V., Evans, N. W., et al. 2006a, ApJL, 650, L41

Zucker, D. B., Belokurov, V., Evans, N. W., et al. 2006b, ApJL, 643, L103