# Identifying RR Lyrae Variable Stars in Six Years of the Dark Energy Survey

K.M. Stringer, A. Drlica-Wagner, L. Macri, C.E. Martínez-Vázquez, A.K. Vivas, P. Ferguson, A.B. Pace, A.R. Walker, E. Neilsen, K. Tavangar, et al.

► To cite this version:

K.M. Stringer, A. Drlica-Wagner, L. Macri, C.E. Martínez-Vázquez, A.K. Vivas, et al. Identifying RR Lyrae Variable Stars in Six Years of the Dark Energy Survey. The Astrophysical Journal, 2021, 911 (2), pp.109. 10.3847/1538-4357/abe873. hal-03070614

HAL Id: hal-03070614

https://hal.science/hal-03070614v1

Submitted on 1 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Identifying RR Lyrae Variable Stars in Six Years of the Dark Energy Survey

K. M. Stringer¹, A. Drlica-Wagner²,³,⁴, L. Macri¹, C. E. Martínez-Vázquez⁵, A. K. Vivas⁵, P. Ferguson¹, A. B. Pace⁶, A. R. Walker⁵, E. Neilsen², K. Tavangar³,⁴, W. Wester², T. M. C. Abbott⁵, M. Aguena⁷,⁸, S. Allam², D. Bacon⁹, K. Bechtol¹⁰, E. Bertin¹¹,¹², D. Brooks¹³, D. L. Burke¹⁴,¹⁵, A. Carnero Rosell¹⁶,¹⁷, M. Carrasco Kind¹⁸,¹⁹, J. Carretero²⁰, M. Costanzi²¹,²², M. Crocce²³,²⁴, L. N. da Costa⁸,²⁵, M. E. S. Pereira²⁶, J. De Vicente²⁷, S. Desai²⁸, H. T. Diehl², P. Doel¹³, I. Ferrero²⁹, J. García-Bellido³⁰, E. Gaztanaga²³,²⁴, D. W. Gerdes²⁶,³¹, D. Gruen¹⁴,¹⁵,³², R. A. Gruendl¹⁸,¹⁹, J. Gschwend⁸,²⁵, G. Gutierrez², S. R. Hinton³³, D. L. Hollowood³⁴, K. Honscheid³⁵,³⁶, B. Hoyle³⁷,³⁸,³⁹, D. J. James⁴⁰, K. Kuehn⁴¹,⁴², N. Kuropatkin², T. S. Li⁴³,⁴⁴, M. A. G. Maia⁸,²⁵, J. L. Marshall¹, F. Menanteau¹⁸,¹⁹, R. Miquel²⁰,⁴⁵, R. Morgan¹⁰, R. L. C. Ogando⁸,²⁵, A. Palmese²,³, F. Paz-Chinchón¹⁹,⁴⁶, A. A. Plazas⁴³, A. Roodman¹⁴,¹⁵, E. Sanchez²⁷, M. Schubnell²⁶, S. Serrano²³,²⁴, I. Sevilla-Noarbe²⁷, M. Smith⁴⁷, M. Soares-Santos²⁶, E. Suchyta⁴⁸, G. Tarle²⁶, D. Thomas⁹, C. To¹⁴,¹⁵,³², T. N. Varga³⁸,³⁹, R. D. Wilkinson⁴⁹, and Y. Zhang²

(DES Collaboration)

¹ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA; stringkm@gmail.com

² 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

⁵ Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile

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

⁷ Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

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

⁹ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

¹⁰ Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706-1390, 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

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

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

¹⁶ Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

¹⁷ Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

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

¹⁹ National Center for Supercomputing Applications, 1205 West Clark 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

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

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

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

²⁶ Department of Physics, University of Michigan, Ann Arbor, 48109, USA

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

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

²⁹ Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway

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

³¹ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

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

³³ School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

³⁴ Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

³⁵ 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

³⁷ Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany

³⁸ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany

³⁹ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany

⁴⁰ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁴¹ Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

⁴² Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA

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

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

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

⁴⁶ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁴⁷ 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

⁴⁹ Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

Received 2020 November 27; revised 2021 January 29; accepted 2021 January 31; published 2021 April 22
Abstract

We present a search for RR Lyrae stars using the full six-year data set from the Dark Energy Survey covering ~5000 deg² of the southern sky. Using a multistage multivariate classification and light-curve template-fitting scheme, we identify RR Lyrae candidates with a median of 35 observations per candidate. We detect 6971 RR Lyrae candidates out to ~335 kpc, and we estimate that our sample is >70% complete at ~150 kpc. We find excellent agreement with other wide-area RR Lyrae catalogs and RR Lyrae studies targeting the Magellanic Clouds and other Milky Way satellite galaxies. We fit the smooth stellar halo density profile using a broken-power-law model with fixed halo flattening (q = 0.7), and we find strong evidence for a break at R₀ = 32.1<sup>+1.1</sup><sub>-0.9</sub> kpc with an inner slope of n₁ = -2.54<sup>+0.09</sup><sub>-0.09</sub> and an outer slope of n₂ = -5.42<sup>+0.13</sup><sub>-0.14</sub>. We use our catalog to perform a search for Milky Way satellite galaxies with large sizes and low luminosities. Using a set of simulated satellite galaxies, we find that our RR Lyrae-based search is more sensitive than those using resolved stellar populations in the regime of large (r<sub>h</sub> ≳ 500 pc), low-surface-brightness dwarf galaxies. A blind search for large, diffuse satellites yields three candidate substructures. The first can be confidently associated with the dwarf galaxy Eridanus II. The second has a distance and proper motion similar to the ultrafaint dwarf galaxy Tucana II but is separated by ~5 deg. The third is close in projection to the globular cluster NGC 1851 but is ~10 kpc more distant and appears to differ in proper motion.

Unified Astronomy Thesaurus concepts: RR Lyrae variable stars (1410); RRab variable stars (1413); Milky Way stellar halo (1060); Milky Way Galaxy (1054)

Supporting material: machine-readable table

**1. Introduction**

Studies of the Milky Way stellar halo provide unique insights into the formation and evolution of our Galaxy (e.g., Johnston et al. 2002; Helmi 2008). Over the past several decades, wide-area digital sky surveys have shown that the Galactic halo hosts a large population of stellar substructures that can be classified as satellite galaxies, star clusters, and stellar streams (e.g., Belokurov et al. 2006; McConnachie 2012; Shipp et al. 2018; Simon et al. 2020). The abundance, distribution, and properties of halo substructures can be used to inform models for the assembly, chemical evolution, and star formation history of our Galaxy (e.g., White & Frenk 1991; Johnston et al. 2008; Tolstoy et al. 2009; Mo et al. 2010; Sharma et al. 2011a; Gallart et al. 2019). In addition, halo structure and substructure are valuable tools for estimating the matter density profile and total mass of the Milky Way (e.g., Deason et al. 2012; Kafle et al. 2012; Bonaca & Hogg 2018). A wide variety of luminous tracers have been used to map the structure and substructure of the Milky Way halo, including main-sequence turnoff stars (e.g., Belokurov et al. 2006; Shipp et al. 2018), blue horizontal branch stars (e.g., Deason et al. 2014), and red giant branch stars (e.g., Sharma et al. 2011b; Sheffield et al. 2014). However, among these stellar tracers, pulsating variable RR Lyrae stars (RRL) are especially useful, due to their distinct temporal signature and standardizable luminosity.

RRL are low-mass stars in the core helium-burning phase of evolution that radially pulsate when they fall within the instability strip (e.g., Walker 1989; Smith 1995; Bono et al. 2011; Marconi 2012). They are found in the horizontal branches of old stellar systems (>10 Gyr) and follow a well-understood period-luminosity-metallicity (PLZ) relation (e.g., Cáceres & Catelan 2008; Marconi et al. 2015). Their age, relatively high luminosity ($M_V = 0.59$ at $[Fe/H] = -1.5$; Cacciari & Clementini 2003), and well-understood PLZ relation make them excellent distance indicators for old, low-metallicity stellar populations in the outer halo of the Milky Way (e.g., Catelan et al. 2004; Vivas et al. 2004; Cáceres & Catelan 2008; Sesar et al. 2010; Stetson et al. 2014; Fiorentino et al. 2015). RRL are sufficiently luminous to be detected at

large distances and are sufficiently numerous to trace the halo substructures with good spatial resolution (e.g., Sesar et al. 2010, 2014; Baker & Willman 2015; Martínez-Vázquez et al. 2019; Torrealba et al. 2015; Vivas & Zinn 2006; Vivas et al. 2020a).

The most numerous subtype of RRL are RRab, which pulsate in the fundamental mode and have light-curve shapes resembling a sawtooth curve with short periods (0.4 ≲ *P* ≲ 0.9 days) and large amplitudes (0.5 ≲ *A*<sub>*g*</sub> ≲ 1.5 mag). In contrast, RRc pulsate in the first overtone and have smoother, more sinusoidal-shaped light curves with shorter periods (0.2 ≲ *P* ≲ 0.5 days) and smaller amplitudes (0.2 ≲ *A*<sub>*g*</sub> ≲ 0.5 mag). The detection and classification of RRL that contain additional pulsation modes require very well-sampled light curves over a long and continuous baseline. For example, RRd pulsate simultaneously in the fundamental and first overtone (e.g., Jerzykiewicz & Wenzel 1977), while RRL may be subject to the Blazhko effect (Blažko 1907; Buchler & Kolláth 2011), that is, a modulation of period and amplitude of unknown origin that can span several to hundreds of days.

Thanks to their distinct light curves and well-defined PLZ relation, RRL overdensities have been shown to be a good tracer of halo substructure (e.g., Vivas et al. 2001; Ivezic et al. 2004; Sesar et al. 2014; Baker & Willman 2015; Sanderson et al. 2017). Indeed, RRL have been detected in nearly all Milky Way satellite galaxies (e.g., see the recent compilation in Martínez-Vázquez et al. 2019) and are abundant in Milky Way stellar streams (e.g., Mateu et al. 2018; Price-Whelan et al. 2019; Ramos et al. 2020; Koposov et al. 2019). While Martínez-Vázquez et al. (2019) and Vivas et al. (2020a) showed that galaxies fainter than $M_V \sim -4.5$ are expected to contain fewer than three of these variables, even a few tightly clustered RRL in the outer halo could indicate the presence of an ultrafaint galaxy (Baker & Willman 2015). Such a technique was recently used to aid in the discovery of the ultradiffuse satellite Antlia 2 (Torrealba et al. 2019).

As the sensitivity of wide-field optical imaging surveys has increased, it has become possible to use RRL to probe the outer halo of our Milky Way at increasingly large distances. However, the temporal coverage of these surveys can be
sparse and nonuniform, requiring the development of statistical algorithms to detect and measure RRL (Hernitschek et al. 2016; Sesar et al. 2017; Medina et al. 2018). In Stringer et al. (2019, hereafter S19), we showed that a substantial number of RRab can be detected even in extremely sparsely sampled multiband light curves from the first three years of the Dark Energy Survey (DES; DES Collaboration 2005; DES Collaboration et al. 2016). Here, we extend this work to use the full six-year DES data set (DES Y6) to assemble a catalog of RRL over 5,000 deg² in the southern Galactic cap with sensitivity out to a heliocentric distance of ~500 kpc.

On average, DES Y6 has approximately twice as many observations of each astronomical source as the three-year data (DES Y3) explored in S19. This larger data set allows us to perform better identification and characterization of RRL candidates. While we show that our catalog agrees well with other overlapping surveys, DES Y6 only provides a median of 35 observations (combining all filters) per RRL candidate. High-cadence follow-up observations will be able to confirm and better characterize candidates in our sample. The catalog resulting from our analysis of the DES Y6 data consists of the locations, periods, and estimated distances of 6971 RRL, with the most distant candidate residing at ~335 kpc. We clearly resolve RRL structures associated with classical Milky Way satellite galaxies (i.e., the Large Magellanic Cloud, Small Magellanic Cloud, Fornax, and Sculptor), we detect previously known RRL associated with Milky Way ultrafaint satellite galaxies (Tucana II, Phoenix II, and Grus I), and we report the first candidate RRL associated with the ultrafaint satellites Eridanus II, Cetus III, and Tucana IV. Based on the successful detection of RRL associated with known ultrafaint satellites, we use our catalog to perform a search for previously undiscovered satellite galaxies in the DES footprint. No high-confidence satellite galaxy candidates are discovered, and we interpret the sensitivity of our search in the context of a suite of satellite galaxy simulations from Drlica-Wagner et al. (2020).

This paper is organized as follows. In Section 2, we present the DES Y6 single-epoch catalog data, our criteria for selecting stellar objects, and our calibration of the photometric uncertainties of steady sources. In Section 3, we describe the color and variability criteria used to select a set of objects for further analysis. In Section 4, we describe the RRL light-curve template-fitting procedure, which yields our catalog of candidate RRL. In Section 5, we discuss our resulting catalog of candidate RRab. We estimate the total efficiency of our identification techniques (Section 5.1), and we associate our RRab catalog with known classical dwarf satellite galaxies, ultrafaint satellites, and globular clusters residing in the DES footprint (Sections 5.2–5.5). In Section 5.6 we use our catalog to estimate the halo density profile. In Section 5.7, we perform a search for low-surface-brightness substructures using our RRab catalog. We state our conclusions in Section 6. The catalog of DES Y6 RRab candidates is available online.<fnref n="50" />

## 2. Data Preparation

### 2.1. DES Y6 Quick Catalog

DES (DES Collaboration 2005; DES Collaboration et al. 2016) was a six-year optical/near-infrared imaging survey covering ~5000 deg² 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). Observations were completed in 2019 January. DES obtained ~10 × 90 s exposures in five broadband filters, grizY (Neilsen et al. 2019).<fnref n="51" /> DES observed with gri in dark time, iz in gray time, and zY in bright time, with each field of the footprint being observed two to three times per year (Diehl et al. 2016, 2019; Neilsen et al. 2019). The 5σ point-source depth of the DES exposures is estimated to be grizY ~ (24.3, 24.1, 23.5, 22.9, 21.5) (Morganson et al. 2018).

As in S19, the light curves for this work were assembled using the internal DES “Quick” release pipeline. The DES Y6 Quick Release catalog (hereafter Y6Q) was constructed using survey exposures processed with the “Final Cut” pipeline from the DES Data Management system (DESDM; Morganson et al. 2018). This pipeline applies instrumental calibrations and detrending corrections to the images, then creates photometric source catalogs for each exposure using SourceExtractor (Bertin & Arnouts 1996). The full details of the DESDM image-reduction and catalog-creation pipelines are summarized in Morganson et al. (2018); we note that the Y6 processing used a lower source detection threshold (Y6 single-epoch catalogs have a detection threshold of ~3σ compared to a Y3 threshold of ~5σ), and it has an improved astrometric calibration based on Gaia DR2. All coordinates used in this paper are in the (J2000) equinox. The photometric calibration is performed with the Forward Global Calibration Module (FGCM; Burke et al. 2018). The relative photometric calibration accuracy in griz is estimated to be better than 3 mmag across the footprint (Sevilla-Noarbe et al. 2020), while the absolute photometric calibration in these bands is estimated to be ~3 mmag based on a comparison with the Hubble Space Telescope standard star C26202 (DES Collaboration et al. 2018). Complete details of the Y6 data processing, calibration, and validation will be released in forthcoming publications by the DES Collaboration.

After the images were reduced through the Final Cut pipeline, several quality cuts were applied to select exposures for the single-epoch catalog. Any images with insufficient depth, poor seeing, poor sky subtraction, astrometric errors, or that contained artifacts such as ghosts, bleed trails, and airplane streaks were excluded. Additionally, only exposures with FGCM zero-point solutions (Burke et al. 2018) were included in the catalog. These selections were applied to images from years one through six of survey operations (Y6), yielding a total of 78,364 exposures.<fnref n="52" />

We assembled a Y6Q unique catalog of astronomical objects by matching sources in a given image to the nearest neighboring detections (within 1" in radius) in all other exposures using CKDTree as implemented in scipy.spatial (Virtanen et al. 2020). When multiple detections from one exposure were located within 1", these were split into multiple objects in the Y6Q catalog. All objects with at least one detection in any band were included in the Y6Q catalog to ensure the inclusion of transient and moving objects. Additionally, the Y6Q catalog was cross-matched with objects detected in the Y6A1_COADD images to provide easy reference

<footnote n="50">⁵⁰ https://des.ncsa.illinois.edu/releases/other/y6-rrl</footnote>

<footnote n="51">⁵¹ DES took 45 s exposures in the Y band in the first three years.</footnote>

<footnote n="52">⁵² The DES Y6Q catalog does not include observations from the DES Science Verification period.</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: reduced median chi-squared in the *r* band, log₁₀(χ²*ν*,*r*), versus median *r*-band magnitude, *m*<sub>*r*</sub>, for objects classified as stars in a single DES HEALPix pixel (*l*, *b* ~ 302, ~ 52). Unlike DES Y3, the photometric uncertainties are found to be slightly underestimated for brighter objects in the DES Y6 data. Blue points show objects that were used to fit the quadratic curve (black line), while red points were objects excluded by a 3σ outlier clipping. Right: We correct this trend by rescaling the photometric errors using the quadratic fit so that the remaining trend in log₁₀(χ²*ν*,*r*) is flat. Although this trend is very minor, it is necessary to correct because the trends differ as a function of band and sky position.</figcaption></figure>

to quantities only available from the DES processing of the coadded images.

The resulting Y6Q catalog contained ~610 million objects distributed across the DES survey footprint. A large number of these objects possess only a single detection, which can occur due to spurious background fluctuations or transient objects. Overall, objects in the Y6Q catalog had a median of six observations spread over the grizY bands. If we require that objects be detected at least once in every band, then the catalog is reduced to ~109 million objects and the median total number of observations for all objects across all bands is 34. This nearly doubles the median total number of observations for RRL identified in the DES Y3 release (S19). Because of the lower signal-to-noise threshold for detections, this catalog is deeper than the DES Y3 one by ~0.75 mag in each band.

Unless explicitly stated otherwise, all magnitudes in this paper are point-spread function (PSF) magnitudes derived by SourceExtractor and corrected for interstellar extinction. Interstellar extinction is calculated following the prescription in DES Collaboration et al. (2018) using the Schlegel et al. (1998) dust maps with the normalization correction from Schlafly & Finkbeiner (2011) and the Fitzpatrick (1999) reddening law.

## 2.2. Rescaling Photometric Uncertainties

Similar to S19, we found magnitude-dependent residuals in the reduced chi-squared of the photometric measurements of objects classified as stars using the criterion described in the following section. The residuals found in DES Y6 were significantly smaller than those found in DES Y3, but they were still large enough to bias the identification of variable sources if left uncorrected. Thus, we followed the same procedure as in S19 to rescale the magnitude errors of each observation according to trends in the reduced chi-squared. Within each region of the survey (HEALPix pixel with nside=32), we define the reduced median chi-squared in band *b* as

$$ \tilde{\chi}_{\nu,b}^2 = \frac{1}{N_b - 1} \sum_{i=1}^{N_b} \frac{(m_{i,b} - \text{med}(m_b))^2}{\sigma_{i,b}^2}, \quad (1) $$

where $m_{i,b}$ is the observed PSF magnitude of an object in observation *i*, $\sigma_{i,b}$ is the reported magnitude uncertainty on that observation, and $N_b$ is the total number of observations of that

object in that band. We fit a quadratic function of the form

$$ \log_{10}(\tilde{\chi}_{\nu,b}^2) = c_{0,b} + c_{1,b}(\text{med}(m_b) - 20) + c_{2,b}(\text{med}(m_b) - 20)^2. \quad (2) $$

Figure 1 shows a noticeable trend in log₁₀(χ²*ν*,*r*), similar to those seen in S19; however, we find that the uncertainties are slightly underestimated for bright objects in DES Y6, as shown by the negative slope in Figure 1. Although this trend is far less pronounced than the trend in DES Y3, we perform this correction since the trends vary slightly over the wide-field footprint. We independently fit the coefficient for each HEALPix region in each band and corrected the uncertainties of the observed MAG_PSF quantities using the appropriate scale factor calculated from these relations. This process effectively rescaled the errors and flattened the trends in χ²*ν*,*b*.

# 3. Selection for Template Fitting

## 3.1. Stellar Source Selection

Since the faint end of the DES object sample is dominated by galaxies, we perform an initial star–galaxy separation to select stellar sources. We use the SPREAD_MODEL_I parameter from the exposure with the largest effective exposure time (Neilsen et al. 2016). The *i* band is preferred for star–galaxy separation because it typically has the best seeing of the *gri* bands observed during dark time (see Section 2.3 in DES Collaboration et al. 2018 and Figure 8 in Diehl et al. 2019). This procedure differs from S19, which considered all objects that passed this SPREAD_MODEL criterion in any of the griz bands to avoid omitting objects that were missing observations in a single band. While this enabled the catalog in S19 to be more complete, many extended sources leaked into the sample and had to be removed through visual inspection. In comparison, a much larger fraction of DES Y6 sources have at least one measurement in the *i* band. Any object for which |SPREAD_MODEL_I| < (0.003 + SPREADERR_MODEL_I) was considered (~310 million sources from the Y6Q catalog pass this cut). Additionally, only objects that were associated with an object detected in the Y6A1 coadded images were considered (0.7 match radius).
## 3.2. External Catalogs and Simulated Data Used to Define the Sample

After rescaling the photometric uncertainties and applying the star–galaxy separation, we remove any objects with fewer than 10 total observations. Such a small number of observations skews their variability statistics and makes template fitting challenging. We then further reduce the size of our catalog by selecting objects that match the colors and temporal variability that are characteristic of RRL.

To determine optimal color and variability cuts to select RRL, we cross-match objects in Y6Q with RRL identified by external surveys, variable objects that are not RRL, and stars used for the photometric calibration that are largely nonvariable. Our sample of external RRL includes objects from Sloan Digital Sky Survey (SDSS) Stripe 82 (hereafter S10; Sesar et al. 2010), the Catalina Sky Surveys DR2 (Drake et al. 2013a, 2013b, 2014; Torrealba et al. 2015; Drake et al. 2017), Pan-STARRS PS1 (Sesar et al. 2017), variables from the Sculptor dSph (Martínez-Vázquez et al. 2016a), RRL from the Fornax dSph (Bersier & Wood 2002), and RRL from Gaia DR2 with measured periods (Holl et al. 2018; Clementini et al. 2019; Rimoldini et al. 2019). As a likely contaminant class, we also select ~16,000 quasars (QSOs) from the KiDS DR3 survey (Nakoneczny et al. 2019; de Jong et al. 2017) and the SDSS-POSS southern sample (MacLeod et al. 2012). Finally, to compare against other (likely nonvariable) contaminants, we match to an internal DES catalog of ~17 million stars that were used in the FGCM zero-point calibration (Burke et al. 2018). We remove calibration stars located less than 10 arcmin from the centers of the Fornax and Sculptor dwarf galaxies since the DES photometry suffers from crowding in these regions. We randomly downsample the QSO and calibration star catalogs to match the size of our external RRab sample. The resulting comparison samples contain 5055 RRab, 472 RRc, 46 RRd, 4 Blazkho RRL, 5055 QSOs, and 5055 calibration stars.

To further guide our selection criteria, specifically at faint magnitudes, we simulate a set of mock RRL light curves. This process follows the procedure described in S19, and we only provide a brief summary below. Our simulated light curves are based on well-sampled light curves of RRL from S10. First, we construct smoothed light-curve shapes using the best-fitting templates and observational parameters for each of the 483 RRL (379 RRab and 104 RRc) identified by S10 and convert their magnitudes into the DES filter system.<fnref n="53" /> We then subtract the S10 estimated distance moduli from each light curve to transform to absolute magnitude. For a set of magnitude bins in the range $15.5 \le g \le 24.5$ with a bin width of 0.5, we shift the light curve for each of the 483 light-curve shapes to an average g magnitude randomly drawn from a uniform distribution within that magnitude bin. As each light curve is shifted to a random distance, any simulated measurement fainter than the Y6 single-epoch limiting magnitude in that band is removed.<fnref n="54" /> We then assign photometric uncertainties to each observation using the rescaled values from Section 2.2 and use them to introduce scatter into the observations. The light curves are then downsampled to the DES observing cadence as determined from a random set of bright stars in DES. The total

simulated sample includes 8211 light curves, of which 6443 were RRab and 1768 were RRc. We do not specifically search for RRc in this work, but include their simulated curves to assess the RRc contamination in our final sample.

## 3.3. Color Selections

RRL inhabit a well-defined region of color–color space (e.g., Ivezić et al. 2005). We therefore apply a color selection to further reduce the number of objects that are passed to the light-curve template-fitting stage. This selection specifically excludes variable stars that are too red to be RRL (e.g., low-mass main-sequence stars). Since the colors of RRL change over the course of their pulsation cycle (e.g., Guldenschuh et al. 2005; Vivas et al. 2017), calculations based on observations obtained at multiple phases will degrade the separation power of this method. Thus, we take advantage of the DES observing cadence to measure “instantaneous colors” based on observations taken within one hour of each other. To reduce the time spent slewing between fields, DES often took sequences of two to three consecutive exposures of the same field in different filters, with the filters chosen according to the seeing, lunar phase, and number of previous observations (see Figure 3 in Diehl et al. 2016). To select these sequences, we group together any observations taken within one hour of each other and, depending on the filters used, calculate $g-r$, $r-i$, or $i-z$ colors. Due to the observing cadence of DES, the median separation time for each color combination is ~120 s. If there are multiple instantaneous colors for an object, we store only the maximum and minimum values.

We develop a multistage color selection to remove objects based on their instantaneous colors (when available), while retaining any objects that did not have a particular instantaneous color available. Using our sample of previously identified RRab, we defined selections as the 99% percentile value of the RRL population (Figure 2). The first two selections use the extinction-corrected minimum instantaneous colors, $(r-i)_{min,0}$ and $(g-r)_{min,0}$, to select objects lying near the blue end of the stellar locus. The first cut selects objects with $(r-i)_{min,0} - \sigma(r-i)_{min} \le 0.068$. Any objects that do not have an instantaneous measurement of $r-i$ but have one in $g-r$ are passed into the next step, which retains all objects with $(g-r)_{min,0} - \sigma(g-r)_{min} \le 0.268$. Any objects that do not have instantaneous measurements of $r-i$ or $g-r$ are retained if they satisfy $(\text{WAVG\_MAG\_PSF\_G} - \text{WAVG\_MAG\_PSF\_R})_0 \le 0.488$. As can be seen from Figure 2, this is the least restrictive of the three cuts. In total, ~51 million stellar sources pass our color cuts, including 98.65% of the sample of previously known RRab.

## 3.4. Variability Selection

We select temporally variable objects by calculating several variability statistics. Many of these quantities, which are summarized in Table 4, are based on the analysis of Sokolovsky et al. (2017) and are described in detail in Appendix A. We calculate these variability statistics from the MAG_PSF measurements and their rescaled uncertainties as described in Section 2.2.

Since numerous color, magnitude, and variability measurements are calculated (many of which are correlated), we use the random forest algorithm to “learn” the optimal boundaries in feature space to separate RRab from non-RRab. Random

<footnote n="53">⁵³ The SDSS-DES filter transformation equations can be found in Appendix F of S19.</footnote>

<footnote n="54">⁵⁴ The single-epoch limiting magnitude was estimated from Table 1 of DES Collaboration et al. (2018) with an adjustment to account for the lower object-detection threshold in Y6.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Multistage color cuts applied to DES stars to remove non-RRL. Each panel shows the distribution of colors for RRab in red, QSOs in blue, and stars used for photometric calibration in gray. The black dashed line indicates the 99th percentile of the RRab distribution and was set as the threshold for inclusion in our sample. Top: distribution of instantaneous $(r - i)_{min,0}$. Middle: distribution of instantaneous $(g - r)_{min,0}$. Bottom: distribution of the extinction-corrected $(g - r)_0$ using weighted-average magnitudes.</figcaption></figure>

forests use a collection of decision trees to predict an object's type. Each decision tree is trained on a random subset of the training population by repeatedly subdividing the sample based on feature values until a user-defined maximum depth or other specified stopping condition is reached. The specific features that are used in each split are chosen randomly and can influence the characteristics of the population that a tree learns to detect. A random forest classifier takes advantage of the learning differences between individual trees by averaging the predictions from a large number of trees to produce an aggregate score (Amit & Geman 1997; Breiman 2001). Random forests are a good choice for this task because they are largely insensitive to uninformative features, so the inclusion of a feature that does not separate the object types will not harm the overall results. For our classifiers, we use

the `RandomForest` implementation in `scikit-learn` (Pedregosa et al. 2012). To reduce our sample of light curves down to a manageable size for template fitting, we divide our pretemplate selection criteria into two phases: (1) remove nonvariable sources, and (2) remove common variable contaminants (i.e., QSOs).

For the first classifier, our training set consists of equal numbers of simulated RRab and calibration stars that are largely nonvariable. We choose to use simulated RRab instead of RRab cross-matched from external catalogs because the simulated RRab cover the entire magnitude range of DES Y6 and thus do a better job of including the decreasing sensitivity to variability for more distant objects with larger photometric uncertainties. Hyperparameters for the random forest classifiers are determined using the `GridSearchCV` function of `scikit-learn`. This first classifier has 35 trees with a depth of eight splits and eight features allowed at each split.

This stage of the selection is intended to remove as many nonvariable sources as possible, so we choose the cutoff classifier score using the $F_β$ score (Baeza-Yates & Ribeiro-Neto 1999),

$$F_{\beta} = \frac{(1 + \beta)^2 \cdot \text{TP}}{(1 + \beta^2) \cdot \text{TP} + \beta^2 \cdot \text{FN} + \text{FP}} \quad (3)$$

where TP (true positives) and FN (false negatives) reflect how many true RRab have classifier scores above or below the cutoff threshold, respectively. Similarly, the FP (false positives) and TN (true negatives) show the number of non-RRab with classifier scores above or below the cutoff threshold. We choose this particular score over other popular metrics like the “informedness” or “$F_1$-score” because we wish to prioritize purity at this stage in the analysis (Powers 2008). The $F_β$ score with $β = 0.5$ allows us to weigh the precision twice as heavily as the recall. The classifier value that maximizes $F_β$ is 0.755. This value yields a sample with an RRab precision of 99.71% and a recall (completeness) of 97.96% when applied to the training set. Approximately 20% of the 51 million input objects pass this classification.

As a second step, we train and apply a classifier optimized to remove variable objects that do not show the strong variability pattern of RRL. In particular, due to the sparse temporal sampling of DES, it can be difficult to distinguish the variability of QSOs from that of RRL (e.g., S19). In addition, QSOs can be unresolved, have blue colors similar to RRL, and are abundant at the faint magnitudes reached by DES (e.g., Tie et al. 2017). Thus, for our training set, we use an equal number of simulated RRab and real QSOs cross-matched from the KiDS DR3 survey (Nakoneczny et al. 2019; de Jong et al. 2017) and the SDSS-POSS southern sample (MacLeod et al. 2012). For this classifier, we use a random forest with 18 deeper trees with nine features allowed for consideration at a total of 16 splits. To prioritize RRL completeness, we use the $F_β$ score with $β = 1.5$ to choose a cutoff score of 0.34. This cutoff score returns a purity of 86.28% and completeness of 96.04% for the training set.

After applying both of these classifiers, ~1.1 million objects remain for light-curve template fitting. The performance curves for both of these initial variability classifiers and their top features are included in Appendix B.
# 4. Template Fitting and Classification

## 4.1. Template Fitting

Following S19, we fit a multiband RRab light-curve template to the extinction-corrected light curves of each object passing our aforementioned selection criteria. Our template is empirically derived from the well-sampled RRab light curves of S10. Our template-fitting procedure is particularly effective when applied to sparsely sampled multiband light curves because it solves for only four independent parameters: period (P), phase (φ), g-band amplitude (A<sub>g</sub>), and distance modulus (μ). In S19, we showed that these parameters can be effectively constrained with a small number of observations, and we refer the reader to that paper for a thorough discussion of this method.

We apply the same template-fitting procedure as S19 with only minor modifications. In particular, we expanded the period range from 0.44–0.89 days in S19 to 0.2–1 days. As in S19, the periods, amplitudes, phases, distance moduli, and residual sum of squares (RSS) per degree of freedom of the fits are kept for the three best-fitting templates. Fitting one light curve requires ~2–6 minutes, depending on the CPU load at the time of processing.

Although both RRab and RRc can be used to estimate distances (Cáceres & Catelan 2008; Marconi et al. 2015, and references therein), RRc are more easily confused with other types of variable objects due to their lower amplitude of variation and more sinusoidal light curves. This, combined with their lower rates of occurrence in halo RRL populations (Martínez-Vázquez et al. 2017), makes them less attractive targets for our analysis of Galactic structures in Section 5. Thus, our template fit and subsequent classification are optimized to identify and fit RRab light curves and properties. Some RRc do pass all of the steps of our identification and fitting procedure and are misidentified as RRab. We explore the recovery rate for RRc with simulated light curves in Section 5.1.

## 4.2. Classification

Even though our color and variability selections decrease by >100× the number of objects that require template fitting, there are still too many for individual visual inspection. Thus, we train another random forest classifier to select potential RRL. We input the RSS, amplitudes, and von Mises-Fisher concentration parameter κ (see Section 4.3 in S19 for more details) from the top three best template fits for each candidate. As in S19, we also calculate the distances of each set of periods and amplitudes to the scaled Oosterhoff relations (Oosterhoff 1939) parameterized by Fabrizio et al. (2019) and scaled to the g band by Vivas et al. (2020b, see their Figure 7).

For this training set, we used the previously known cross-matched RRab described in Section 2, instead of the simulated light curves. Although the former do not span the full magnitude range of DES Y6, they provide a more realistic representation of the performance of template fits to RRab beyond those from S10. The template-fitting procedure performs artificially too well on simulated light curves because the templates and simulations are based on the same data (S19). For the non-RRab set, we include objects likely to contaminate our sample, such as the cross-matched QSOs from Section 3.2 as well as objects that were visually rejected as extended sources and artifacts in a previous iteration of the catalog.

The cutoff score is selected using Matthew’s correlation coefficient (MCC), which balances true and false positives and negatives for binary classification problems, even in the case of imbalanced class sizes (Matthews 1975). The MCC is defined as

$$MCC = \frac{TP \times TN - FP \times FN}{\sqrt{(TP + FP)(TP + FN)(TN + FP)(TN + FN)}}, \quad (4)$$

where TP is the true positives, TN is the true negatives, FP is the false positives, and FN is the false negatives. We choose a classifier score cutoff of 0.605 that maximizes the MCC for the training sample. For the training set, this cutoff score recovers 96.07% of the input RRab with 98.92% purity. When we apply this classifier to the full results of the template fitting, 10,812 objects pass the cutoff as potential candidates. The performance curves and top features are shown in Appendix B. We discuss the recovery fraction for a broader range of RRL distances in Section 5.1.

Since the purity of the DES star–galaxy separation decreases dramatically at fainter magnitudes (e.g., Shipp et al. 2018), we visually inspect every candidate that had not been previously identified as an RRab by another survey (see Figure 3). During this visual inspection, we remove any obvious extended sources, bright oversaturated objects, and extremely poor-fitting light curves. After this visual inspection, our final sample consists of 6971 objects.

# 5. Results

Our selection process results in a catalog of 6971 objects that are consistent with the colors, variability, and light curves of RRab. We hereafter refer to objects in this sample as DES Y6 “RRab candidates.” We show a selection of light curves covering the full magnitude range of our sample in Figure 3. DES Y6 provides a median of 35 observations per RRab candidate. The equatorial positions and heliocentric distances of DES Y6 RRab candidates are shown in the top and bottom panels of Figure 4, respectively. The RRab candidates that were previously identified by other studies, including the DES Y3 RRab search of S19, are shown in Figure 5. Because the boundaries of the DES Y3 and Y6 catalogs differ (most significantly around the LMC outskirts), this is not a one-to-one comparison.

Although image cutouts and light-curve fits are visually inspected for each of these candidates, they still require additional observations to confirm that they are RRL. This is especially true of objects with sparsely sampled light curves or located beyond 150 kpc. The properties of visually accepted and previously identified RRab candidates are summarized in Table 1 (full version available online).

## 5.1. Completeness for Faint Objects

Due to the limited depth and area of external RRL catalogs, we are required to use simulated RRL to assess the completeness of our selection for faint RRL (D ≳ 100 kpc). To assess the recovery of our template-fitting and classification process and to assess how frequently RRc are misclassified as RRab, we follow the procedure described in Section 3.2 to simulate an independent set of light curves for ~3100 RRab and ~900 RRc. We subject these simulated light curves to the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Sample of RRab candidates selected to represent the magnitude range of our catalog. Light curves are arranged vertically by increasing magnitude with best-fit heliocentric distances ranging from 7 kpc (top left) to 280 kpc (lower right). Observed magnitudes and RRab template light curves are colored by band (g: green, r: red, i: orange, z: blue, Y: purple). If not visible, photometric errors are smaller than the plotting symbols.</figcaption></figure>

same minimum number of observations, color cuts, variability classifier, template fitting, and final classifier that are applied to the real data. The recovery rate and precision of the parameter estimates are shown in Figure 6.

We find that our ability to identify RRab and recover periods within 1% of the input value improves dramatically as the number of observations increases and degrades dramatically at distances ≳200 kpc. The uncertainties in period and amplitude for both RRab and RRc using the template-fitting technique improve when there are more observations available and decline for objects at larger distances since their photometric uncertainties are larger. The period recovery for RRab approaches 100% when the light curve has ≥30 observations. The fraction of RRc misidentified as RRab by the classification pipeline is overall very low (≤10% for light curves with ≥25 observations and <20% at all distances). Given that RRc make up ≲30% of RRL (Soszyński et al. 2019; Martínez-Vázquez et al. 2017), the low efficiency suggests that RRc are unlikely to contribute substantially to our final catalog.

Our recovery rate for distant RRab (≥50% for $\mu \le 22$) opens an exciting new discovery space for stellar structures beyond 100 kpc. In the following sections, we compare the DES Y6 RRab candidate sample to external catalogs in classical and ultrafaint dwarf galaxies. In addition, we use this catalog to fit the Milky Way stellar halo profile and search for previously undiscovered Milky Way satellites.

## 5.2. The Magellanic Clouds

The variable star content of the Magellanic Clouds (MCs) has been extensively studied, with the OGLE-IV catalog of Soszyński et al. (2019) containing 47,828 RRL in both MCs. Although the DES footprint only covers the outskirts of the

MCs, a comparison with the OGLE catalog still yields a relatively large number of objects in common (792) and allows us to assess the recovery of light-curve properties of the DES Y6 RRab candidates. The OGLE catalog is ideal for this comparison because its light curves have <100 epochs (and consequently, period estimates and classifications are very robust) and it contains many types of variables (allowing us to investigate the contamination rate of our catalog).

The variables in common between the two catalogs include 53 objects classified as RRc or RRd by OGLE, implying a ~7% contamination by these types in our sample. This is consistent with estimates from our simulations (Figure 6) under the assumption that stellar populations are similar in the Milky Way and LMC halos, and with the “cloud” seen in the period-amplitude diagram (left panel of Figure 7) at short periods and low amplitudes.

A total of 720 (91%) of the matches have an excellent agreement in period (right panel of Figure 7). The median difference in period for objects on the 1:1 line in Figure 7 is only 0.001%. Most period mismatches are due to 1-day aliasing.

We also search for matches between our catalog and other types of variables in OGLE, finding one object that they identified as an anomalous Cepheid (AC). The light curves of ACs and RRL are similar, and it is difficult to distinguish them in the field. ACs are rare, so the contamination of ACs in the DES Y6 RRab catalog is expected to be low. We found no matches with the OGLE catalog of 40,204 eclipsing variables in the MCs (Soszyński et al. 2019).

We search for new RRab in the MCs using the DES Y6 coverage beyond the OGLE footprint. We start by selecting a range of distance moduli for each MC based on stars in
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: sky location of DES Y6 RRab candidates colored by distance modulus. This figure uses a McBryde-Thomas flat polar quartic projection with the DES footprint shown in black. Bottom: heliocentric distance distribution of the DES Y6 RRab candidates as a function of R.A. Overdensities corresponding to the LMC ($\alpha \sim 80^{\circ}$), Fornax ($\alpha \sim 39^{\circ}$), and Sculptor ($\alpha \sim 15^{\circ}$) can be seen in both panels. Stars associated with Eridanus II can be seen in the bottom panel at $\alpha \sim 56^{\circ}$ and $D > 300$ kpc.</figcaption></figure>

common with OGLE. For the LMC, the distance modulus distribution showed a clear peak at $17.9 < \mu < 18.8$, containing 608 RRL. The SMC overlap with OGLE is significantly smaller, but there is still a clear peak in the distance modulus distribution at $18.1 < \mu < 19.2$, with 27 RRL matching OGLE in that distance range. Matched OGLE-DES stars are shown as black diamonds in Figure 8. We then select DES Y6 RRab candidates that were not in OGLE but have distance moduli

consistent with the selection defined above and have an angular separation of less than 25° and 15° from the centers of the LMC and SMC, respectively. These stars are potential new members of the MCs. Most of them are located outside the OGLE footprint, although there are also some stars that may have been missed by OGLE.

We use Gaia DR2 data (Gaia Collaboration et al. 2018) to check whether our candidate MC members exhibit proper
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Overlap with external RRab catalogs. The distribution of this RRab catalog as a function of the average extinction-corrected *r* magnitudes, ⟨*r*⟩, are plotted in red. The distributions of previously identified RRab from external catalogs that were recovered in our catalog are plotted in black. Bottom right: We note that the combination of the increased number of observations, decreased signal-to-noise threshold, and alternative selection cuts in Y6 increased the depth of our catalog by ~1 magnitude, but also extended the detection range to brighter magnitudes as well.</figcaption></figure>

<table><caption>**Table 1**
DES Y6 RRab Candidates</caption><thead><tr><th>DES Y6 ID</th><th>α<br>(deg, J2000)</th><th>δ</th><th>⟨g⟩</th><th>⟨r⟩</th><th>⟨i⟩<br>(mag)</th><th>⟨z⟩</th><th>⟨Y⟩</th><th>P<br>(days)</th><th>A<sub>g</sub><br>(mag)</th><th>μ</th></tr></thead><tbody><tr><td>871765223</td><td>300.6754</td><td>-53.9670</td><td>16.647</td><td>16.495</td><td>16.489</td><td>16.512</td><td>16.561</td><td>0.5956</td><td>0.881</td><td>15.98</td></tr><tr><td>871820547</td><td>300.4591</td><td>-50.7457</td><td>15.352</td><td>15.201</td><td>15.221</td><td>15.225</td><td>15.039</td><td>0.5441</td><td>1.107</td><td>14.55</td></tr><tr><td>872060389</td><td>301.1574</td><td>-54.3232</td><td>19.217</td><td>18.935</td><td>18.841</td><td>18.647</td><td>18.647</td><td>0.6158</td><td>1.181</td><td>18.14</td></tr><tr><td>872212127</td><td>301.4859</td><td>-57.0505</td><td>17.941</td><td>17.671</td><td>17.770</td><td>17.715</td><td>17.810</td><td>0.5638</td><td>1.059</td><td>17.13</td></tr><tr><td>872268895</td><td>301.5402</td><td>-52.5811</td><td>14.797</td><td>14.721</td><td>...</td><td>14.637</td><td>14.741</td><td>0.5807</td><td>0.670</td><td>14.17</td></tr><tr><td>872444946</td><td>300.5736</td><td>-56.2696</td><td>16.417</td><td>16.272</td><td>16.270</td><td>16.223</td><td>16.241</td><td>0.6750</td><td>0.974</td><td>15.68</td></tr></tbody></table>

**Note.** DES Y6 ID: DES Y6A1 COADD_OBJECT_ID number. α: R.A. δ: decl. ⟨grizY⟩: Mean extinction-corrected magnitude. P: Best-fit period. A<sub>g</sub>: Best-fit amplitude in DES g. μ: Best-fit distance modulus. The full version of this catalog, including feature values and cross-matching information, is available in the online data products at this URL.

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

motions (PMs) consistent with those of the MCs. In Figure 9 we show the PMs of the stars matched to OGLE in each of the MCs and compare them with the new member candidates. For clarity, we show two plots for the LMC: one containing new candidate members within 20° of its center (top panel), and another for those lying between 20° and 25° (center panel). We find that new LMC member candidates in the former group have a distribution of PMs that is consistent with OGLE stars, while the latter do not. We find 202 LMC candidate RRL out to 20°, excluding 16 members of the LMC globular cluster Reticulum (which we discuss separately in Section 5.5). Beyond this limit and out to 25°, the density of LMC candidate RRL drops appreciably with only 25 possible members.

The old stellar population of the LMC extends to large angular separations, as shown in previous studies. Saha et al. (2010) traced the LMC population with main-sequence turnoff stars out to 16° from its center. Using a similar technique, Nidever et al. (2019) further extended the detection to 21°. In addition, Belokurov & Koposov (2016) used publicly available DES images to find a lumpy distribution of blue horizontal

branch (BHB) stars extending out to 30° and possibly up to 50°. We note in passing that some of our candidates match the location and distance of their “S1” substructure. Our strong detection of RRL in the LMC periphery confirms these previous results and strengthens the case for an extended, old (≥10 Gyr) stellar population in the LMC.

The southern decl. limit of DES is δ = −65°3. Consequently, only the very external parts of the SMC (which is centered at δ ~ −72°8) are within the footprint of the survey. Nonetheless, the OGLE catalog contains SMC RRL as far north as δ ~ −62°7, 27 of which are in common with our sample (Figure 8). We find 18 other RRL within the distance range of the SMC out to an angular separation of 15 deg (a region beyond the OGLE coverage). There is good agreement in the PMs of these new member candidates, as seen in the bottom panel of Figure 9.

An interesting feature associated with the SMC in this part of the sky is the so-called Small Magellanic Cloud Northern Over-Density (SMCNOD; Pieres et al. 2017), whose area is indicated by a blue ellipse in the left panel of Figure 8. Pieres et al. (2017)
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Recovery of simulated RRL light curves. RRab results are plotted in red, while RRc are plotted in purple. All quantities in the left column are plotted as functions of the number of light-curve observations, $N_{\text{obs}}$, while the right column shows these same quantities plotted as a function of distance modulus, $\mu$. We do not explicitly search for RRc in this work, and these values can be used to assess the RRc contamination in our catalog. The quantities shown are as follows. First row: fraction of simulated curves identified as RRab by the full classification pipeline. Second row: fraction of periods correctly estimated within 1% of their input values. The period accuracy approaches 100% near $N_{\text{obs}} \ge 30$ and rapidly decreases beyond $\mu \sim 21$. Third row: standard deviation of period estimate precision. Bottom row: standard deviation of amplitude estimate precision. We see a degradation in precision of the amplitude measurement for light curves with $N_{\text{obs}} \lesssim 30$ and $\mu \gtrsim 21$, similar to that of the period estimates. Uncertainties on the recovery efficiency in the first two rows are assessed through the Bayesian technique of Paterno (2004), while the bottom two rows show uncertainties from bootstrap resampling (e.g., Efron 1982).</figcaption></figure>

estimated that the SMCNOD is primarily composed of an intermediate-age population of ~6 Gyr, which is too young to contain RRL. Indeed, the variable star content in the SMCNOD was examined by Prudil et al. (2018) with an earlier release of the OGLE catalog (Soszyński et al. 2017). Although they found eight OGLE RRab inside the SMCNOD area, they concluded their density was compatible with the expected SMC background at that distance. In contrast, we find almost twice as many (15 RRab candidates) in the same region.

### 5.3. Sculptor and Fornax

In addition to the MCs, two other prominent overdensities of RRab candidates are within the DES footprint: the classical dwarf spheroidal (dSph) galaxies Sculptor ($\alpha$, D = 15°02.86 kpc) and Fornax ($\alpha$, D = 39°96, 147 kpc). Both can be clearly seen in Figure 4.

Martínez-Vázquez et al. (2016a, hereafter MV16) presented the most complete and extensive study of the variable star population in Sculptor, reporting 536 RRL, of which 289 were
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: period-amplitude diagram for DES Y6 RRab candidates. Overplotted are the Oosterhoff relations parameterized by Fabrizio et al. (2019) and scaled to the g band by Vivas et al. (2020b). Right: comparison of the periods of 792 RRL in common with OGLE. The red line is a 1:1 relation (91% of objects), while the cyan lines correspond to the ±1-day aliases.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Aitoff projection of the sky in equatorial coordinates showing the region around the LMC (left) and SMC (right). In both panels, the gray points show the density of OGLE RRL, the black diamonds show the DES RRab candidates matched with OGLE, while the red circles are DES RRab candidates that are not matched to OGLE but have distance moduli in agreement with the distance of each galaxy. Left: Teal points show possible LMC member stars that have distance moduli matching the LMC ($17.9 < \mu < 18.8$) but are separated from the LMC by $20^\circ-25^\circ$. The green square shows the location of the Reticulum globular cluster, which has 16 DES Y6 RRab candidates associated with it. Right: RRab candidates potentially associated with the SMC ($18.1 < \mu < 19.2$); the location of the SMCNOD (Pieres et al. 2017) is indicated with a blue ellipse.</figcaption></figure>

RRab. We search for potential Sculptor members by selecting all RRL in our catalog out to $1.5 \times$ its tidal radius of 69.1 (Munoz et al. 2018) and with distance moduli within $19.67 \pm 0.75$. We find 116 RRL within these limits, all but four matching the catalog of MV16. The lower-than-average completeness in this region ($39\%$ at $\mu \sim 19.5$) is due to crowding effects near the center of the galaxy.

MV16 classified six of our RRL candidates as RRc/RRd. This represents a 5% contamination, in agreement with the value derived from the MCs (Section 5.2). Removing these misclassified stars, we find that the periods for the rest of our sample agree very well with those from MV16, with a median difference of only 0.002%.

Three of the four new candidate members associated with this galaxy lie outside the footprint explored by MV16. One is located at 93', significantly beyond the tidal radius of Sculptor but with a compatible Gaia DR2 proper motion. The mean distance modulus of our Sculptor RRL is $\langle\mu\rangle = 19.49 \pm 0.09$ mag, while Martínez-Vázquez et al. (2015) found a value of $19.62 \pm 0.04$ mag. Although small, the difference can be attributed to the fact that our analysis assumes all RRL have the mean metallicity of the Milky Way halo, while the old stellar

population in Sculptor has a wide range of metallicities extending down to [Fe/H] $\sim -2.4$ (Martínez-Vázquez et al. 2016b).

While more distant than the MCs and Sculptor, Fornax also displays a prominent overdensity of DES Y6 RRL candidates. Again, we select candidate members as stars within $1.5 \times$ of its tidal radius of 77.5 (Wang et al. 2019) and distance moduli within $20.82 \pm 0.75$, yielding 1385 RRL. The mean distance modulus of these stars is $20.67 \pm 0.08$, which is again slightly smaller than the best-available value of $20.82 \pm 0.02$ (Karczmarek et al. 2017), likely due to the higher metallicity of our template.

Surprisingly, 51 of our Fornax RRL candidates are located beyond its tidal radius. These stars are uniformly distributed around the galaxy (Figure 10), reaching out to 114' from its center (note that our search radius only extended to 116'.2). In a recent study based on data from DES Y3, Wang et al. (2019) concluded that no significant extratidal disturbances are observed down to a surface brightness limit of $\sim 32.1$ mag arcsec$^{-2}$. Fornax is known to have several bursts of star formation and is dominated by a population of age $\sim 5$ Gyr (de Boer et al. 2012; Rusakov et al. 2021). While this dominant population is too young to have produced the observed RRL,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Gaia DR2 proper motions of DES Y6 RRab candidates around the LMC (top and middle panels) and in the SMC (bottom panel). In all panels, the black points are candidates that match the OGLE catalog. For the LMC, we show separately the candidates closer to the center of the galaxy in red (top panel) and candidates with larger separation in teal (middle panel). For the SMC, we show new candidate members within 15° in red.</figcaption></figure>

an old stellar population (>10 Gyr) is also present in Fornax. This older stellar population has been found to be more spatially extended than the younger populations (Wang et al. 2019). The detection of RRL outside the tidal radius suggests a very low surface brightness, old extratidal population.

Fornax is known to be rich in RRL, but a complete census is not readily available yet. Bersier & Wood (2002) found 525

<figure><img src="image_12.png" /><figcaption>**Figure 10.** RRL from DES in the Fornax dSph galaxy. The black ellipse shows the tidal radius and structural parameters as derived by Wang et al. (2019). The red crosses indicate the RRL that match Bersier & Wood (2002). Magenta triangles are possible anomalous Cepheids in Fornax that were confused as RRL in the DES catalog.</figcaption></figure>

RRL in the central part of Fornax, although the quality of their light curves was not good enough to do a proper classification into RRab and R Rc types. Based exclusively on their periods, they estimated that 396 of those stars may be RRab. We identified 275 of those stars in our catalog, shown in red in Figure 10, which suggests a 69% recovery rate, higher than in Sculptor. A more complete search for RRL was most recently made by Fiorentino et al. (2017). In this work, they found 990 RRab and 436 RRc in a region of 54' × 50' centered on the galaxy. Unfortunately, this catalog is not publicly available, preventing a direct comparison. However, our results suggest that the total population of RRL in Fornax is not yet known, especially at large angular separations.

We notice that within the Fornax search area, there is a group of seven stars with $19.5 < g < 20.7$, ~1 mag brighter than the Fornax RRL. It is unlikely that these are field halo stars, since RRL are rare at such large distances from the Galactic center. A more likely explanation for this group is that they are actually AC stars associated with Fornax. As discussed previously, the light curves of RRL and ACs are easily confused. In this case, these stars reside in the region of the Fornax color-magnitude diagram (CMD) that is expected for AC stars, and such objects are known to exist in Fornax (Bersier & Wood 2002; Greco et al. 2005), although none of these stars match previously known variables. The spatial distribution of these candidates is shown in Figure 10.

## 5.4. Cepheids in Local Group Galaxies and Beyond

The Local Group galaxies Phoenix, IC 1613, and Tucana are located within the DES footprint at approximate distances of 415 kpc, 755 kpc, and 887 kpc, respectively (McConnachie 2012). Both Phoenix and IC 1613 are spatially coincident with overdensities of objects in our catalog, while Tucana is not.

In the case of Phoenix, we found three RRL candidates within 5' from the center of the galaxy, which has a half-light radius $r_h = 3.76'$, with $22.3 < g < 22.7$. True RRL stars in this galaxy are expected to be ~1 mag fainter. Since it is unlikely to find distant halo field stars in the line of sight of Phoenix, we
believe instead that these may be misclassified AC stars; such objects are known to exist in this galaxy (Gallart et al. 2004).

In the case of IC1613, we found a group of 11 RRL candidates with $22 < g < 23$ and within $3 r_h$ from the center of the galaxy. Given the larger distance of IC1613, this range of observed magnitudes is appropriate for classical Cepheids, as even ACs will be too faint for our catalog in this galaxy. The sample of classical Cepheids in IC1613 by Bernard et al. (2010) shows numerous stars in this magnitude range, a few of which have periods as short as 0.6 days. Udalski et al. (2001) found 138 Cepheids within the central $14'.2 \times 14'.2$ region of IC1613. We thus suspect some of the 11 objects in our catalog may be short-period classical Cepheids in this galaxy, even though we do not find any matches with previous catalogs. It is also possible that at these faint magnitudes crowding and the misclassification of background galaxies may lead to spurious detections.

Other nearby galaxies (beyond the Local Group) in the DES footprint are ESO 410-G005, ESO 294-G010, NGC 55, NGC 300, and IC 5152. Since these galaxies have distance moduli of ~26.5 mag (McConnachie 2012), it is unlikely that our analysis would detect any variable stars. Indeed, no overdensities in the DES Y6 RRab catalog are associated with these objects.

## 5.5. Ultrafaint Dwarf Galaxies and Globular Clusters

RRL in ultrafaint dwarf galaxies (UFDs) can provide independent estimates of the distances to these systems, which are particularly important given their low surface brightnesses and sparsely populated CMDs. Recent censuses of RRL in UFDs can be found in Martínez-Vázquez et al. (2019) and Vivas et al. (2020a). We use the DES Y6 RRab candidate catalog to search for RRL in the vicinities of 19 UFDs in the DES footprint and find strong evidence of variables associated with Eridanus II, as well as tentative identifications in Cetus III and Tucana IV.

Eridanus II is among the most luminous ($M_V \sim -7.1$) and distant ($D \sim 366$ kpc) UFDs (Crnojević et al. 2016). We search for RRL in our catalog within $10 r_h = 10'.1$ and find five stars with a very narrow range of distance moduli, $\mu_0 = 22.45 \pm 0.04$ mag. All variables are tightly concentrated within $2.4 r_h = 3'.4$ and are confirmed to be RRL in this system by C. E. Martínez-Vázquez et al. (2021, in preparation). Their mean period is 0.663 days, consistent with the Oo II group found in other UFDs (Martínez-Vázquez et al. 2019). Given the low completeness of our survey at these faint magnitudes, the RRL population in this system should be significantly larger. There are two additional stars in our search area with $g = 22.1 - 22.3$, or ~0.8 mag brighter than the RRL. These may be anomalous Cepheids in Eridanus II.

Cetus III is somewhat closer ($D \sim 251$ kpc) and has an absolute magnitude of only $M_V \sim -2.5$ (Homma et al. 2018). Little is known about this small ultrafaint dwarf, and spectroscopic confirmation of its nature is not yet available. We found one RRab candidate located $1'$ from its center (~1 $r_h$) with $\mu = 21.33$ mag, which would put the galaxy at 185 kpc from the Sun, somewhat closer than the aforementioned estimate. However, the difference is consistent with the bias introduced by our adoption of the mean halo metallicity for our RRL template and the lower metallicity expected for this UFD. Further studies of Cetus III and this RRL candidate are needed to confirm their association.

We identify one RRL in Tucana IV with $g = 18.78$, consistent with spectroscopically confirmed HB members of this galaxy (Simon et al. 2020). The variable, however, is located at $56'$ from the center, equivalent to $6 r_h$. Thus, the physical association is unclear, although the Gaia DR2 proper motion is consistent. Another possibility is that this may be an extratidal star, similar to those seen in Tucana III, Eridanus III, and Reticulum III (Vivas et al. 2020a).

We also note that we recover known RRL in other UFDs, including star V1 in Tucana II (Vivas et al. 2020a), V1 in Phoenix II, and V2 in Grus I (Martínez-Vázquez et al. 2019). However, since more comprehensive analyses of these objects exist in the referenced publications, we do not examine them in detail here.

Finally, there are several globular clusters within the DES footprint. RRL in Milky Way globular clusters closer than 10 kpc often saturate the DES images. Nonetheless, we recover star V21 in M2 (NGC 7089; Clement et al. 2001, 2017 July version), which is one of the 23 RRab known in that cluster. On the other hand, the LMC globular cluster Reticulum is a better target for the range of magnitudes of our catalog. We recover 14 of its 22 known RRab (Kuehn et al. 2013), a 63% recovery rate. We find two additional RRab candidates spatially coincident with the cluster that have the appropriate magnitude and proper motions (from Gaia DR2) to be members. The mean distance modulus of these 16 variables is $\mu = 18.36 \pm 0.04$ mag, in agreement with the value obtained by Kuehn et al. (2013) of $18.40 \pm 0.20$ mag. The good agreement in these estimates is due to the similar metallicity of the RRL in the cluster ([Fe/H] $\sim -1.6$) and the mean metallicity of the Milky Way halo adopted for our templates. Additionally, there is a third new possible member that has consistent magnitude and PM, but it is located farther away, at $14'$ from the center of the cluster.

## 5.6. Halo Density Profile

The structure of the Milky Way stellar halo encapsulates information about the formation and evolution of our Galaxy (e.g., Johnston et al. 2002). Several lines of observational evidence suggest that the halo density profile exhibits a break at Galactocentric distances of 20–35 kpc, with star counts falling off more rapidly beyond this radius (e.g., Watkins et al. 2009; Sesar et al. 2010; Deason et al. 2011; Sesar et al. 2011; Zinn et al. 2014; Pila-Díez et al. 2015; Xue et al. 2015; Pieres et al. 2020). Such a broken-power-law profile may be produced through the accretion of a massive satellite galaxy (Bullock & Johnston 2005; Deason et al. 2013, 2018), which corroborates recent claims of the Gaia–Enceladus merger (Belokurov et al. 2018; Helmi et al. 2018). While quantitative estimates vary by analysis, studies across a wide range of stellar tracers find that shallower power-law slopes ($n_1 \sim -2$ to $-3.5$) are preferred in the inner region of the halo, while steeper values ($n_2 \sim -3.8$ to $-5.8$) are preferred at larger distances (see Pila-Díez et al. 2015 for a recent compilation). RRL have provided an important probe of the stellar density profile, with evidence for the broken-power-law model first claimed by Saha (1985) and more recently by Zinn et al. (2014) and Medina et al. (2018). The deep, wide-area catalog of DES Y6 RRab candidates offers an excellent opportunity to measure the density profile of the Milky Way stellar halo over a wide region of the southern sky.

We estimate the halo density profile by first removing RRab candidates around the LMC, SMC, Fornax, and Sculptor (see
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Left: period-amplitude diagram of DES Y6 RRab candidates (black) and high-confidence candidates associated with the Oosterhoff I sequence and possessing amplitude variations >0.75 mag (blue). Right: density of RRab candidates as a function of heliocentric distance. The gray points show the observed distribution of all RRab candidates, while the black curve corrects for the detection efficiency estimated in Section 5.1. The blue curve shows the density of high-confidence RRab associated with the Oo-I sequence in the left panel, corrected by detection efficiency (which is >95% out to 150 kpc).</figcaption></figure>

Appendix C), which yields a sample of 3603 candidate RRab. We group these stars into 51 bins of heliocentric distance 9 < D < 400 kpc, and we calculate the number density by correcting for the detection efficiency of our catalog (Figure 11). We find that the density profile of RRab exhibits a break at a heliocentric distance of D ~ 25 kpc, with an excess relative to a simple broken-power-law model at heliocentric distances of 80 kpc < D < 300 kpc. The excess of candidates at large distances is not explained by a change in detection efficiency, which is found to be decreasing monotonically with distance and is a fractionally smaller effect than the observed excess (Section 5.1). To investigate this excess in more detail, we select a subpopulation of high-confidence RRab candidates close to the Oosterhoff I sequence with a measured amplitude change of >0.75 mag. We find that the excess is greatly reduced in the high-confidence subpopulation, suggesting that contamination from faint sources may be responsible for the excess. To further investigate possible contamination, we cross-match RRab candidates in the S82 region with spectroscopically confirmed QSOs from SDSS DR7 (Schneider et al. 2010) and SDSS DR16 (Lyke et al. 2020). We find that the QSO contamination of our catalog in the S82 region is ~2%; however, we find that all contaminating QSOs are fit with distance moduli 20 < μ < 22 (100 kpc < D < 250 kpc), leading to a contamination rate of ~16% in the same distance range as the observed excess.

While suggestive, this contamination is still far less than would be necessary to account for the observed excess. Furthermore, we note that the RRab candidates in this distance range are not uniformly distributed over the footprint, as would be expected from extragalactic contamination and the uniform recovery efficiency estimated from our RRL simulations. Rather, RRab candidates in this distance range are preferentially distributed at α ~ 0, coincidentally close to where the Magellanic Stream crosses the DES footprint.

Large-scale anisotropies in the halo RRL distribution have been claimed by Iorio et al. (2018) using a combined catalog of Gaia+2MASS RRL. Boubert et al. (2019) provided supporting evidence for this structure using observations from the Catalina Surveys (e.g., Torrealba et al. 2015; Drake et al. 2017) and the

sample of RR Lyrae variables identified in PS1 (Hernitschek et al. 2016). While this structure is significantly closer than the excess observed here (Galactocentric distance of ~20 kpc), Boubert et al. (2019) claim a Magellanic origin for this overdensity, which could extend to greater distances where the infall track of the Magellanic Clouds crosses the DES footprint. Further investigation of these distant candidates is necessary to better understand possible anisotropies in the RRL distribution at distances ≳80 kpc.

Given the significant uncertainties in the contamination rate and possible anisotropies at heliocentric distances ≳80 kpc, we constrain our study of the Milky Way halo to smaller distances, where we estimate that our completeness is ≳75%. To measure the Milky Way stellar halo density, we transform each of our RRab candidates into Galactocentric coordinates, (x, y, z), assuming that the solar Galactic center distance is 8.178 kpc (Gravity Collaboration et al. 2019). We further calculate the elliptical Galactocentric radius, $r_e = \sqrt{x^2 + y^2 + (z/q)^2}$. Following Faccioli et al. (2014), we perform our fit assuming a fixed halo flattening of $q = 0.7$ (e.g., Sesar et al. 2011). We group candidate RRab into 41 bins in elliptical Galactocentric radius 9 kpc < $r_e < 100$ kpc, and we fit the halo profile with an elliptical broken power law following the description of Pila-Díez et al. (2015):

$$ \rho(x, y, z) = \rho_0 \begin{cases} r_e^{n_1}, & r_e < R_0 \\ r_e^{n_2} \cdot R_0^{n_1-n_2}, & r_e \ge R_0, \end{cases} \quad (5) $$

where $ρ_0$ is the density normalization, $n_1$ is the inner power-law index, $n_2$ is the outer power-law index, and $R_0$ is the break radius (in elliptical Galactocentric coordinates). We perform a binned Poisson maximum-likelihood fit for the observed counts in each bin, k, as a function of the model parameters, $θ = \{ρ_0, R_0, n_1, n_2\}$. Our likelihood analysis is described in more detail in Appendix C. Within each bin, we correct the observed number of RRL for the detection completeness of our catalog estimated from simulated RRab described in Section 5.1. We fit the model parameters using a Markov Chain Monte Carlo ensemble sampler (emcee; Foreman-Mackey et al. 2013), and
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Density of RRab candidates binned in elliptical Galactocentric radius ($q = 0.7$). The black markers show the full RRab candidate sample, while blue markers show high-confidence RRab associated with the Oosterhoff I sequence and having measured amplitude of variations of >0.75 mag. Errors represent the square root of the number of stars in each bin and the bin width.</figcaption></figure>

we report the median, 16th, and 84th percentiles of the marginalized posterior distributions of each parameter in Table 2. We perform this analysis for both the full DES Y6 RRab candidate sample and the high-confidence RRab associated with the Oosterhoff I sequence (Figure 12).

In Figure 13 we compare the best-fit parameters from our two RRab candidate samples to broken-power-law fits to the halo in the literature. We generally find that our best-fit break radius of ~32 kpc is slightly larger than many other analyses; however, this could be brought into better agreement through a smaller halo flattening. Our measured inner slope values are consistent within $2\sigma$ of each other and are broadly consistent with other values in the literature. Due to the saturation threshold of the DES images, fits for $n_1$ are largely driven by RRL with $r_e \gtrsim 14$ kpc and are highly correlated with the overall normalization parameter, $\rho_0$. Due to the large area and sensitivity of DES, we have several thousand RRab candidates in the range $30 \text{ kpc} < r_e < 100 \text{ kpc}$. This allows the DES data to provide tight constraints on the outer power-law slope, $n_2$. The value of $n_2 = -5.42_{-0.14}^{+0.13}$ measured using DES RRL candidates is steeper than that measured by many other tracers and surveys (e.g., Sesar et al. 2011; Xue et al. 2015; Pieres et al. 2020), but it is consistent with RRL measurements by Sesar et al. (2010) and Zinn et al. (2014).

## 5.7. Search for New Satellite Galaxies

Clusters of RRL can be used to identify faint satellite galaxies that may have evaded detection by other methods. In particular, Baker & Willman (2015) argued that the combination of three-dimensional information and the sparsity of halo RRL at distances >50 kpc make RRL particularly useful for identifying Milky Way satellites with half-light radii $r_h > 500$ pc residing in the outer halo. The RRL catalogs from Gaia DR2 have demonstrated the viability of this search technique through the discovery of the ultradiffuse satellite Antlia II (Torrealba et al. 2019) and the detection of several candidate stellar streams (e.g., Mateu et al. 2018).

We search for satellite galaxies coincident with each of the DES Y6 RRab candidates beyond the masked regions around

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Comparison of best-fit broken-power-law halo parameters from this work (S21) and those collected from the literature. Circle markers represent measurements from RRL, squares represent measurements from main-sequence turnoff stars, while triangles represent measurements from giant branch stars. Parameters are the break radius (left), the inner power-law index (center), and the outer power-law index (right). References are W09 (Watkins et al. 2009), S10 (Sesar et al. 2010), S11 (Sesar et al. 2011), D11 (Deason et al. 2011), F14 (Faccioli et al. 2014), Z14 (Zinn et al. 2014), P15 (Pila-Díez et al. 2015), X15 (Xue et al. 2015), D16 (Das et al. 2016), M18 (Medina et al. 2018), and P20 (Pieres et al. 2020).</figcaption></figure>

<table><caption>**Table 2**
Broken-power-law Halo Density Parameters</caption><thead><tr><th>Sample</th><th>ρ<sub>0</sub><br>(kpc<sup>-3</sup>)</th><th>R<sub>0</sub><br>(kpc)</th><th>n<sub>1</sub></th><th>n<sub>2</sub></th></tr></thead><tbody><tr><td>Full</td><td>990<sup>+300</sup><sub>-230</sub></td><td>32.1<sup>+1.1</sup><sub>-0.9</sub></td><td>-2.54<sup>+0.09</sup><sub>-0.09</sub></td><td>-5.42<sup>+0.13</sup><sub>-0.14</sub></td></tr><tr><td>Oo-I</td><td>240<sup>+110</sup><sub>-80</sub></td><td>31.3<sup>+1.3</sup><sub>-1.0</sub></td><td>-2.32<sup>+0.14</sup><sub>-0.13</sub></td><td>-5.59<sup>+0.17</sup><sub>-0.19</sub></td></tr></tbody></table>

**Note.** The halo density profile fit assumes a fixed halo flattening of $q = 0.7$.

the LMC, SMC, Fornax, and Sculptor. Our algorithm is based on a simple binned Poisson likelihood, which is qualitatively similar to searches for resolved stellar populations (Bechtol et al. 2015; Drlica-Wagner et al. 2015), but optimized to the detection of satellites with half-light radii $r_h \gtrsim 500$ pc. At the location and distance of each RRab in our masked catalog, we define a search cylinder with a fixed radius of $2r_h = 1$ kpc. The depth of our search cylinder is calculated from the quadrature sum of $2r_h$ and the systematic uncertainty on the measured distance to the RRL, $2\sigma(\mu) = 0.2$ mag. Each cylinder is expected to contain ~90% of the RRL population of a satellite with $r_h = 500$ kpc centered at the search location.

We determine the local expected density of field RRL, $\rho_f$, from a cylindrical annulus centered on our search location with inner and outer radii of $4r_h$ and $8r_h$, respectively. In cases where no RRL are contained within our background annulus, we assume the global background rate as estimated from the field density at that distance (Figure 11). We multiply the background density by the volume of our search cylinder to predict the expected number of field RRL within our search region, $\lambda$.

We calculate the significance of a putative satellite at each search location (i.e., at the location of each RRab candidate in our catalog) as the Poisson probability of detecting $k$ or more
<figure><img src="image_16.png" /><figcaption>Figure 14. Detection efficiency of our search for Milky Way satellites using DES Y6 RRab candidates. Detection efficiency ranges from 0% (blue) to 100% (yellow) and is shown as a function of azimuthally averaged, projected 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 satellite galaxies are indicated in gray. The light blue dashed line shows the 50% detectability contour for RRL searches at distances >50 kpc predicted by Baker & Willman (2015), while the gray dashed line shows an analytic approximation to the 50% detectability contour from isochrone-matched filter searches (Drlica-Wagner et al. 2020).</figcaption></figure>

RRL given an expectation of $\lambda$:

$$P(k|\lambda) = \frac{\lambda^k e^{-\lambda}}{k!} \qquad (6)$$

To select candidate satellites, we apply a significance threshold of $p < 3 \times 10^{-7}$, which corresponds to a one-sided Gaussian significance of $5\sigma$. We also require that satellite candidates consist of at least three RRab candidates.

We quantify the sensitivity of our search using a suite of $10^5$ simulated satellite galaxies generated by Drlica-Wagner et al. (2020). These satellites span a range of stellar mass, heliocentric distance, size, ellipticity, and position angle (see Table 1 of Drlica-Wagner et al. 2020). The simulated satellites are distributed uniformly over the DES footprint and uniformly in distance modulus. For each simulated satellite, we predict the expected number of RRL as a function of $M_V$ using a fit to the RRL population of Milky Way satellites provided in Equation (4) of Martínez-Vázquez et al. (2019):

$$\log(N_{\text{RRL}}) = -0.29 \times M_V - 0.80 \qquad (7)$$

To predict the number of RRab, we multiply the number of RRL by the fraction of RRab, $N_{\text{RRab}} = f_{\text{ab}}N_{\text{RRL}}$, where $f_{\text{ab}} = 0.71$ for Milky Way satellites (see Table 6 of Martínez-Vázquez et al. 2017). The expected number of RRab observed by DES is corrected for the detection efficiency of our search and fitting procedure, which depends on the distance of the simulated satellite. The spatial distribution of simulated RRab was drawn from an elliptical Plummer profile (Plummer 1911). Distances were assigned from a Gaussian distribution centered on the distance of the simulated satellite matched to the 3D half-light radius of the satellite. An additional systematic Gaussian scatter in distance modulus, $\Delta(\mu) = 0.1$ mag, was applied to each simulated RRab. We inject each simulated satellite into the DES Y6 RRab candidate catalog individually and attempt to recover it with our search algorithm. The

resulting detection efficiency as a function of $M_V$ and physical half-light radius $r_h$ is shown in Figure 14.

As expected, our search is more sensitive than isochrone-matched filter searches (Drlica-Wagner et al. 2020, dashed line) for satellites with large sizes. Our RRL satellite search is significantly less sensitive for compact satellites, due to the large assumed kernel ($2r_h = 1$ kpc). At small heliocentric distances, this large search kernel leads to an expected background contribution from halo RRL that is comparable with the RRL signal from a satellite with $M_V \sim -6$. At larger distances, the density of halo RRL decreases, and the sensitivity of our search increases until RRL detection efficiency starts to dominate. As a test, we rerun these simulations using a kernel matched to the true size of each simulated satellite, and we find a significant improvement in the sensitivity for small satellites. However, the isochrone-matched filter searches remain more sensitive for satellites with $M_V > -5$, due to the small number of RRab expected from these satellites.

With the sensitivity of our search characterized on simulations, we apply our search to the DES Y6 candidate RRab catalog. We find three satellite candidates that pass our detection criteria of significance >5$\sigma$ and $N_{\text{RRab}} > 3$. The characteristics of each RRab candidate member associated with these satellite candidates are listed in Table 3.

One of these candidates (SubId = 1) is associated with the known satellite Eridanus II, located at a distance of 360 kpc (Section 5.5). Based on the simulations described above (Figure 14), we expect our search to be ~60% efficient for satellites with the size, luminosity, and distance of Eridanus II. This could suggest that Eridanus II may have a larger-than-expected number of RRL.

A second overdensity (SubId = 2) consists of three RRab candidates located at $(\alpha, \delta) \sim 349.8, -55.6$ at a distance of 59.5 kpc. This overdensity is in the same region of the sky and at roughly the same distance as the Tucana II dwarf galaxy,
<table><caption>**Table 3**
Halo Substructure Candidates</caption><thead><tr><th>SubID</th><th>Star ID</th><th>α<br>deg, J2000</th><th>δ<br>mag</th><th>μ<br>kpc</th></tr></thead><tbody><tr><td>1</td><td>1410940364</td><td>56.10118</td><td>-43.50492</td><td>22.5</td><td>316.2</td></tr><tr><td>1</td><td>1410941177</td><td>56.07512</td><td>-43.51257</td><td>22.4</td><td>301.8</td></tr><tr><td>1</td><td>1410942171</td><td>56.09949</td><td>-43.52116</td><td>22.4</td><td>305.1</td></tr><tr><td>1</td><td>1410943700</td><td>56.07939</td><td>-43.53486</td><td>22.5</td><td>314.1</td></tr><tr><td>1</td><td>1410944738</td><td>56.01504</td><td>-43.54419</td><td>22.4</td><td>307.6</td></tr><tr><td>2</td><td>994558814</td><td>348.37430</td><td>-55.77342</td><td>18.7</td><td>55.0</td></tr><tr><td>2</td><td>1002871838</td><td>349.83519</td><td>-55.33808</td><td>18.9</td><td>60.5</td></tr><tr><td>2</td><td>1007027294</td><td>351.32392</td><td>-55.76393</td><td>19.0</td><td>63.0</td></tr><tr><td>3</td><td>1533696855</td><td>77.66305</td><td>-38.90769</td><td>17.0</td><td>24.6</td></tr><tr><td>3</td><td>1539669959</td><td>77.95832</td><td>-38.64577</td><td>17.1</td><td>25.8</td></tr><tr><td>3</td><td>1548634844</td><td>80.17935</td><td>-39.58687</td><td>16.8</td><td>23.2</td></tr><tr><td>3</td><td>1616940370</td><td>77.32652</td><td>-37.11542</td><td>16.7</td><td>21.9</td></tr></tbody></table>

$D = 57.5$ kpc; however, it is separated from Tucana II by ~5 deg on the sky. Interestingly, the Gaia DR2 proper motions of the RRab candidates in this overdensity are similar to the proper motions of confirmed members of Tucana II, albeit with large uncertainties (Figure 15). It has been suggested that the diffuse structure of Tucana II could be an indication of tidal disruption (e.g., Bechtol et al. 2015). Belokurov & Koposov (2016) showed evidence for extended stellar structure around Tucana II using BHB stars (i.e., their “S2a” cloud at $\mu \sim 18.8$), and Chiti et al. (2020) used Gaia proper motions combined with photometric metallicities to identify probable members >1 deg from Tucana II. Further observations are required to confirm an association between this overdensity of RRab candidates and the Tucana II dwarf galaxy.

The third candidate (SubId = 3) consists of four RRab candidates located at $(\alpha,\delta) \sim 78.3, -38.6$ at a distance of 23.8 kpc. This is close in projection to the globular cluster NGC 1851, but at a larger distance. Shipp et al. (2018) found evidence for extended tidal features surrounding this cluster, but the proper motions of the stars associated with this candidate structure are not aligned with the motion of NGC 1851. We followed the procedure of Shipp et al. (2018) to use the DES Y6 coadd object catalog to search for correlated overdensities of main-sequence turnoff stars associated with our candidate substructures using isochrone-filtered stellar density maps. However, we do not detect any previously unknown overdensities in turnoff stars, suggesting that radial velocities will be critical for confirming the second two candidate structures.

## 6. Conclusion

We have used various statistical techniques to derive a catalog of candidate RRab detected in six years of deep, wide-area imaging from DES. The DES Y6 data offers significant improvements on prior results using only three years of data (Stringer et al. 2019), resulting in a catalog of 6971 RRab candidates. At the bright end, our catalog has significant overlap with surveys such as Gaia DR2, Pan-STARRS, and the Catalina surveys, providing strong evidence of the effectiveness of this method. We publicly release the best-fit parameters and light curves for our RRL candidates.

We find good agreement in the measured properties of our RRab candidates when matched against external catalogs from

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Proper motions of RRab candidate stars associated with satellite candidate 2 (Table 3). Stars associated with this satellite candidate are shown in blue, spectroscopically confirmed member stars and proper motion candidate members from the Tucana II satellite are shown in green, the systemic motion of Tucana II is shown in orange, and foreground stars are shown in red (Pace & Li 2019).</figcaption></figure>

the MCs, Fornax, and Sculptor. In addition, we recover RRL detected in the ultrafaint dwarfs Tucana II, Phoenix II, and Grus I (Martínez-Vázquez et al. 2019). DES extends significantly deeper than any of these surveys, allowing us to detect RRab candidates out to a distance of ~335 kpc, about 1 mag deeper than Stringer et al. (2019). Indeed, we discover a group of five RRab candidates associated with the distant ultrafaint dwarf galaxy Eridanus II. We also report tentative RRL members of the ultrafaint systems Cetus III (Homma et al. 2018) and Tucana IV (Drlica-Wagner et al. 2015). We fit the stellar density profile of the Milky Way halo in the range of elliptical Galactocentric distances $9 < r_e < 100$ kpc. Assuming a halo flattening of $q = 0.7$, we find that the halo is well fit by a broken-power-law model with a break radius of $R_0 = 32.1_{-0.9}^{+1.1}$, an inner slope of $n_1 = -2.54_{-0.09}^{+0.09}$, and an outer slope of $n_2 = -5.42_{-0.14}^{+0.13}$. These values agree between analyses of the full RRab candidate sample and a high-confidence sample of RRab candidates with large measured amplitude variations associated with the Oosterhoff I sequence. We further use our catalog of DES Y6 RRab candidates to search for halo substructures, with characteristic sizes of $r_h \sim 500$ pc. This search confidently detects the Eridanus II dwarf galaxy and two other candidate overdensities that are not confidently associated with known halo substructures.

RRL have long been recognized as a powerful probe of the Milky Way’s outer stellar halo. However, it is only recently that surveys have been able to combine wide area coverage, deep imaging, and sufficient cadence to confidently identify RRL at distances >100 kpc. While DES was not optimized to search for RRL, it nonetheless provides an exceptional catalog of RRab candidates at distances beyond what is achievable by surveys on smaller telescopes (i.e., Gaia, Catalina, and PS1). In the near future, the Vera C. Rubin Observatory Legacy Survey of Space and Time will provide hundreds of observations over the entire southern sky, promising to provide a high-quality catalog of RRL extending to the edge of the Milky Way halo.

K.M.S., L.M.M., and P.S.F. acknowledge support from the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University. This material is based upon work that was supported by the Fermilab Visiting Scholars Award
Program of the Universities Research Association. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. A.B.P. acknowledges support from NSF grant AST-1813881.

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

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NFS's NOIRLab, 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 at NSF's NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS website is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under program IDs 177.A-3016, 177.A-3017, and 177.A-3018, and on data products produced by Target/OmegaCEN, INAF-OACN, INAF-OAPD, and the KiDS production team, on behalf of the KiDS consortium. OmegaCEN and the KiDS production team acknowledge support by NOVA and NWO-M grants. Members of INAF-OAPD and INAF-OACN also acknowledge the support from the Department of Physics & Astronomy of the University of Padova and from the Department of Physics of Univ. Federico II (Naples).

Facilities: Blanco (DECam), Gaia.

Software: Astropy (Astropy Collaboration et al. 2018), HEALPix (Górski et al. 2005), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SourceExtractor (Bertin & Arnouts 1996) scikit-learn (Pedregosa et al. 2012), SciPy (Virtanen et al. 2020).

## Appendix A: Variability Statistics

Here, we describe the variability statistics used to select objects for template fitting. A glossary of terms is provided in Table 4. For each object, we multiplied the reported photometric uncertainties by scaling factors based on the best-fit values from Equation (2). We use these rescaled errors in the calculation of all variability statistics. To avoid confusion, we will refer to the error-weighted mean as $\bar{m}$ and the median as med($m$).
<table><caption>**Table 4**

Glossary of Terms</caption><thead><tr><th>Abbreviation</th><th>Full Name</th><th>Bands Used</th><th>Ref.</th></tr></thead><tbody><tr><td>IQR<sub>b</sub></td><td>Interquartile range for band <i>b</i></td><td>g, r, i, z</td><td>1</td></tr><tr><td>J<sub>b</sub></td><td>Stetson's <i>J</i> statistic for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td>2</td></tr><tr><td>MAD<sub>b</sub></td><td>Median absolute deviation for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td>1</td></tr><tr><td>NAPD<sub>b</sub></td><td>N absolute pairwise distances for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td>1</td></tr><tr><td>σ<sup>2</sup><sub>NXS,b</sub></td><td>Normalized excess variance for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td>3,4</td></tr><tr><td>Q<sub>n,b</sub></td><td><i>k</i>th value of absolute pairwise distances for band <i>b</i></td><td>g, r, i, z</td><td>5</td></tr><tr><td>Δ(mag)<sub>b</sub></td><td>Magnitude range for band <i>b</i></td><td>g, r, i, z</td><td></td></tr><tr><td>χ<sup>2</sup><sub>ν,b</sub></td><td>Reduced chi-squared statistic for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td></td></tr><tr><td>χ̃<sup>2</sup><sub>ν,b</sub></td><td>Reduced median chi-squared statistic for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td></td></tr><tr><td>RoMS<sub>b</sub></td><td>Robust median statistic for band <i>b</i></td><td>g, r, i, z, (grizY)</td><td>6</td></tr><tr><td>S<sub>n,b</sub></td><td>Median value of absolute pairwise distances for band <i>b</i></td><td>g, r, i, z</td><td>5</td></tr><tr><td>W_Range<sub>b</sub></td><td>Weighted magnitude range for band <i>b</i></td><td>g, r, i, z</td><td>A</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,b</sub>)<sub>th</sub></td><td>Common logarithm of the reduced chi-squared above threshold for band <i>b</i></td><td>g, r, i, z</td><td>A</td></tr><tr><td>rss<sub>ν,j</sub></td><td>Residual sum of squares from template curves. <i>j</i> ∈ [0, 2] denotes the rank order minima of this function</td><td>(grizY)</td><td>S19</td></tr><tr><td>μ</td><td>Distance modulus of the best-fitting template</td><td>...</td><td>S19</td></tr><tr><td>A<sub>g</sub></td><td>g-band amplitude for the best-fitting template</td><td>...</td><td>S19</td></tr><tr><td>amp<sub>j</sub></td><td>g-band amplitude for the <i>j</i>th best-fitting template, <i>j</i> ∈ [0, 2]</td><td>...</td><td>S19</td></tr><tr><td>P</td><td>Best-fitting template period in units of days</td><td>...</td><td>S19</td></tr><tr><td>φ</td><td>Template estimated phase offset in units of days</td><td>...</td><td>S19</td></tr><tr><td>D(Oo-I)<sub>j</sub></td><td>Distance from the Oosterhoff I sequence for <i>j</i>th best-fit period and g-band amplitude, <i>j</i> ∈ [0, 2]</td><td>...</td><td>S19</td></tr><tr><td>D(Oo-II)<sub>j</sub></td><td>Distance from the Oosterhoff II sequence for <i>j</i>th best-fit period and g-band amplitude, <i>j</i> ∈ [0, 2]</td><td>...</td><td>S19</td></tr><tr><td>D(Oo-int)<sub>j</sub></td><td>Distance from the intermediate Oosterhoff sequence for <i>j</i>th best-fit period and g-band amplitude, <i>j</i> ∈ [0, 2]</td><td>...</td><td>S19</td></tr><tr><td>RFi_score</td><td>Output score from the <i>i</i>th random forest classifier, <i>i</i> ∈ [1, 2]</td><td></td><td></td></tr></tbody></table>

**Note.** References: (1) Sokolovsky et al. (2017), (2) Stetson (1996), (3) Nandra et al. (1997), (4) Simm et al. (2015), (5) Rousseeuw & Croux (1993), (6) Enoch et al. (2003), (A) Appendix A, (S19) Stringer et al. (2019).

In addition to the reduced median chi-squared, $\tilde{\chi}_{\nu,b}^2$ (Equation (1)), we calculate the reduced chi-squared from the mean magnitude in each band, $\chi_{\nu,b}^2$, and its common logarithm

$$ \log_{10}(\chi_{\nu,b}^2) = \log_{10}\left(\frac{1}{N_b - 1} \sum_{i=1}^{N_b} \frac{(m_{i,b} - \bar{m}_b)^2}{\sigma_{i,b}^2}\right). \quad (A1) $$

As these quantities use the rescaled errors determined in Section 2.2, the nonvarying sources have a distribution centered around zero, with positive outliers denoting true variable objects.

We also measure the range in magnitude, which we call $\Delta(mag)_b$, in each single-band light curve to relay information about the amplitude of an object's variation. The uncertainties on the maximum and minimum magnitudes used to calculate this quantity are recorded as well.

Another proxy for variability amplitude is the normalized excess variance ($\sigma_{NXS}^2$). This statistic was first defined by Nandra et al. (1997) as

$$ \sigma_{\text{NXS}}^2 = \frac{1}{N\bar{m}^2} \sum_{i=1}^{N} [(m_i - \bar{m})^2 - \sigma_i^2]. \quad (A2) $$

Although this metric is commonly used for X-ray analyses of active galactic nuclei, it has successfully been deployed on the sparsely sampled Pan-STARRS 3π optical light curves (Simm et al. 2015).

We measure the overall scatter of a light curve using the median absolute deviation (MAD):

$$ \text{MAD} = \text{med}(|m_i - \text{med}(m)|). \quad (A3) $$

MAD is slightly more stable in the presence of outliers than the standard deviation as it does not amplify the effects of an outlier by squaring it. This metric should be sensitive to repeated variations; however, a real RRL with observations sampled at close to the same phase value will not appear variable in MAD.

The robust median statistic (RoMS) is a more robust analog of $\chi_{\text{red}}^2$, which is less sensitive to bias in the presence of non-Gaussian uncertainties. This metric was first defined in Enoch et al. (2003) as

$$ \text{RoMS} = \frac{1}{N-1} \sum_{i=1}^{N} \frac{|m_i - \text{med}(m)|}{\sigma_i}. \quad (A4) $$

In a single band, RoMS tends toward values of one for nonvarying sources.

Alternative measures of deviation $S_n$ and $Q_n$ were first proposed by Rousseeuw & Croux (1993). These metrics seek to measure the midpoint of a data set like MAD, but do not rely on a central reference value and are thus better estimators for asymmetric distributions. To account for the scatter caused by the large uncertainties of faint observations, we will apply these statistics to the pairwise distances between all individual observations divided by their combined uncertainties. Here, $S_n$ measures the median of the median of these weighted pairwise differences and is defined as

$$ S_n = 1.1926 \text{ med}_i \left( \text{med}_j \left( \frac{|m_i - m_j|}{\sqrt{\sigma_{m_i}^2 + \sigma_{m_j}^2}} \right) \right) \quad (A5) $$
where $m_i$ and $m_j$ are separate observations. A similar metric $Q_n$ is defined as

$$Q_n = 2.2219 \left( \frac{|m_i - m_j|}{\sqrt{\sigma_{m_i}^2 + \sigma_{m_j}^2}} ; i < j \right)_k \quad (A6)$$

Note that $Q_n$ records the kth value in the weighted absolute pairwise differences between all of the observations, and k is the binomial coefficient $\binom{h}{2}$ with $h = \lfloor n/2 \rfloor + 1$ for n observations. Thus, $Q_n$ records roughly the midpoint of these values. Another way to measure this is to calculate quantiles for the differences in separate observations. We use the 90% quantile values of these error-weighted N absolute pairwise distances (NAPD) as a proxy for the weighted range that is less sensitive to outliers.

Stetson's J variability index (Stetson 1996) has been widely used to identify pulsating variables such as Cepheids and RRL. It builds upon the Welch-Stetson variability index I (Welch & Stetson 1993), which measures the correlation between n subsequent pairs of observations. Stetson's J index measures this correlation using single observations as well as pairs, so we can apply it to both the single-band and multiband light curves:

$$J = \frac{\sum_{k=1}^{n} w_k \operatorname{sgn}(P_k) \sqrt{|P_k|}}{\sum_{k=1}^{n} w_k} \quad (A7)$$

where sgn is the sign function and

$$P_k = \begin{cases} \text{pair} & \left( \sqrt{\frac{n_v}{n_v - 1} \frac{v_i - \bar{v}}{\sigma_{v_i}}} \right) \left( \sqrt{\frac{n_b}{n_b - 1} \frac{b_i - \bar{b}}{\sigma_{b_i}}} \right) \\ \text{single} & \frac{n_v}{n_v - 1} \left( \frac{v_i - \bar{v}}{\sigma_{v_i}} \right)^2 - 1 \end{cases} \quad (A8)$$

for pairs of observations in bands v and b and for single observations. For $\bar{v}$ and $\bar{b}$, we calculate the weighted mean for the observations in that band. Sources with purely noisy light curves have values close to zero.

The DES observing strategy often yields sequences of two to three exposures taken of the same field in different filters, with the filters chosen according to the seeing, Moon phase, and number of missing observations in that field (see Figure 3 in Diehl et al. 2016). Thus, to take advantage of this, we group together any observations taken within an hour. Any group with three observations abc within the same hour was treated as three pairs ab, bc, ac, each weighted with a factor of 2/3 as prescribed in Stetson (1996). All other groups with one or two observations are weighted equally with a factor of 1. For this measurement, we excluded the generally noisier observations in Y.

We defined two new variability metrics intended to avoid overly penalizing faint objects, which can manifest measurement variability beyond that expected from their statistical uncertainties. The first of these metrics was calculated from the difference between the minimum and maximum measured magnitude in each band, $mag_b$, weighted by the photometric

<figure><img src="image_18.png" /><figcaption>Figure 16. Distribution of $log_{10}(\chi^2_{\nu,b})$ for simulated RRab light curves in each band binned by magnitude. The mean $log_{10}(\chi^2_{\nu,b})$ in each band are plotted with solid curves, while the first percentile values for $log_{10}(\chi^2_{\nu,b})$ in each band are plotted with dashed curves. The black curve shows the best-fit quadratic curve to the first percentile value of the i-band curve, which was used to transform the $log_{10}(\chi^2_{\nu,b})$ features for target selection.</figcaption></figure>

uncertainties on these values added in quadrature:

$$W_Range_b = \frac{\text{mag}_{b,\text{max}} - \text{mag}_{b,\text{min}}}{[\sigma_{\text{mag}_{b,\text{max}}}^2 + \sigma_{\text{mag}_{b,\text{min}}}^2]^{1/2}} \quad (A9)$$

Fainter stars have larger photometric errors and are harder to identify as variable, as reflected in smaller values of $log(\chi^2_{\nu,b})$ (Figure 16). To avoid rejecting faint variable objects while retaining high purity for bright variables, we define a threshold on $log(\chi^2_{\nu,b})$ that changes with magnitude. To derive this threshold, we binned the simulated RRab by their weighted-average magnitudes in bins of 1 magnitude and fit a quadratic curve to the 1% percentile value of the $log_{10}(\chi^2_{\nu,i})$, shown by the black curve in Figure 16, which follows the form

$$\log_{10}(\chi^2_{\nu,b})_{\text{th}} = \log_{10}(\chi^2_{\nu,b}) - (a_1 \text{mag}_b - a_2)^2 - a_3, \quad (A10)$$

where $mag_b$ is the weighted-average PSF magnitude of the individual measurements (WAVG MAG PSF) in band b. The best-fit values of the constants were found to be $a_1 = -0.07305$, $a_2 = 17.5165$, and $a_3 = 1.6036$. The magnitude independent quantity, $log_{10}(\chi^2_{\nu,b})_{\text{th}}$, can be thought of as broadening the $log_{10}(\chi^2_{\nu,b})$ criteria to retain high efficiency for faint variable sources. RRL have positive values or values near zero across their entire magnitude range, while nonvariable objects have more negative values.

We also define a metric for discriminating RRL based on their locations in period-amplitude space. We define a distance metric for the separation between the template-fit period and amplitude of objects with the Oosterhoff I, Oosterhoff II, and Oosterhoff intermediate sequences as parameterized by Fabrizio et al. (2019) and scaled to the g-band amplitude using $A_g/A_V = 1.29$ as determined by Vivas et al. (2020b). We calculate the distance between each object and each Oosterhoff sequence from the logarithm of the period in days and g-band amplitude in magnitudes using a rectangular approximation. These distances are denoted $D(Oo-I)_j$, $D(Oo-II)_j$, and $D(Oo-int)_j$, where $j \in [0, 2]$ indicates the jth best template fit.
**Appendix B**
**Performance Curves and Top Features for Random Forest**
**Classifiers**

Tables 5–7 present the top 10 features for the initial variable selection, second-stage variability classifier, and light-curve

selection, respectively. The performance curves for these selections are shown in Figures 17–19.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Performance curves for the initial variability classifier. Left: $F_\beta$ scores over all possible cutoff score choices. The black star denotes where $F_\beta$ is maximized. Right: receiver operating characteristic curve of the RF classifier trained on cross-matched RRab and calibration stars. The black star shows the location of the preferred cutoff score.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Performance curves for the second-stage classifier.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Performance curves for the third-stage classifier after template fitting.</figcaption></figure>
<table><caption>**Table 5**
Top 10 Features for Initial Variability Selection</caption><thead><tr><th>Feature Name</th><th>Importance</th></tr></thead><tbody><tr><td>RoMS<sub>g</sub></td><td>0.1637</td></tr><tr><td>J<sub>grizY</sub></td><td>0.1604</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,i</sub>)<sub>th</sub></td><td>0.1357</td></tr><tr><td>W_Range<sub>g</sub></td><td>0.1055</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,g</sub>)<sub>th</sub></td><td>0.0842</td></tr><tr><td>RoMS<sub>r</sub></td><td>0.0815</td></tr><tr><td>J<sub>g</sub></td><td>0.0581</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,r</sub>)<sub>th</sub></td><td>0.0571</td></tr><tr><td>J<sub>r</sub></td><td>0.0544</td></tr><tr><td>σ<sup>2</sup><sub>NXS,g</sub></td><td>0.0292</td></tr></tbody></table>

<table><caption>**Table 6**
Top 10 Features for Second Variability Classifier</caption><thead><tr><th>Feature Name</th><th>Importance</th></tr></thead><tbody><tr><td>RoMS<sub>g</sub></td><td>0.1977</td></tr><tr><td>J<sub>g</sub></td><td>0.1097</td></tr><tr><td>W_Range<sub>g</sub></td><td>0.0769</td></tr><tr><td>NAPD<sub>90,g</sub></td><td>0.0562</td></tr><tr><td>χ<sup>2</sup><sub>ν,g</sub>/χ<sup>2</sup><sub>ν,r</sub></td><td>0.0485</td></tr><tr><td>RoMS<sub>r</sub></td><td>0.0358</td></tr><tr><td>NAPD<sub>90,i</sub></td><td>0.0562</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,r</sub>)<sub>th</sub></td><td>0.0340</td></tr><tr><td>RoMS<sub>z</sub></td><td>0.0312</td></tr><tr><td>log<sub>10</sub>(χ<sup>2</sup><sub>ν,g</sub>)<sub>th</sub></td><td>0.0321</td></tr></tbody></table>

<table><caption>**Table 7**
Top 10 Features for Candidate Light-curve Selection</caption><thead><tr><th>Feature Name</th><th>Importance</th></tr></thead><tbody><tr><td>rfss<sub>ν,1</sub> - rfss<sub>ν,0</sub></td><td>0.1603</td></tr><tr><td>rfss<sub>ν,0</sub></td><td></td></tr><tr><td>rfss<sub>ν,2</sub> - rfss<sub>ν,0</sub></td><td>0.1559</td></tr><tr><td>rfss<sub>ν,0</sub></td><td></td></tr><tr><td>rfss<sub>ν,1</sub> - rfss<sub>ν,0</sub></td><td>0.0660</td></tr><tr><td>rfss<sub>ν,0</sub></td><td></td></tr><tr><td>amp<sub>0</sub></td><td>0.0644</td></tr><tr><td>rfss<sub>ν,0</sub></td><td>0.0543</td></tr><tr><td>rfss<sub>ν,2</sub> - rfss<sub>ν,0</sub></td><td>0.0448</td></tr><tr><td>1/2 (RF1_score + RF2_score)</td><td>0.0282</td></tr><tr><td>D(Oo-int)<sub>0</sub></td><td>0.0267</td></tr><tr><td>RF1_score</td><td>0.0262</td></tr><tr><td>amp<sub>0</sub></td><td>0.0241</td></tr></tbody></table>

## Appendix C
### Halo Profile Fit

We fit the halo profile using a standard binned Poisson maximum-likelihood approach. However, as with many maximum-likelihood analyses, the calculation of the predicted number of counts from our model is sufficiently complex that it merits a devoted discussion.

We define our likelihood, $\mathcal{L}$, as the product over bins, $i$, of the Poisson likelihood for observing $k$ RRab candidates given a model prediction of $\lambda$ counts:

$$ \mathcal{L} = \prod_i \frac{\lambda_i^{k_i} e^{-\lambda_i}}{k_i!}. \qquad (C1) $$

<table><caption>**Table 8**
Masked Regions for Halo RRL Study</caption><thead><tr><th>Galaxy</th><th>(α, δ)<br>(deg, J2000)</th><th>Mask Radius<br>(deg)</th></tr></thead><tbody><tr><td>LMC</td><td>(80.8938, -69.7561)</td><td>27.1</td></tr><tr><td>SMC</td><td>(13.1867, -72.8286)</td><td>12.6</td></tr><tr><td>Fornax</td><td>(39.9583, -34.4997)</td><td>2.0</td></tr><tr><td>Sculptor</td><td>(15.0183, -33.7186)</td><td>1.0</td></tr></tbody></table>

<table><caption>**Table 9**
Priors on the Broken-power-law Model Parameters</caption><thead><tr><th>Parameter Name</th><th>Prior</th><th>Range</th></tr></thead><tbody><tr><td>Normalization, ρ<sub>0</sub> (kpc<sup>-3</sup>)</td><td>uniform</td><td>[0, 10<sup>6</sup>]</td></tr><tr><td>Break radius, R<sub>0</sub> (kpc)</td><td>uniform</td><td>[10, 50]</td></tr><tr><td>Inner slope, n<sub>1</sub></td><td>uniform</td><td>[-3.5, -1.0]</td></tr><tr><td>Outer slope, n<sub>2</sub></td><td>uniform</td><td>[-3.0, -6.5]</td></tr></tbody></table>

Note. The halo density profile fit assumes a fixed halo flattening of $q = 0.7$.

It is more computationally feasible to work with the logarithm of the likelihood, and thus we define

$$ \log \mathcal{L} = \sum_i k_i \log(\lambda_i) - \lambda_i - \log(k_i!), \qquad (C2) $$

where the term $-\log(k!)$ does not depend on the model parameters and can be discarded. The model-predicted number of counts in a bin, $\lambda_i$, is a function of the elliptical Galactocentric radius of the bin, $r_{e,i}$, and the model parameters $\theta$ of the broken-power-law density model described in Equation (5):

$$ \lambda_i(r_{e,i}, \theta) = \int \rho(r_{e,i}, \theta) f(D) dV \qquad (C3) $$

where $f(D)$ is the efficiency of detecting RRab at the heliocentric distance, $D$, of each volume element. This formulation naturally incorporates the geometric effects of the solar offset from the Galactic center, which are most noticeable for distances of $D \lesssim 15$ kpc (i.e., the heliocentric distance cut translates to a cut in $r_e$). Numerically, we perform the integration in Equation (C3) over HEALPix pixels of area $\sim 0.84 \text{ deg}^2$ (nside = 64), incorporating the coverage fraction of each pixel, which is estimated at scales of $\sim 166 \text{ arcsec}^2$ (nside = 16384).

With our likelihood thus defined, we seek to sample the posterior probability distribution as defined by Bayes' theorem. We explore the posterior probability distribution with Markov Chain Monte Carlo using the affine-invariant ensemble sampler emcee (Foreman-Mackey et al. 2013). We exclude the regions listed in Table 8 and assign uniform priors to each of the model parameters following the range described in Table 9. We sample the posterior using 100 walkers with 5000 samples each, discarding the first 1000 samples as burn-in. The resulting posterior distributions from the RRab and Oo-I subselection can be found in Figure 20. The best-fit parameters are assigned from the median of the posterior, and the errors are derived from the 16th and 84th percentiles of the posterior distribution.
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Posterior probability distributions for the parameters of the broken-power-law fits to the full RRab candidate sample (top) and the high-confidence sample associated with the Oosterhoff I sequence (Section 5.6).</figcaption></figure>

As can be seen in Figure 20, there is significant correlation between the normalization parameter, $\rho_0$, and the inner power-law slope, $n_1$. We note that in this regime, the analysis is especially sensitive to the geometric corrections described above for calculating the model-predicted counts (i.e., Equation (C3)).

## ORCID iDs

K. M. Stringer https://orcid.org/0000-0002-4624-2772

A. Drlica-Wagner https://orcid.org/0000-0001-8251-933X

L. Macri https://orcid.org/0000-0002-1775-4859

C. E. Martínez-Vázquez https://orcid.org/0000-0002-9144-7726

A. K. Vivas https://orcid.org/0000-0003-4341-6172

P. Ferguson https://orcid.org/0000-0001-6957-1627

A. B. Pace https://orcid.org/0000-0002-6021-8760

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

E. Neilsen https://orcid.org/0000-0002-7357-0317

K. Tavangar https://orcid.org/0000-0001-6584-6144

W. Wester https://orcid.org/0000-0003-0072-6736

S. Allam https://orcid.org/0000-0002-7069-7857

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

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

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

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

M. Costanzi https://orcid.org/0000-0001-8158-1449

M. Crocce https://orcid.org/0000-0002-9745-6228

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

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

E. Gaztanaga https://orcid.org/0000-0001-9632-0815

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

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

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

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

K. Honscheid https://orcid.org/0000-0002-6550-2023

B. Hoyle https://orcid.org/0000-0002-2571-1357

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

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

N. Kuropatkin https://orcid.org/0000-0003-2511-0946

T. S. Li https://orcid.org/0000-0002-9110-6163

M. A. G. Maia https://orcid.org/0000-0001-9856-9307

J. L. Marshall https://orcid.org/0000-0003-0710-9474

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

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

R. Morgan https://orcid.org/0000-0002-7016-5471

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

A. Palmese https://orcid.org/0000-0002-6011-0530

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

A. Roodman https://orcid.org/0000-0001-5326-3486

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

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

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

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

D. Thomas https://orcid.org/0000-0002-6325-5671

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

R. D. Wilkinson https://orcid.org/0000-0002-3908-7313

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

## References

Amit, Y., & Geman, D. 1997, Neural Comput., 9, 1545

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Baeza-Yates, R. A., & Ribeiro-Neto, B. 1999, Modern Information Retrieval (Boston, MA: Addison-Wesley)
Baker, M., & Willman, B. 2015, AJ, 150, 160

Bechtol, K., Drlica-Wagner, A., Balbinot, E., et al. 2015, ApJ, 807, 50

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Belokurov, V., & Koposov, S. E. 2016, MNRAS, 456, 602

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJL, 642, L137

Bernard, E. J., Monelli, M., Gallart, C., et al. 2010, ApJ, 712, 1259

Bersier, D., & Wood, P. R. 2002, AJ, 123, 840

Bertin, E., & Arnouts, S. 1996, A&As, 117, 393

Blažko, S. 1907, AN, 175, 325

Bonaca, A., & Hogg, D. W. 2018, ApJ, 867, 101

Bono, G., Dall'Ora, M., Caputo, F., et al. 2011, in RR Lyrae Stars, Metal-Poor Stars, and the Galaxy, ed. A. McWilliam, Vol. 5 (Pasadena, CA: Carnegie Observatories), 1

Boubert, D., Belokurov, V., Erkal, D., & Iorio, G. 2019, MNRAS, 482, 4562

Breiman, L. 2001, Mach. Learn., 45, 5

Buchler, J. R., & Kolláth, Z. 2011, ApJ, 731, 24

Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931

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

Cacciari, C., & Clementini, G. 2003, in Stellar Candles for the Extragalactic Distance Scale, Lecture Notes in Physics, ed. D. Alloin & W. Gieren, Vol. 635 (Berlin: Springer), 105

Cáceres, C., & Catelan, M. 2008, ApJS, 179, 242

Catelan, M., Pritzl, B. J., & Smith, H. A. 2004, ApJS, 154, 633

Chiti, A., Frebel, A., Jerjen, H., Kim, D., & Norris, J. E. 2020, ApJ, 891, 8

Clement, C. M., Muzzin, A., Dufton, Q., et al. 2001, AJ, 122, 2587

Clementini, G., Ripepi, V., Molinaro, R., et al. 2019, A&A, 622, A60

Crnojević, D., Sand, D. J., Zaritsky, D., et al. 2016, ApJL, 824, L14

Das, P., Williams, A., & Binney, J. 2016, MNRAS, 463, 3169

de Boer, T. J. L., Tolstoy, E., Hill, V., et al. 2012, A&A, 544, A73

de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134

DES Collaboration, Abbott, T., Abdalla, F. B., et al. 2016, MNRAS, 460, 1270

DES Collaboration, Abbott, T. M. C., Abdalla, F. B., et al. 2018, ApJS, 239, 18

Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903

Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840

Deason, A. J., Belokurov, V., Evans, N. W., & Johnston, K. V. 2013, ApJ, 763, 113

Deason, A. J., Belokurov, V., Koposov, S. E., & Lancaster, L. 2018, ApJL, 862, L1

Deason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30

DES Collaboration 2005, arXiv:astro-ph/0510346

Diehl, H. T., Neilsen, E., Gruendl, R., et al. 2016, Proc. SPIE, 9910, 99101D

Diehl, H. T., Yanny, B., Tucker, D. L., Paz-Chinchón, F., & Neilsen, E. 2019, The Dark Energy Survey and Operations: Year 6 – The Finale, OSTI Technical Report, 1596042

Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013a, ApJ, 763, 32

Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013b, ApJ, 765, 154

Drake, A. J., Djorgovski, S. G., Catelan, M., et al. 2017, MNRAS, 469, 3688

Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014, ApJS, 213, 9

Drlica-Wagner, A., Bechtol, K., Mau, S., et al. 2020, ApJ, 893, 47

Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al. 2015, ApJ, 813, 109

Efron, B. 1982, The Jackknife, the Bootstrap and other Resampling Plans (Philadelphia, PA: SIAM)

Enoch, M. L., Brown, M. E., & Burgasser, A. J. 2003, AJ, 126, 1006

Fabrizio, M., Bono, G., Braga, V. F., et al. 2019, ApJ, 882, 169

Faccioli, L., Smith, M. C., Yuan, H. B., et al. 2014, ApJ, 788, 105

Fiorentino, G., Bono, G., Monelli, M., et al. 2015, ApJL, 798, L12

Fiorentino, G., Monelli, M., Stetson, P. B., et al. 2017, A&A, 599, A125

Fitzpatrick, E. L. 1999, PASP, 111, 63

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

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

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gallart, C., Aparicio, A., Freedman, W. L., et al. 2004, AJ, 127, 1486

Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, NatAs, 3, 932

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10

Greco, C., Clementini, G., Held, E. V., et al. 2005, arXiv:astro-ph/0507244

Guldenschuh, K. A., Layden, A. C., Wan, Y., et al. 2005, PASP, 117, 721

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Helmi, A. 2008, A&Ar, 15, 145

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Hernitschek, N., Schlafly, E. F., Sesar, B., et al. 2016, ApJ, 817, 73

Holl, B., Audard, M., Nienartowicz, K., et al. 2018, A&A, 618, A30

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

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

Iorio, G., Belokurov, V., Erkal, D., et al. 2018, MNRAS, 474, 2142

Ivezić, Z., Lupton, R., Schlegel, D., et al. 2004, in ASP Conf. Ser. 317, Milky Way Surveys: The Structure and Evolution of our Galaxy, ed. D. Clemens, R. Shah, & T. Brainerd (San Francisco, CA: ASP), 179

Ivezić, Ž, Vivas, A. K., Lupton, R. H., & Zinn, R. 2005, AJ, 129, 1096

Jerzykiewicz, M., & Wenzel, W. 1977, AcA, 27, 35

Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936

Johnston, K. V., Spergel, D. N., & Haydn, C. 2002, ApJ, 570, 656

Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2012, ApJ, 761, 98

Karczmarek, P., Pietrzyński, G., Górski, M., Gieren, W., & Bersier, D. 2017, AJ, 154, 263

Koposov, S. E., Belokurov, V., Li, T. S., et al. 2019, MNRAS, 485, 4726

Kuehn, C. A., Dame, K., Smith, H. A., et al. 2013, AJ, 145, 160

Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8

MacLeod, C. L., Ivezić, Ž, Sesar, B., et al. 2012, ApJ, 753, 106

Marconi, M. 2012, MSAIS, 19, 138

Marconi, M., Coppola, G., Bono, G., et al. 2015, ApJ, 808, 50

Martínez-Vázquez, C. E., Monelli, M., Bernard, E. J., et al. 2017, ApJ, 850, 137

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

Martínez-Vázquez, C. E., Monelli, M., Gallart, C., et al. 2016b, MNRAS, 461, L41

Martínez-Vázquez, C. E., Stetson, P. B., Monelli, M., et al. 2016a, MNRAS, 462, 4349

Martínez-Vázquez, C. E., Vivas, A. K., Gurevich, M., et al. 2019, MNRAS, 490, 2183

Mateu, C., Read, J. I., & Kawata, D. 2018, MNRAS, 474, 4112

Matthews, B. W. 1975, BBA, 405, 442

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

Medina, G. E., Muñoz, R. R., Vivas, A. K., et al. 2018, ApJ, 855, 43

Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge: Cambridge Univ. Press)

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

Munoz, R. R., Cote, P., Santana, F. A., et al. 2018, ApJ, 860, 66

Nakoneczny, S., Bilicki, M., Solarz, A., et al. 2019, A&A, 624, A13

Nandra, K., George, I. M., Mushotzky, R. F., Turner, T. J., & Yaqoob, T. 1997, ApJ, 476, 70

Neilsen, E. H. J., Annis, J. T., Diehl, H. T., et al. 2019, arXiv:1912.06254

Neilsen, E. H., Jr., Bernstein, G., Gruendl, R., & Kent, S. 2016, Limiting Magnitude, τ, teff, and Image Quality in DES Year 1, Fermilab Technical Report, FERMILAB-TM-2610-AE-CD

Nidever, D. L., Olsen, K., Choi, Y., et al. 2019, ApJ, 874, 118

Oosterhoff, P. T. 1939, Obs, 62, 104

Pace, A. B., & Li, T. S. 2019, ApJ, 875, 77

Pattoero, M. 2004, Calculating efficiencies and their uncertainties, OSTI Technical Report, 15017262

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2012, arXiv:1201.0490

Pieres, A., Girardi, L., Balbinot, E., et al. 2020, MNRAS, 497, 1547

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

Pila-Díez, B., de Jong, J. T. A., Kuijken, K., van der Burg, R. F. J., & Hoekstra, H. 2015, A&A, 579, A38

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

Powers, D. 2008, J. Mach. Learn. Technol., 2, 37

Price-Whelan, A. M., Mateu, C., Iorio, G., et al. 2019, AJ, 158, 223

Prudil, Z., Grebel, E. K., Dékány, I., & Smolec, R. 2018, MNRAS, 480, 669

Ramos, P., Mateu, C., Antoja, T., et al. 2020, A&A, 638, A104

Rimoldini, L., Holl, B., Audard, M., et al. 2019, A&A, 625, A97

Rousseeuw, P. J., & Croux, C. 1993, J. Am. Stat. Assoc., 88, 1273

Rusakov, V., Monelli, M., Gallart, C., et al. 2021, MNRAS, 502, 642

Saha, A. 1985, ApJ, 289, 310

Saha, A., Olszewski, E. W., Brondel, B., et al. 2010, AJ, 140, 1719

Sanderson, R. E., Secunda, A., Johnston, K. V., & Bochanski, J. J. 2017, MNRAS, 470, S04

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

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

Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360

Sesar, B., Banholzer, S. R., Cohen, J. G., et al. 2014, ApJ, 793, 135

Sesar, B., Hernitschek, N., Mitrović, S., et al. 2011, AJ, 153, 204

Sesar, B., Ivezić, Ž, Grammer, S. H., et al. 2010, ApJ, 708, 717

Sesar, B., Jurić, M., & Ivezić, Ž. 2011, ApJ, 731, 4

Sevilla-Noarbe, I., Bechtol, K., & Kind, M. C. 2020, arXiv:2011.03407
Sharma, S., Bland-Hawthorn, J., Johnston, K. V., & Binney, J. 2011a, ApJ, 730, 3

Sharma, S., Johnston, K. V., Majewski, S. R., Bullock, J., & Muñoz, R. R. 2011b, ApJ, 728, 106

Sheffield, A. A., Johnston, K. V., Majewski, S. R., et al. 2014, ApJ, 793, 62

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

Simm, T., Saglia, R., Salvato, M., et al. 2015, A&A, 584, A106

Simon, J. D., Li, T. S., Erkal, D., et al. 2020, ApJ, 892, 137

Smith, H. A. 1995, RR Lyrae Stars, Cambridge Astrophysics Series, Vol. 27 (Cambridge: Cambridge Univ. Press), 27

Sokolovsky, K. V., Gavras, P., Karampelas, A., et al. 2017, MNRAS, 464, 274

Soszyński, I., Udalski, A., Szymański, M. K., et al. 2017, AcA, 67, 103

Soszyński, I., Udalski, A., Szymański, M. K., et al. 2019, AcA, 69, 87

Stetson, P. B. 1996, PASP, 108, 851

Stetson, P. B., Fiorentino, G., Bono, G., et al. 2014, PASP, 126, 616

Stringer, K. M., Long, J. P., Macri, L. M., et al. 2019, AJ, 158, 16

Tie, S. S., Martini, P., Mudd, D., et al. 2017, AJ, 153, 107

Tolstoy, E., Hill, V., & Tosi, M. 2009, ARA&A, 47, 371

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

Torrealba, G., Catelan, M., Drake, A. J., et al. 2015, MNRAS, 446, 2251

Udalski, A., Wyrzykowski, L., Pietrzyński, G., et al. 2001, AcA, 51, 221

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vivas, A. K., Martínez-Vázquez, C., & Walker, A. R. 2020a, ApJS, 247, 35

Vivas, A. K., Saha, A., Olsen, K., et al. 2017, AJ, 154, 85

Vivas, A. K., Walker, A. R., Martínez-Vázquez, C. E., et al. 2020b, MNRAS, 492, 1061

Vivas, A. K., & Zinn, R. 2006, AJ, 132, 714

Vivas, A. K., Zinn, R., Abad, C., et al. 2004, AJ, 127, 1158

Vivas, A. K., Zinn, R., Andrews, P., et al. 2001, ApJL, 554, L33

Walker, A. R. 1989, PASP, 101, 570

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

Watkins, L. L., Evans, N. W., Belokurov, V., et al. 2009, MNRAS, 398, 1757

Welch, D. L., & Stetson, P. B. 1993, AJ, 105, 1813

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

Xue, X.-X., Rix, H.-W., Ma, Z., et al. 2015, ApJ, 809, 144

Zinn, R., Horowitz, B., Vivas, A. K., et al. 2014, ApJ, 781, 22