Identification of RR Lyrae stars in multiband,
sparsely-sampled data from the Dark Energy Survey
using template fitting and Random Forest classification

Katelyn M. Stringer, J.P. Long, L.M. Macri, J.L. Marshall, A. Drlica-Wagner, C.E. Martínez-Vázquez, A.K. Vivas, K. Bechtol, E. Morganson, M. Carrasco Kind, et al.

► To cite this version:

Katelyn M. Stringer, J.P. Long, L.M. Macri, J.L. Marshall, A. Drlica-Wagner, et al. Identification of RR Lyrae stars in multiband, sparsely-sampled data from the Dark Energy Survey using template fitting and Random Forest classification. Astron.J., 2019, 158 (1), pp.16. 10.3847/1538-3881/ab1f46 . hal-02136317

HAL Id: hal-02136317

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

Submitted on 29 Sep 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
<figure><img src="image_2.png" /></figure>

# Identification of RR Lyrae Stars in Multiband, Sparsely Sampled Data from the Dark Energy Survey Using Template Fitting and Random Forest Classification

K. M. Stringer¹, J. P. Long², L. M. Macri¹, J. L. Marshall¹, A. Drlica-Wagner³,⁴, C. E. Martínez-Vázquez⁵, A. K. Vivas⁵, K. Bechtol⁶,⁷, E. Morganson⁸,⁹, M. Carrasco Kind⁸,⁹, A. B. Pace¹, A. R. Walker⁵, C. Nielsen¹,¹⁰, T. S. Li¹¹, E. Rykoff¹²,¹³, D. Burke¹²,¹³, A. Carnero Rosell¹⁴,¹⁵, E. Neilsen¹¹, P. Ferguson¹, S. A. Cantu¹, J. L. Myron¹, L. Strigari¹, A. Farahi¹⁶, F. Paz-Chinchón⁸,⁹, D. Tucker³, Z. Lin², D. Hatt³, J. F. Maner¹, L. Plybon¹, A. H. Riley¹, E. O. Nadler¹², T. M. C. Abbott⁵, S. Allam¹¹, J. Annis¹¹, E. Bertin¹⁷,¹⁸, D. Brooks¹⁹, E. Buckley-Geer¹¹, J. Carretero²⁰, C. E. Cunha¹², C. B. D'Andrea²¹, L. N. da Costa¹⁴,¹⁵, J. De Vicente²², S. Desai²³, P. Doel¹⁹, T. F. Eifler²⁴,²⁵, B. Flaugher¹¹, J. Frieman³,¹¹, J. García-Bellido²⁶, E. Gaztanaga²⁷,²⁸, D. Gruen¹²,¹³,²⁹, J. Gschwend¹⁴,¹⁵, G. Gutierrez¹¹, W. G. Hartley¹⁹,³⁰, D. L. Hollowood³¹, B. Hoyle³²,³³, D. J. James³⁴, K. Kuehn³⁵, N. Kuropatkin¹¹, P. Melchior³⁶, R. Miquel²⁰,³⁷, D. L. C. Ogando¹⁴,¹⁵, A. A. Plazas²⁵, E. Sanchez²², B. Santiago¹⁴,³⁸, V. Scarpine¹¹, M. Schubnell³⁹, S. Serrano²⁷,²⁸, I. Sevilla-Noarbe²², M. Smith⁴⁰, R. C. Smith⁵, M. Soares-Santos⁴¹, F. Sobreira¹⁴,⁴², E. Suchyta⁴³, M. E. C. Swanson⁸, G. Tarle³⁹, D. Thomas⁴⁴, V. Vikram⁴⁵, and B. Yanny¹¹

(DES Collaboration)

¹ George P. and Cynthia W. Mitchell Institute for Fundamental Physics and Astronomy, Department of Physics and Astronomy, Texas A & M University, College Station, TX 77843, USA; kstringer@tamu.edu

² Department of Statistics, Texas A&M University, 3143 TAMU, College Station, TX 77843, USA

³ 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, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

⁶ Wisconsin IceCube Particle Astrophysics Center (WIPAC), Madison, WI 53703, USA

⁷ Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA

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

⁹ Department of Astronomy, University of Illinois, Urbana, IL 61801, USA

¹⁰ Department of Physics, Purdue University, 1396 Physics Building, West Lafayette, Indiana 47907-1396, USA

¹¹ Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA

¹² Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA

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

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

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

¹⁶ Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15312, 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

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

²¹ Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

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

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

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

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

²⁶ Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain

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

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

³⁰ Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

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

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

³⁴ Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA

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

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

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

³⁸ Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS—91501-970, Brazil

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

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

⁴¹ Physics Department, Brandeis University, 415 South Street, Waltham, MA 02453, USA

⁴² Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil

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

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

⁴⁵ Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA

Received 2018 December 29; revised 2019 April 29; accepted 2019 April 30; published 2019 June 14

## Abstract

Many studies have shown that RR Lyrae variable stars (RRL) are powerful stellar tracers of Galactic halo structure and satellite galaxies. The Dark Energy Survey (DES), with its deep and wide coverage ($g \sim 23.5$ mag in a single
exposure; over 5000 deg²) provides a rich opportunity to search for substructures out to the edge of the Milky Way halo. However, the sparse and unevenly sampled multiband light curves from the DES wide-field survey (a median of four observations in each of *grizY* over the first three years) pose a challenge for traditional techniques used to detect RRL. We present an empirically motivated and computationally efficient template-fitting method to identify these variable stars using three years of DES data. When tested on DES light curves of previously classified objects in SDSS stripe 82, our algorithm recovers 89% of RRL periods to within 1% of their true value with 85% purity and 76% completeness. Using this method, we identify 5783 RRL candidates, ~28% of which are previously undiscovered. This method will be useful for identifying RRL in other sparse multiband data sets.

**Key words:** catalogs – galaxy: halo – galaxy: structure – methods: statistical – stars: variables: RR Lyrae

*Supporting material:* machine-readable tables

# 1. Introduction

RR Lyrae variable stars (RRL) are old (age >10 Gyr) horizontal branch stars that pulsate with short periods (0.2–1.2 days). They have become one of the most widely used stellar tracers in Milky Way and Local Group studies. Thanks to the discovery of RR Lyrae itself (Pickering et al. 1901) and the subsequent studies of their pulsation (see King & Cox 1968 and Catelan & Smith 2015 for a review of the pioneering and current works in this field), these stars have well-understood period–luminosity–metallicity (P–L–Z) relations<fnref n="46" /> (e.g., Cáceres & Catelan 2008; Marconi et al. 2015), making them excellent distance indicators, especially in the near-infrared bands. This, combined with their bright luminosities ($M_V \sim 0.6$) and advanced ages make RRL well suited to trace discrete stellar populations (satellite galaxies, star clusters, and streams) within the Milky Way halo (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).

Locating these stellar populations is crucial for testing the $\Lambda$CDM hierarchical model, which predicts that the halos of large galaxies like the Milky Way are formed through the accretion and disruption of lower mass halos (Bullock & Johnston 2005). Recent re-examinations of these simulations predict that the outer reaches of the stellar halo ($d \gtrsim 100$ kpc) are primarily composed of the most recently accreted satellites and that thousands of RRL should be present in them (Sanderson et al. 2017). Once satellite galaxies and their disrupted remains are found, their distribution and properties can reveal valuable clues about the formation history, dark matter density profile, and mass of the Milky Way. While these objects are interesting in their own right, the statistical information about this sample is vital in placing the Milky Way in a broader cosmological context.

Numerous Milky Way substructures have already been discovered. Eleven “classical” dwarf galaxies were known to orbit the Milky Way before 2005 (McConnachie 2012).<fnref n="47" /> Thanks to the advent of wide-field surveys such as the Sloan Digital Sky Survey (SDSS; York et al. 2000), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016), and the Dark Energy Survey (DES; Dark Energy Survey Collaboration et al. 2016), over 40 new dwarf satellite candidates have been discovered (Willman et al. 2005a, 2005b; Belokurov et al. 2006b, 2007, 2008, 2009, 2010; Grillmair 2006, 2009; Sakamoto & Hasegawa 2006; Zucker

et al. 2006a, 2006b; Irwin et al. 2007; Walsh et al. 2007; Bechtol et al. 2015; Drlica-Wagner et al. 2015; Kim et al. 2015; Koposov et al. 2015, 2018; Laevens et al. 2015a, 2015b; Luque et al. 2016, 2017; Torrealba et al. 2016a, 2016b, 2018). In addition to galaxies that are still intact, tidal streams, the disrupted remains of satellite galaxies and globular clusters, have been discovered to be prevalent within the Milky Way halo (e.g., Odenkirchen et al. 2001; Newberg et al. 2002; Belokurov et al. 2006a; Bernard et al. 2014, 2016; Koposov et al. 2014; Balbinot et al. 2016; Grillmair & Carlin 2016; Mateu et al. 2018; Shipp et al. 2018). Aside from these, additional large stellar overdensities populate the Milky Way stellar halo with origins still unknown (e.g., Vivas et al. 2001; Newberg et al. 2002; Majewski et al. 2004; Rocha-Pinto et al. 2004; Belokurov et al. 2007; Sesar et al. 2007; Sharma et al. 2010; Deason et al. 2014; Li et al. 2016; Pieres et al. 2017; Bergemann et al. 2018; Prudil et al. 2018).

Most of these satellites and streams were discovered as stellar overdensities in photometric catalogs (Willman 2010 and references therein). However, this detection method is biased against diffuse objects with low surface brightness ($\mu_{V,0} \gtrsim 29$ mag/arcsec²; Baker & Willman 2015), so an alternative method is needed to locate other faint structures that may have evaded detection. RRL are sufficiently rare so as to not randomly form in pairs outside of stellar structures, so searching for groups of spatially close RRL provides an independent method to detect new structures (Ivezić et al. 2004a; Sesar et al. 2014; Baker & Willman 2015; Medina et al. 2017, 2018). Indeed, at least one RRL has been found in almost every satellite galaxy with available time series data<fnref n="48" /> (Boettcher et al. 2013; Vivas et al. 2016; Martínez-Vázquez et al. 2017, and references therein). Thus, identifying RRL in the halo can increase the census of old, metal-poor satellite galaxies, streams, and overdensities, and improve our understanding of the Milky Way.

The two most common subtypes of RRL are those pulsating in the fundamental mode, RRab, and those pulsating in the first overtone, RRc. When their light curves are adequately sampled, RRab are easily identified by their short periods (0.4 $\lesssim$ P $\lesssim$ 1 day), relatively large pulsation amplitudes (0.5 $\lesssim$ A_g $\lesssim$ 1.5 mag), and a characteristic sawtooth shape. RRc have shorter periods (0.2 $\lesssim$ P $\lesssim$ 0.45 days), smaller amplitudes (0.2 $\lesssim$ A_g $\lesssim$ 0.8 mag), more sinusoidal-shaped light curves, and are generally less numerous than RRab. The fraction of RRab to RRc and the average periods of each are highly dependent on the metallicity of the stellar population in which they formed and is still not fully understood (see Catelan 2009 and

<footnote n="46">⁴⁶ These are sometimes presented as period–luminosity–color (P–L–C) relations.</footnote>

<footnote n="47">⁴⁷ The nature of the Canis Major Overdensity as a satellite galaxy is in doubt, due to a lack of an RRL excess and a potential warp in the Milky Way disk (Mateu et al. 2009).</footnote>

<footnote n="48">⁴⁸ One notable exception is the satellite galaxy candidate Carina III, which currently has no detected RRL in its vicinity (Torrealba et al. 2018).</footnote>
references therein). Most populations of RRL in the Milky Way are commonly subdivided into Oosterhoff I, II, and III groups based on these observational properties (named after the first dichotomy applied to globular clusters by Oosterhoff 1939). We refer the interested reader to Table 6 in Martínez-Vázquez et al. (2017) for a summary of these properties for a selection of Local Group dwarf galaxies.

Period-finding algorithms have long been used in conjunction with visual inspection to identify RRL from their time series photometry. However, with the dramatic increases in available data in recent years, the need for automated detection algorithms has grown significantly. Stetson (1996) made great strides in this regard when he introduced an automated method to identify Cepheid variables using template light curves to estimate their periods and a scoring system based on calculated variability indices. Recent studies have extended period-finding techniques to multiple filters (Mateu et al. 2012; Mondrik et al. 2015; VanderPlas & Ivezić 2015; Saha & Vivas 2017; Huijse et al. 2018). However, even these algorithms suffer in performance when applied to extremely sparsely sampled data. Hernitschek et al. (2016) and Sesar et al. (2017) developed separate techniques to identify RRL in the sparsely sampled multiband Pan-STARRS data (Chambers et al. 2016) and found thousands of such variables.

We add to this census by presenting new RRL candidates discovered in the first three years of the DES data. DES is a five-year multiband (*grizY*) imaging survey using the Dark Energy Camera (Flaugher et al. 2015) on the 4m Blanco Telescope at Cerro Tololo Inter-American Observatory (CTIO). After the conclusion of its observations, DES will provide a deep (~25 mag in the final coadded images) and wide (~5000 deg²) data set near the Southern Galactic cap (The Dark Energy Survey Collaboration 2005; Diehl et al. 2016). By the end of the survey, the entire footprint will have been imaged ~10 times in each band. While the main goal of the survey is to better constrain certain cosmological parameters, the deep and wide survey data provide an excellent test bed for probing Milky Way substructure with RRL. However, like Pan-STARRS, the DES light curves are multiband and poorly sampled. In this paper, we detail how we overcome these challenges by creating an empirically derived light curve template and a computationally efficient fitting algorithm to determine periods and other light curve parameters. We use these methods to identify 5783 RRL candidates, 28% of which are new discoveries, including three with a heliocentric distance >220 kpc. This novel technique will prove useful for other sparsely sampled multiband data sets.

The rest of the paper is organized as follows: Section 2 describes how we extracted star-like objects from the DES Y3 data release; Section 3 explains how we rescaled the photometric uncertainties and applied simple metrics to select variable objects; Section 4 presents the multiband RRL template, its application to DES light curves, and the construction and performance of our random forest classifier; Section 5 presents our RRL catalog, a comparison to overlapping surveys, and parameter uncertainties; and Section 6 discusses possible biases, the spatial distribution of the candidates, and potential future application for the Large Synoptic Survey Telescope (LSST).

## 2. Data

### 2.1. DES Year 3 Quick Release

This work is based on the DES internal Year 3 Quick Release catalog (hereafter Y3Q2), which contains all the single-epoch data from years 1–3 that formed the basis for the coadded DES first public data release<fnref n="49" /> (Abbott et al. 2018, hereafter DR1). The Y3Q2 data set was developed in the same manner as the Y2Q1 data release used for the stellar overdensity searches in Drlica-Wagner et al. (2015) (see their Section 2.1 for a detailed description of how the Quick Release catalogs were generated), with one major change. Instead of using stellar locus regression (SLR; Ivezić et al. 2004b; MacDonald et al. 2004; High et al. 2009; Gilbank et al. 2011; Coupon et al. 2012; Desai et al. 2012; Kelly et al. 2014) to determine zero points for the absolute photometric calibration, the Y3Q2 release utilizes the Forward Global Calibration Module (FGCM) photometric zero points (Burke et al. 2018). Y3Q2 contains single-epoch catalogs generated from the reduced FINALCUT DES images (Morganson et al. 2018) and a cross-matched “coadded catalog” generated from these single-epoch measurements. This catalog does not contain information from exposures in which an object was not detected at approximately a 5σ level. The DES Y3Q2 coadded catalog contains nearly 2.9 × 10⁸ unique objects and spans the entire survey footprint with S/N ~ 10 at a median depth of 23.5, 23.3, 22.8, 22.1, and 20.7 mag in *grizY*, respectively (Abbott et al. 2018).

As DES images are collected, the filter to be used and the location to be imaged are prioritized according to the time of the year, the sky conditions (Moon phase, seeing, weather), and how many times that particular area has already been imaged (Neilsen & Annis 2014; see Figure 3 in Diehl et al. 2016). While this strategy ensures uniform depth and the best use of the observing time, objects in the wide-field survey are sampled with a highly unpredictable cadence. In the Y3Q2 data set, individual objects can have from two to over 50 observations depending on their location. We ensure that each light curve only contains photometric observations by requiring that each observation has a SExtractor warning value FLAG ≤ 4, is sufficiently far away from masked regions in the images (IMAFLAG 4), and has a zero-point correction available (FGCM_FLAG ≤ 4). After these cuts, the median number of total observations for a given object is 10, while the median number of observations in each band across the survey region is 4. The effects of the survey coverage and these cuts are discussed more in Section 6.1.

### 2.2. Object Selection

We selected our objects using the coadded catalog before examining the time series data, because the former contained most of the information needed to identify candidates (such as the number of times each object was imaged in each band and the star–galaxy classification). We further restricted the sample to stellar-like objects by following a prescription similar to Bechtol et al. (2015), based on the SExtractor (Bertin & Arnouts 1996) SPREAD MODEL parameter, which selects stars well down to r ~ 23 (Drlica-Wagner et al. 2015). Lastly, we required at least five total observations to be able to search for variability.

<footnote n="49">⁴⁹ https://des.ncsa.illinois.edu/releases/dr1</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: variation of $log_{10}(\chi^2_{\nu,r})$ vs. the median $r$ magnitude, $\bar{m}_r$, demonstrating that photometric errors are slightly overestimated for brighter objects in the DES pipeline. Red points were excluded using an iterative $3\sigma$ clipping procedure. The black curve shows the quadratic fit that was used to rescale the errors. Right: distribution after the photometric errors were rescaled.</figcaption></figure>

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

We selected objects that are bright enough to be detected in multiple images by requiring the coadded PSF (WAVG\_MAG\_PSF) or the aperture magnitudes in the exposure with the best seeing in that band (MAG\_AUTO) to be brighter than the median depth of the Y3Q2 single-epoch exposures across the entire survey region (see Section 2.1). We excluded all objects for which the coadded photometry errors exceed 0.3 mag in each of griz to reduce the number of spurious detections.

To ensure that we did not discard stars with missing data in a single band, we considered these quantities separately for each of griz. We did not use the Y data for these initial selections because those exposures are generally taken under worse seeing conditions than the other bands (Diehl et al. 2016) and are thus a poor choice to use for star–galaxy separation. The star–galaxy separation we used performs best in riz due to the better seeing conditions for those observations, as discussed in the previous section. Including objects which passed this cut in g likely allowed some extended sources into our sample, which we discuss further in Section 5.1.

Although RRL have well-characterized colors, we did not employ a color cut in this early stage of the analysis because we did not want to exclude any potential RRL with poor coverage across filters or pulsation phase. Simultaneous colors were not available for some objects in the DES footprint, so we would have to calculate colors using coadded magnitudes or magnitudes from arbitrary phases in the star’s variation to calculate colors, which would expand the range of possible colors for RRL in our data. The RRL template we describe in Section 4.1 provides the color information we need to identify RRL candidates. In the future, when more epochs of DES data are available, color cuts will be a more reliable RRL indicator prior to the template fitting.

In summary, our combined selection criteria were:

1. $\ge 2$ observations in g, r, i, or z;

2. $|SPREAD\_MODEL| \le (0.003 + SPREADERR\_MODEL)$ in g, r, i, or z;

3. $(16 \le \text{WAVG\_MAG\_PSF} \le \text{median depth})$ or $(16 \le \text{MAG\_AUTO} \le \text{median depth})$ in g, r, i, or z;

4. $(\text{WAVG\_MAGERR\_PSF} < 0.3)$ or $(\text{MAGERR\_AUTO} < 0.3)$ in g, r, i, or z.

A sample of $\sim 1.5 \times 10^8$ objects passed all of these combined selection criteria. We used their time series data instead of their coadded values for the remainder of our analysis.

## 3. Variability Analysis

### 3.1. Error Rescaling

Photometric uncertainties can have a large impact on the success of our variability classification. We first account for both the photometric uncertainties reported by the DES pipeline and the uncertainties in the FGCM zero-point solution<fnref n="50" /> for each exposure by adding them in quadrature. Because photometric uncertainties can be over- or underestimated for different magnitude ranges (Kaluzny et al. 1998), we calculated the reduced chi-squared statistic, $\chi^2_{\nu,b}$, from the median magnitude $\bar{m}_b$ in a given band $b$ for each light curve:

$$ \chi^2_{\nu,b} = \frac{1}{N_b - 1} \sum_{i=1}^{N_b} \frac{(m_{i,b} - \bar{m}_b)^2}{\sigma^2_{i,b}}, \quad (1) $$

where $N_b$ is the number of observations for a unique object in band $b$, $m_{i,b}$ is the $i$th observation in that band, and $\sigma_{i,b}$ is the photometric uncertainty combined with the zero-point uncertainty for that observation. As this statistic measures the goodness of fit to a constant value of $\bar{m}_b$, one would expect $\chi^2_{\nu,b} \approx 1$ for a nonvariable source and $\chi^2_{\nu,b} > 1$ for variable sources.

Because the majority of objects within a given field will have constant (nonvarying) light curves, any overall trend of $\chi^2_{\nu,b}$ versus magnitude will be indicative of incorrect estimations of photometric uncertainty. For ease of calculations, we subdivided the single-epoch data by HEALPix (nside = 32) (Górski et al. 2005)<fnref n="51" /> and filter. This resulted in 1772 unique DES HEALPix regions. We fit a quadratic function:

$$ \log_{10}(\chi^2_{\nu,b}) = c_{0,b} + c_{1,b}(\bar{m}_b - 20) + c_{2,b}(\bar{m}_b - 20)^2 \quad (2) $$

for each filter $b$ in each of these regions, excluding variables and outliers by applying an iterative $3\sigma$ clip from the median value using the `sigmae clip` function in `astropy.stats`. We show the initial trend in $\chi^2_{\nu,r}$ for the objects in HEALPix 11678 in the $r$ band in the left portion of Figure 1. We multiplied the reported photometric uncertainties of each object by scaling factors based on the best-fit value of $\chi^2_{\nu,b}$ for each of

<footnote n="50">⁵⁰ At the time of this analysis, only a prerelease version of the FGCM zero points was available. A later version of these zero points was used for other DES Year 3 analyses.</footnote>

<footnote n="51">⁵¹ These HEALPix indices are also provided in the DES DR1 products (Abbott et al. 2018).</footnote>
<table><caption>**Table 1**
Error Rescaling Coefficients for Equation (2)</caption><thead><tr><th>HEALPix</th><th>Band <i>b</i></th><th><i>c</i><sub>0,b</sub></th><th><i>c</i><sub>1,b</sub></th><th><i>c</i><sub>2,b</sub></th></tr></thead><tbody><tr><td rowspan="5">11678</td><td><i>g</i></td><td>-0.2672</td><td>0.0816</td><td>-0.0152</td></tr><tr><td><i>r</i></td><td>-0.2572</td><td>0.0793</td><td>-0.0134</td></tr><tr><td><i>i</i></td><td>-0.1180</td><td>0.0395</td><td>-0.0156</td></tr><tr><td><i>z</i></td><td>-0.2160</td><td>0.0827</td><td>-0.0198</td></tr><tr><td><i>Y</i></td><td>-0.1942</td><td>0.0505</td><td>-0.0250</td></tr></tbody></table>

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

its magnitudes in a given band. Once the uncertainties were suitably rescaled, we repeated the calculation to verify that no trends remained. The right portion of Figure 1 shows the resulting lack of trend in $\chi^2_{\nu,r}$ after the rescaling procedure. Table 1 lists the coefficients used to rescale the errors in each band for this example region (a full version is available online). Figure 2 shows the final rescaled photometric uncertainties for the entire survey region as a function of magnitude for the r band.

## 3.2. Variability Cuts

Once photometric uncertainties were rescaled, we assessed the variability of each light curve using two simple metrics. The first was $\chi^2_{\nu,b}$ described in Section 3.1. The second was a metric we called “significance,” consisting of a weighted range of the magnitudes in one band that acts as a proxy for light-curve amplitude:<fnref n="52" />

$$ \text{significance}_b = \frac{(m_{\text{max},b} - m_{\text{min},b})}{\sqrt{\sigma_{\text{max},b}^2 + \sigma_{\text{min},b}^2}} \quad (3) $$

To test the effectiveness of these metrics and determine the threshold values to separate variables from constant stars, we assembled a labeled training set of previously classified objects in the SDSS stripe 82 region (hereafter, “S82”). S82 is a ~300 deg² area spanning 300° ≲ α ≲ 60° and |δ| ≲ 1.25° that was observed 70–90 times by SDSS in ugriz over a period of 10 yr. Numerous authors used the resultant well-sampled multiband light curves to identify thousands of variable stars in the region with high confidence (Ivezić et al. 2007; Sesar et al. 2010; Süveges et al. 2012). These labeled objects are extremely useful for studies of variables from both hemispheres thanks to their equatorial location. Although the magnitude range of DES is deeper than that of SDSS, there is sufficient overlap to create a well-populated training set for our study. Using the calibration and variable catalogs from Ivezić et al. (2007), we cross-matched 641,710 “standard” (i.e., constant) stars and 16,752 variables in common between SDSS and DES objects. We also identified 296 RRL in common between DES and either Sesar et al. (2010) or Süveges et al. (2012), consisting of 238 and 58 objects of subtype RRab and RRc, respectively.

As an example, we show the cumulative distributions of the $\chi^2_{\nu,b}$ values and “significance” values for the cross-matched objects in the r band in Figure 3. For both metrics, the threshold

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Uncertainties in the r-band magnitudes for the entire survey region after the error rescaling described in Section 3.1. The dashed black line denotes the median survey depth in this band, and the solid black line shows the third-degree polynomial fit to the uncertainties.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Initial variability metric values for previously identified objects in S82. Top: cumulative distribution of $\log_{10}(\chi^2_{\nu,r})$. The magenta dashed line denotes the chosen threshold value of 0.5; objects with larger values in any band are considered variable. Bottom: histogram of $\log_{10}(\text{significance}_r)$. The vertical magenta dashed line denotes the chosen threshold value of 1.0; objects with larger values in any band are considered variable. Note that although some real RRab are excluded by this cut, most of the non-RRab variables are also excluded.</figcaption></figure>

values were chosen to minimize the number of nonvariable stars that would be subject to subsequent analysis. Any objects that showed $\log_{10}(\chi^2_{\nu,b}) \ge 0.5$ and significance $\ge 1$ in any one of grizY<fnref n="53" /> were kept for subsequent analysis. When these cuts were applied across all five filter bands, 234 (~98%) RRab, 57 (~98%) RRc, 5196 (~31%) variable, and 3004 (≤0.05%) standard light curves from S82 met these criteria. These results for our training set are summarized in Table 2. Over the entire survey region, approximately ~7 × 10⁵ light curves passed these variability cuts. We caution that passing this criterion is

<footnote n="52">⁵² While other metrics such as the Welch-Stetson *I* (Welch & Stetson 1993) and Stetson *J* (Stetson 1996) indices are widely used and very effective at detecting variability, due to the sparsity of our observations, we chose to use metrics that were agnostic of observation time.</footnote>

<footnote n="53">⁵³ Unlike the initial cuts we applied to the coadded catalog, we included the *Y* band in these variability cuts because the *Y*-band values were weighted by the photometric uncertainties.</footnote>
<table><caption>**Table 2**
Training Set of Cross-matched S82 Objects</caption><thead><tr><td>SDSS Label</td><td>Present in DES</td><td>Passed Cuts</td></tr></thead><tbody><tr><td>RRab</td><td>238</td><td>234</td></tr><tr><td>RRc</td><td>58</td><td>57</td></tr><tr><td>Variables</td><td>16,752</td><td>5196</td></tr><tr><td>Standards</td><td>641,710</td><td>3004</td></tr></tbody></table>

Note. Objects were originally identified in Ivezić et al. (2007), Sesar et al. (2010), and Süveges et al. (2012).

simply an initial cut and does not imply that these sources are truly variable.

Despite these cuts, a small but non-negligible fraction of the objects identified as “standard” by Ivezić et al. (2007) were still selected. It is possible that some of these objects are truly variable sources that did not display significant changes in the previous studies. Another possibility is erroneous photometry that, while rare, occurs sometimes in the Y3Q2 data set, due to incorrectly attributing observations from separate sources to one object or imperfect masking of observations obtained in very poor weather conditions. While these objects may have passed the initial variability cuts, their light curves fit the RRL template poorly, and most of them were rejected later in our analysis.

All of the selected light curves were corrected for extinction using reddening values from the maps of Schlegel et al. (1998) multiplied by filter coefficients derived from the Fitzpatrick (1999) reddening law (for $R_V = 3.1$) and the adjustments to the Schlegel et al. (1998) map presented by Schlafly & Finkbeiner (2011; see Section 4.2 in Abbott et al. 2018 for more details). We then used the extinction-corrected light curves as the input for our template-fitting algorithm.

# 4. Candidate Identification

## 4.1. RR Lyrae Template

Our current work introduces a novel method of identifying RRab by fitting an empirically derived periodic model to the sparsely sampled multiband light curves. The model has the form

$$m_b(t) = \mu + M_b(\omega) + a\gamma_b(\omega t + \phi), \quad (4)$$

where $m_b$ is the measured magnitude in a given band b at a given time t, $\mu$ is the distance modulus, $M_b$ is the average absolute magnitude in that band, $\omega = 1/P$ is the frequency of the variability (inverse of the period P), a is the g-band amplitude (the amplitudes of the curves for the other bands are proportional to a), $\gamma_b$ is a periodic shape function, and $\phi$ is the phase. Only the four parameters $\mu$, $a$, $\omega$, and $\phi$ are estimated during the fitting process while the forms of $\gamma_b$ and $M_b(\omega)$ are fixed. These were derived using well-sampled RRab light curves from S82 (Sesar et al. 2010) and shifted to adjust for slight differences between the SDSS and DES filters. A more detailed description of the template construction is included in Appendix A.

Using the same reasoning as Sesar et al. (2017), we chose to exclude RRc from our study because (a) their sinusoidal light curves are difficult to distinguish from light curves from other variable objects such as eclipsing binaries, (b) their small amplitudes would make them difficult to identify in our sparse data, and (c) searching over a larger period range to recover their short periods (~0.3 days) would introduce additional common period aliases into our sample. While excluding RRc weakens

<figure><img src="image_7.png" /><figcaption>**Figure 4.** P–L relations used in our template-fitting procedure. RRL have nearly constant absolute magnitudes in g regardless of period. See Section 4.1 for more details.</figcaption></figure>

our sample size for tracing substructure, it is not prohibitive because RRc are usually less numerous than RRab. Furthermore, this approach allowed us to use only one generalized RRab shape for our template instead of an ensemble of shapes as in Sesar et al. (2017). While they were able to recover a more diverse group of RRL by fitting multiple shapes, our approach makes our algorithm more computationally efficient.

The P–L–Z relations were implemented in our template-fitting procedure as P–L relations evaluated at a starting metallicity of [Fe/H] = -2.85. However, the values of these offsets between template filter curves were later adjusted to fit the S82 RRL light curves, so this metallicity should not be treated as the true value of the template. Hence, the metallicity value of our template is somewhat ambiguous. We discuss the effects of the metallicity further in Section 5.4. The final template P–L relations are shown in Figure 4. More details on how these parameters were derived are presented in Appendix B.

## 4.2. Template Fitting

In this section, we include a high-level overview of the template-fitting procedure. A more detailed description of this process can be found in Appendix C.

To fit the template to a light curve, the algorithm performs a grid search over a specified range of frequencies. To prevent misestimated and aliased periods outside of the true range of periods known for RRab, we restricted this range to values corresponding to periods of 0.44–0.89 days following the period–amplitude relation shown in Figure 16 of Sesar et al. (2010). At each frequency grid point, the algorithm first calculates $M_b(\omega)$ (see Equation (4) and Figure 4) and subtracts these values from the light-curve magnitudes in the appropriate band. Then, with the frequency fixed, the model alternates between estimating the best-fitting $\phi$ using Newton’s method and $(a, \mu)$ using weighted least squares. The $(\phi, a, \mu)$ values that minimize the weighted residual sum of squares (RSS) at each frequency grid point are chosen as the best-fitting parameters and act as the starting point for the parameter search at the next grid point. The $(\omega, \phi, a, \mu)$ at the global minimum RSS value over the entire frequency grid are chosen as the best-fitting parameters.<fnref n="54" />

<footnote n="54">⁵⁴ In his study of Cepheid variables, Stetson (1996) also developed a template-fitting method based on least squares. However, instead of using a string-length minimization technique in a single band, we use all bands simultaneously and used a fixed shape instead of a family of derived template curves calculated for each trial period.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Residual sum of squares (RSS) curve for an RRab originally discovered by Sesar et al. (2010). The red arrow denotes the global minimum of the RSS, which corresponds to the true period of 0.5336 days.</figcaption></figure>

A major strength of this algorithm is that the template shape is fit simultaneously to the light-curve data in all five bands combined. Because there are only four free parameters that must be fit for the entire light curve, unique solutions can be found for sparse light curves with very few measurements in any single band. An example RSS curve for an RRab from the labeled training set is shown in Figure 5 and its corresponding template fit is shown in Figure 6. Because the local minima in the RSS curve are sometimes very close in value, we include the best-fitting values for the top three minima of RSS in our data products for completeness, but do not discuss the results of the second and third minima further in this work.

Our algorithm is also effective at estimating distances. At the best-fit parameter estimates, the template-fitting algorithm correctly estimated ~81% of the S82 RRL distances to within 3% of the values obtained by Sesar et al. (2010) and Süveges et al. (2012; if available) with an overall standard deviation of 2.8% (see Section 5.4 for an extended discussion of the uncertainties in distance modulus). The accuracy of the template estimates of both the period and the distance for the training set of RRab light curves is summarized in Table 3.

Our algorithm is computationally efficient and only takes ~3–5 minutes per light curve on an Intel Xeon E5420 processor. The template-fitting code returned the estimated parameters of the top three best-fitting templates as well as the features used in the random forest classification detailed in Section 4.4 and Table 4. The computation time for fitting the template and calculating features for ~7 × 10⁵ light curves was ~44 K CPU hours. For comparison to a similar analysis, our algorithm is >9× faster than the template-fitting methods used by Sesar et al. (2017), which required ~30 minutes per star.

There are several other well-documented methods available in the literature for identifying RRL in multiband data (e.g., VanderPlas & Ivezić 2015; Hernitschek et al. 2016 and Sesar et al. 2017), which yield excellent results for data sets with an average of 35 or more observations per light curve. We present this algorithm as an alternative to these other methods for especially sparse multiband data sets (see Section 5.3 for a discussion of the observational limitations). The template and the associated fitting code are available at https://github.com/longjp/rr-templates, archived in Zenodo (doi:10.5281/zenodo.1252333), and further described in Appendix D.

## 4.3. Feature Selection

While it is possible to identify RRL by visually inspecting their phase-folded light curves, the sheer volume of light curves

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Top: poorly sampled DES light curve of an RRab originally discovered by Sesar et al. (2010; same star as Figure 5). Bottom: phased light curve of the same source with the period correctly estimated by our algorithm. Note: photometric uncertainties are smaller than the plotting symbols.</figcaption></figure>

<table><caption>**Table 3**
Period and Distance Estimation Accuracy</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">% of RRab within</th><th rowspan="2">σ<sub>parameter</sub></th></tr><tr><th>1%</th><th>3%</th></tr></thead><tbody><tr><td>ΔP/P<sub>prev</sub></td><td>88.89</td><td>89.74</td><td>6.81%</td></tr><tr><td>ΔD/D<sub>prev</sub></td><td>44.64</td><td>81.11</td><td>2.83%</td></tr></tbody></table>

Note. $P_{\text{prev}}$ and $D_{\text{prev}}$ are the values reported in Sesar et al. (2010) and Süveges et al. (2012).

in our data set makes this classification method unfeasible. Instead, we computed numerical features to describe the behavior of the light curves. To assess the specific parameter space occupied by RRab, we compiled a training set consisting of all the cross-matched labeled objects from S82 that passed the initial variability cuts (discussed in Section 3.2). This left us with a training set of 234 RRab, 57 RRc, 5196 other variable objects, and 3004 “standard” sources. Because we only aimed to identify RRab, we chose a simple identification scheme with two classes: RRab and non-RRab. This resulted in an RRab class with 234 objects and a non-RRab class with 8257 objects.

With the goal of identifying RRab, we chose features that were motivated by how well the light curves fit the RRab template and other observed properties of RRab. As demonstrated in Figure 3, RRab have relatively large $log_{10}(\chi^2_{\nu,b})$ compared to most of the other objects in our sample. So, to quantify the base variability of the light curve while ignoring spurious signals or missing observations in any particular band, we included the median of this value calculated across all five filters as a feature. Additionally, most true RRab should fit the template with small residuals, so we quantified the quality of the best template fits using the RSS per degree of freedom, $RSS_{\text{dof}} = RSS/(N_{\text{obs}} - 4)$. Then, to amplify the separation
<table><caption>**Table 4**
Random Forest Features</caption><thead><tr><th>Feature Name</th><th>Description</th><th>Importance</th></tr></thead><tbody><tr><td>lchi_med</td><td>Median log<sub>10</sub>(χ<sup>2</sup><sub>ν,b</sub>) value across all bands</td><td>0.2232</td></tr><tr><td>amp_rss_0</td><td>Best-fitting amplitude divided by the best-fitting RSS/dof</td><td>0.1942</td></tr><tr><td>f_dist1_0</td><td>Closest distance of the best-fitting period/amplitude to the Oosterhoff I relation from Sesar et al. (2010)</td><td>0.1591</td></tr><tr><td>rss_dof_0</td><td>RSS of the best-fitting template per degree of freedom</td><td>0.1571</td></tr><tr><td>f_dist2_0</td><td>Closest distance of the best-fit period/amplitude to the Sesar et al. (2010) curve dividing the Oosterhoff I and II groups</td><td>0.1025</td></tr><tr><td>amp_0</td><td>Best-fitting amplitude</td><td>0.0792</td></tr><tr><td>rss_lchi_med</td><td>RSS of the best-fitting template divided by the median log<sub>10</sub>(χ<sup>2</sup><sub>ν,b</sub>) value across all bands</td><td>0.0502</td></tr><tr><td>κ_0</td><td>Best-fitting von Mises-Fisher concentration parameter of phases in the folded light curve</td><td>0.0345</td></tr></tbody></table>

Note. All of these featured values for each candidate are included in the machine-readable version of Table 6.

provided by these two characteristics, we divided the RSS<sub>dof</sub> by the median log<sub>10</sub>(χ<sup>2</sup><sub>ν,b</sub>) to form another feature.

To take advantage of the distinctive amplitude ranges of RRab, we also selected the amplitude of the best-fitting template as a feature. We then created a new feature by dividing this amplitude by the RSS<sub>dof</sub>, expecting that the large amplitudes and excellent template fits of RRab would clearly distinguish them from other objects. We can take advantage of these amplitudes again by evaluating how closely each object matches the observational trends shown by RRab in the first two Oosterhoff groups (see the introduction for a brief description). To measure how closely the objects’ estimated template parameters matched these trends for RRab, we calculated the distance of the object in period-amplitude space from the Oosterhoff I relation measured in Sesar et al. (2010) and their shifted curve, which separates the Oosterhoff I and II populations (see their Figure 16 and our Figure 12).

Our last feature attempted to quantify the phase distribution of the observations in each light curve. Period-finding algorithms often recover periods at common aliases, sometimes resulting in light curves with many of their observations clustered near a particular phase. Because light-curve phases are periodic, the two-dimensional case of the von Mises-Fisher distribution (Fisher 1953; Jupp & Mardia 1989) is a good approximation. This distribution can be written as

$$f(x) = \frac{e^{\kappa \cos(x - \mu)}}{2\pi I_0(\kappa)}, \quad (5)$$

where $κ$ is the concentration parameter, $μ$ is the mean, and $I_0(κ)$ is the modified Bessel function of the first kind at order 0. The von Mises-Fisher distribution is akin to a Gaussian distribution wrapped around a circle, where the $κ$ parameter is analogous to the inverse of the variance (see Section 2.2 in Sra 2016 for a more detailed description). We calculated this concentration parameter $κ$ for each light curve folded over the best-fitting template period. Light curves with observations highly clustered in phase will have very large values of $κ$, aiding in the rejection of objects with aliased periods.

Although our choice of classifier is generally robust against non-informative features, we limited our features to these to make the classifier results easier to interpret. The features are summarized and ranked by their importance, or how much they contributed to splitting the data across all decision trees in our classifier,<sup>55</sup> in Table 4. These features are shown in Figure 7.

The development of additional features to further separate the classes will be explored in future work.

### 4.4. Random Forest Classifier

To identify likely RRab automatedly and consistently, we trained a random forest classifier (Amit & Geman 1997; Breiman 2001) using these features. The random forest is a machine-learning algorithm that predicts classes for data by combining results from a “forest” of decision trees. Each decision tree consists of a series of nodes where the data are split into subgroups based on the values of a random subset of their features, or characteristics. Before the random forest can make accurate predictions, it must be trained to recognize the trends in features that correspond to different classes. Thus, one needs a labeled training set to build the random forest. Each decision tree uses the labels to build a series of boundaries in feature space that divide the data into their correct classes. Once trained, the random forest algorithm assigns a score to unlabeled data based on the proportion of trees that identify them as a particular class. Random forest classifiers have been extremely successful in variable star classification (see Richards et al. 2011 for a comparison with other machine-learning techniques), even in the case of sparsely sampled Pan-STARRS PS1 light curves (Hernitschek et al. 2016). Thus, the random forest was a natural choice of classifier for this study.

We created the classifier using the `RandomForest` package available in `scikit-learn` (Pedegrosa et al. 2012). To prevent overfitting to our small training set and to ensure repeatability, the classifier contained 500 trees with a maximum depth of 5 and used a random seed of 10.

We assessed the performance of our classifier by estimating the purity (the percentage of objects classified as RRab that were truly so) and the completeness (the percentage of real RRab that were identified as such) as a function of the class score reported by the random forest. The purity and completeness were estimated using a five-fold cross-validation technique, where the data were divided into five test groups and classified based on a classifier trained on the other four groups. The classifier correctly identified 190 of the 234 RRab used to train the random forest as such with a score threshold of ≥35%. As shown in Figure 8, defining RRab as objects with a score ≥35% yields a purity of 85% and a completeness of 76%. A common metric used to assess the performance of a classifier is the area under the curve (AUC) shown in Figure 8, which we find to be 0.864. The purity and completeness calculated at other score thresholds are listed in Table 5. We include all other objects with lower scores in our catalog so that other score

<footnote n="55">55 See Section 1.11.2.5 in the `scikit-learn` documentation for more details.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Features used to identify RRab plotted for the training set. Red points denote previously identified RRab while black X's are non-RRab. While the RRab clearly occupy a specific region in this feature space, they are not linearly separable.</figcaption></figure>

thresholds can be specified by interested readers. Although an incorrect period estimate led to a worse RSS value for the fit, some RRab with incorrect period estimates were still correctly identified as such by the classifier.

Because our training set was mostly composed of nearby RRab confined to a small region in the survey footprint, it is imperative to assess our classifier performance using other samples of RRL. To this end, we cross-matched our sample with external surveys in Section 5.2 and applied our method

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Purity/completeness curve for the random forest classifier trained on cross-matched objects in S82. The black star denotes a classifier-reported score of 0.35, where the purity is ~85% and the completeness is ~76%. The area under the curve is 0.864.</figcaption></figure>

<table><caption>**Table 5**
Estimated RRab Selection Purity and Completeness</caption><thead><tr><td>Score Threshold</td><td>Purity</td><td>Completeness</td></tr></thead><tbody><tr><td>0.00</td><td>0.043</td><td>0.983</td></tr><tr><td>0.05</td><td>0.402</td><td>0.928</td></tr><tr><td>0.10</td><td>0.567</td><td>0.886</td></tr><tr><td>0.15</td><td>0.661</td><td>0.852</td></tr><tr><td>0.20</td><td>0.727</td><td>0.840</td></tr><tr><td>0.25</td><td>0.773</td><td>0.819</td></tr><tr><td>0.30</td><td>0.815</td><td>0.798</td></tr><tr><td>0.35</td><td>0.845</td><td>0.756</td></tr><tr><td>0.40</td><td>0.877</td><td>0.722</td></tr><tr><td>0.45</td><td>0.896</td><td>0.693</td></tr><tr><td>0.50</td><td>0.922</td><td>0.651</td></tr><tr><td>0.55</td><td>0.955</td><td>0.630</td></tr><tr><td>0.60</td><td>0.953</td><td>0.596</td></tr><tr><td>0.65</td><td>0.956</td><td>0.554</td></tr><tr><td>0.70</td><td>0.968</td><td>0.525</td></tr><tr><td>0.75</td><td>0.973</td><td>0.470</td></tr><tr><td>0.80</td><td>0.971</td><td>0.432</td></tr><tr><td>0.85</td><td>0.988</td><td>0.348</td></tr><tr><td>0.90</td><td>0.982</td><td>0.231</td></tr><tr><td>0.95</td><td>1.000</td><td>0.071</td></tr></tbody></table>

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

to simulated RRab light curves at fainter magnitudes in Section 5.3.

# 5. Catalog Description

## 5.1. Visual Validation

We applied the classifier to the ~7 × 10⁵ objects with template fits and found 8026 RRL candidates with a score ≥0.35. Although most of our candidates were indeed RRL found in other surveys, there were still some nonstellar interlopers in the sample, due to the lenient initial cuts on the shape of the photometric point-spread function (Section 2.2). Thus, we visually inspected all RRab candidate light curves and their DR1 coadded images. After visually validating the candidates and removing any objects with non-RRL classifications in the Simbad database (Wenger et al. 2000), 1786 objects were discarded and 5783 RRL candidates remained in the catalog, with the rest too ambiguous to confirm. A sample of
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Sample of DES coadded images and representative light curves of visually accepted RRL candidates with classifier scores exceeding 0.94, labeled with their Y3Q2 ID number. The observations and templates are colored by filter using the same convention as Figure 6.</figcaption></figure>

visually verified candidates with high ($p > 0.94$) and low ($p < 0.36$) classifier scores is shown in Figures 9 and 10, respectively. The most typical contaminants in the sample were extended sources. Given our lenient selection criteria described in Section 2.2, it is not surprising that some of these objects made it into our final sample. Examples of candidates that were rejected by visual inspection as being extended sources are shown in Figure 11. The full catalog of candidates, their best-fit parameters, and their features are available at https://des.ncsa.illinois.edu/releases/other/y3-rrl. A sample view of this catalog is shown in Table 6. Appendix E contains a detailed description of these data products.

Although the light curves have been visually inspected, further photometric observations of some extremely poorly sampled candidate light curves would be useful to confirm their classification. One may wish to remove these more uncertain candidates from their analysis by only considering objects with a larger minimum number of observations. Some of these candidates have gaps in observations near their maximum and minimum brightness, providing poor constraints on their estimated amplitudes and mean magnitudes. Therefore, we assigned a flag to each object based on how its phase-folded light curve is sampled. An object with fewer than two observations near its minimum brightness ($0.55 \le \text{phase} < 0.87$, which we chose to encompass both the near-constant portion of the light curve where other authors chose their minimum phase (e.g., Vivas et al. 2017) and the 10% quantile of template magnitudes) will receive a “flag_minmax” value of 1, while an object with $<2$ observations near its maximum brightness ($0.96 \le \text{phase} < 1$ or $0 \le \text{phase} < 0.05$ corresponding to the 90% quantile in template magnitudes) receives a “flag_minmax” value of 2. Objects missing observations near both of these receive a flag value of 3.

Figure 12 shows a Bailey (period–amplitude) diagram of the candidates rejected by the classifier or our visual checks plotted in black, visually accepted candidates shown in red, and ambiguous candidates in blue. We plot the Oosterhoff I (Oosterhoff 1939) relation and the curve dividing groups I and II from Sesar et al. (2010), which we used in the calculation of our features, in solid and dashed black lines. Many of these ambiguous candidates are likely RRab, but cannot be classified as such with high confidence in this work. Due to the sparse nature of our observations, we cannot detect amplitude modulations such as those arising from the Blažko effect (Blažko 1907), although we did recover 5 out of the 12 Blažko RRL previously identified in the Catalina Surveys (Drake et al. 2014, 2017).

## 5.2. Comparison with Overlapping Catalogs

The DES footprint has significant overlap with other surveys, such as Gaia (Gaia Collaboration et al. 2016), Pan-STARRS (Chambers et al. 2016), the Catalina Surveys (Drake et al. 2009), and the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018). We used our cross-matches with these external RRab catalogs to independently assess the performance of our algorithm at the different magnitude ranges probed by these surveys. We used the SkyCoord package in astropy (The Astropy Collaboration et al. 2018) to select matches within 1" of DES objects while removing duplicates. Details for each individual survey are presented in the following paragraphs and summarized in Table 7, while Figure 13 shows the respective overlaps with DES.

Clementini et al. (2018) found over $1.4 \times 10^5$ RRL in Gaia DR2, including $\sim 5 \times 10^4$ that were previously unknown. These variables were identified from multiband (G, $G_{BP}$, $G_{RP}$)
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Sample of DES coadded images and representative light curves of visually accepted RRL candidates with classifier scores below 0.36, labeled with their Y3Q2 ID number. The observations and templates are colored by filter using the same convention as Figure 6.</figcaption></figure>

light curves that had at least 12 observations in *G* (see Figure 10 in Holl et al. 2018). While the *Gaia* temporal coverage is very uneven, their RRL catalog spans the entire sky (see Figure 26 in Clementini et al. 2018) and has high purity (~91%), making it an excellent independent check of our method at brighter magnitudes. A total of 4609 *Gaia* DR2 RRab were present in our initial stellar catalog (Section 2.2) and 3227 (~70%) were identified as such. To assess this recovery another way, if we create a purity versus completeness curve from these cross-matches like the one shown in Figure 8, we find an AUC of 0.727. As we have significantly fewer single-band observations than *Gaia* DR2, it is not surprising that we do not recover all of their RRab.

We also searched for RRab discovered in Pan-STARRS PS1 by Sesar et al. (2017). Like DES, Pan-STARRS has sparsely sampled multiband light curves and Sesar et al. (2017) employed a similar template-fitting method to identify these variables. However, Sesar et al. (2017) used the final data release of PS1 with an average of 67 observations per object (compared to our median of 18). We adopted their suggested ab_score cut of 0.8 to select only RRab. As Pan-STARRS primarily surveyed the northern hemisphere, we found just 1021 RRab in our initial stellar catalog, but we identified 805 (~79%) as such, with an AUC of 0.681. As the Pan-STARRS light curves are the most similar to the DES Y3 ones out of all the external catalogs under consideration, our similar classification results show that our approach is similarly effective to the methods used by Sesar et al. (2017).

The Catalina Surveys RRL catalog (Drake et al. 2013a, 2013b, 2014; Torrealba et al. 2015; Drake et al. 2017) is based on a wide-field (26,000 deg²) time series survey that probes the variable sky to a depth of V ~ 19–20 mag. The observations are unfiltered and collected in sequences of four images equally

spaced over 30 minutes in each pointing (Drake et al. 2009). After several years of operation, the Catalina Surveys have over 200 observations for most of their variables (Drake et al. 2014), which makes the catalog largely complete. Given the limited magnitude overlap between the Catalina Surveys and DES, we only found 1463 of their 32,775 RRab in our initial stellar catalog, but we identified 1185 (~81%) as such, with an AUC of 0.733.

ATLAS, a planetary defense initiative with a high cadence well suited for variability studies, recently released its first catalog of variable stars (Heinze et al. 2018). Thus far, ATLAS has at least 200 observations across two filters (c, o) over one-fourth of the sky. We select RRab stars from the ATLAS DR1 variable star catalog using the suggested CasJobs query in Appendix 10.2 of Heinze et al. (2018). As ATLAS is based in the northern hemisphere and quite shallow compared to DES (r ≈ 20 mag), we only have 484 of their 21,061 RRab in our initial stellar catalog but identify 391 (~81%) as such, with an AUC of 0.635. This recovery rate is quite similar to the ones for Pan-STARRS and the Catalina Surveys.

In addition to searching for RRab candidates with previous identifications from the aforementioned wide-field surveys, we also checked for overlaps near the Magellanic Clouds (Soszyński et al. 2016), the Fornax dSph (Bersier & Wood 2002), the Sculptor dSph (Martínez-Vázquez et al. 2016), in the General Catalogue of Variable Stars (Ram et al. 2017), and in the SIMBAD database (Wenger et al. 2000). To the best of our knowledge, and based on publicly available catalogs, 1603 (nearly 28% of our sample) are newly discovered RRab candidates. Although the external catalogs under consideration are not complete, the fraction of their RRab recovered by our analysis is consistent with our estimate of ~75% completeness. Our method is just as effective (if not more so) at recovering
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Sample of DES coadded images and representative light curves of visually rejected candidates (extended sources or possible supernova) that passed the classifier score threshold, labeled with their Y3Q2 ID number. The observations and templates are colored by filter using the same convention as Figure 6.</figcaption></figure>

RRL from suboptimally sampled data than the methods used in comparable surveys.

Although we recover most of the RRab in the aforementioned overlapping catalogs, we can see from the AUC of each of these that there is a marked degradation in our algorithm's performance when applied to light curves outside our S82 training set. Thus, we use their AUC values to construct a confidence interval for the performance of our classifier. With the AUC of the training set and all four of these external cross-matches, we find a mean AUC of 0.728 with a standard deviation of 0.077. From this, we can determine that our classification methods have a lower efficiency for fainter RRab. Unfortunately, we do not have well-characterized training data in a comparable filter system for fainter RRab, so we tested this with simulated light curves.

### 5.3. Estimated Recovery Rates and Uncertainties from Simulated Data

To estimate the robustness of our results for the noisier photometry at fainter magnitudes, we followed a method similar to Medina et al. (2018) and applied our method to simulated light curves with known light curve parameters in the DES filter system. We created the simulated light curves by sampling the smoothed templates of Sesar et al. (2010) in *gatspy* (VanderPlas & Ivezić 2015) with the DES cadence from different areas of the survey. We shifted these light curves to various distances by adding the appropriate distance modulus and inserting scatter in the observations based on the magnitude-dependent uncertainty relations we found in Section 3.1 (shown in Figure 2). Appendix F contains further details on the construction of these simulated light curves.

Figure 14 shows the recovery rates of both the classifier and the period as a function of magnitude and total number of

observations. As expected, the recovery rate of our algorithm decreases significantly with increasing distance modulus. This is mostly due to the larger photometric uncertainties and fewer observations due to the brighter limiting magnitudes for the redder bands (see Section 2.1). We see that the accuracy of the period estimation decreases following the trend of increasing photometric errors shown in Figure 14, and dramatically improves with increasing total number of observations up to $N \sim 20$. As expected, the RRab classification accuracy follows a similar trend. We find that our template fitting recovers the true period to within 1% for 95% of the simulated light curves with $N = 20$ observations.

Beyond assessing our classifier performance with these simulated light curves, we can also use them to estimate the uncertainties of the best-fitting template parameters. To make sure we treat light curves with especially poor phase coverage separately, we divided the simulated light curves into groups based on their “flag_minmax” values (described in Section 5.1). Then, we subdivided those into bins of two observations and 0.5 mag wide in $N$ and $\mu$, respectively. In each of these bins, we calculate the fraction of light curves with period estimates within 1% of their input values for each phase-sampling group. To quantify the uncertainty of the period estimates, we calculated the standard deviation of $\Delta P/P = (P_{est} - P_{true})/P_{true}$, where the “est” subscript represents the parameter estimate from the template fitting and “true” represents the input value of the simulated light curve. Likewise, we calculated $\Delta a = a_{est} - a_{true}$ to quantify the uncertainty of the amplitude estimates. The number of light curves included in each bin differs widely, so we estimate the spread of these uncertainty values within each subgroup with jackknife resampling. These results are shown in Figure 14.
<table><caption>**Table 6**
DES RRab Candidates</caption><thead><tr><th>DES Y3Q2 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><th>p</th></tr></thead><tbody><tr><td>11136400113264</td><td>0.000107</td><td>-59.559187</td><td>17.657</td><td>17.500</td><td>17.547</td><td>17.535</td><td>17.588</td><td>0.6415</td><td>0.551</td><td>16.90</td><td>0.436</td></tr><tr><td>10646400013129</td><td>0.013042</td><td>-2.430057</td><td>17.234</td><td>17.052</td><td>17.099</td><td>17.103</td><td>17.252</td><td>0.4938</td><td>1.151</td><td>16.66</td><td>0.733</td></tr><tr><td>11004800140792</td><td>0.131498</td><td>-41.482218</td><td>19.229</td><td>19.096</td><td>19.125</td><td>19.195</td><td>19.145</td><td>0.5893</td><td>1.137</td><td>18.69</td><td>0.568</td></tr><tr><td>11108800089990</td><td>0.134204</td><td>-54.295118</td><td>15.754</td><td>15.625</td><td>15.557</td><td>15.554</td><td>15.571</td><td>0.6698</td><td>0.469</td><td>14.97</td><td>0.659</td></tr><tr><td>10595200009863</td><td>0.283437</td><td>1.178535</td><td>17.986</td><td>17.922</td><td>17.908</td><td>17.868</td><td>17.869</td><td>0.5480</td><td>1.033</td><td>16.96</td><td>0.886</td></tr></tbody></table>

Notes. DES Y3Q2 ID: DES Y3Q2 QUICK_OBJECT_ID number. α: R.A. δ: decl. (grizY): mean extinction-corrected magnitude. P: best-fit period. A_g: best-fit amplitude in DES g. μ: best-fit distance modulus. p: RRab score assigned by the classifier. The full version of this catalog, including feature values and cross matching information, is available in the online data products at https://des.ncsa.illinois.edu/releases/other/y3-rrl.

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

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Bailey diagram of template-estimated amplitudes and periods for objects that passed the initial variability cuts in black and visually accepted RRab identified by our classifier in red. Ambiguous candidates that could not be visually accepted are plotted in blue. We overplot the Oosterhoff I relation and the curve dividing the Oosterhoff I and II populations from Sesar et al. (2010) in black solid and dashed lines, respectively. The abundance of objects with periods of P = 0.5 days denotes a common alias of the 1 day rotation period of Earth.</figcaption></figure>

Other than fluctuations due to the small sample sizes in some of the bins, these values follow expected trends. When there are fewer observations to constrain the parameter values during the template fitting, both the period and the amplitude are more uncertain, with these values beginning to stabilize around N = 20 observations. In distance space, the parameter estimates are generally low until μ ≈ 20, where the brighter detection limits of the redder filters decrease the number of observations in the light curves. We have very few simulated light curves that are missing observations near their maximum only or both their maximum and minimum (the blue and orange points in Figure 14), so we cannot draw any definitive conclusions about the effect of phase sampling on the estimation of these parameters. We assign these parameter uncertainties to the real RRab candidates based on the best-fitting third-degree polynomial to the trends in N observations for all simulated light curves. We do not assign uncertainties to objects with N > 43 observations, due to a lack of simulated data with that sampling. We also do not report these uncertainties for objects not identified as RRab by the classifier because these simulated light curves do not accurately represent the behavior of non-RRab. The coefficients of the best-fitting polynomials are included in Table 8, and the uncertainties are included in the full catalog described in Appendix E.

The uncertainty of the remaining parameter φ is significantly more difficult to constrain. Phases for individual observations

in the folded light curves are calculated using phase = (MJD/P)mod 1. Any small offset in the period will compound over successive pulsations and result in a phase offset that varies over time. Even simulated light curves with ΔP/P < 0.0005 (a difference of <1 minute) can yield Δφ ≈ 0.5 after three years when compared to the phases calculated using the input period. Thus, we do not report these uncertainties in φ as they require a level of period precision we do not attain even in light curves with N > 20. We caution against using the phases reported here for purposes other than plotting the template curves.

## 5.4. Uncertainties in the Distance Moduli

Because the absolute magnitudes of RRL depend on their metallicities (see Figure 14 in Marconi et al. 2015), the fixed [Fe/H] of our RRL model contributes systematic uncertainty to our distance estimates. Although the abundances of the individual RRL in our catalog are unknown, we can approximate the size and direction of this effect by comparing our results to those from an external catalog with metallicity measurements. The Catalina Surveys catalog of Torrealba et al. (2015; hereafter T15) is convenient for this purpose because it has photometric estimates of [Fe/H] and a significant overlap with the DES survey footprint (although it has a brighter magnitude limit; see Section 5.2). We calculated the difference in distance moduli for 521 RRL in common between both catalogs, Δμ = μ<sub>DES</sub> − μ<sub>T15</sub>, which we plot as a function of [Fe/H] in Figure 15.

We split the sample into bins of 0.1 dex in metallicity and perform an iterative 3σ clip from the median value using the `sigmaClip` function in `astropy.stats`. We fit a linear relation between Δμ and [Fe/H]:

$$ \Delta\mu(\text{DES} - \text{T15}) = (-0.058 \pm 0.003) + (0.168 \pm 0.009) ([\text{Fe/H}] + 1.5). \quad (6) $$

The root-mean-square error of the fit is 0.06 mag, consistent with the standard deviation of ΔD/D between this work and Sesar et al. (2010) listed in Table 3. Thus, we estimate our statistical uncertainty in distance moduli to be σ<sub>stat</sub> = 0.06 mag.

As is evident from Figure 15, our algorithm systematically underestimates distances for very metal-poor RRL and over-estimates distances for more metal-rich RRL. While this cross-matched sample illustrates how much of an effect an RRab’s metallicity has on the accuracy of our distance estimates, we caution against using this sample to derive a metallicity
<table><caption>**Table 7**
Description of Selected External RRL Catalogs and Their Overlap with DES</caption><thead><tr><th rowspan="2">Survey</th><th rowspan="2">Area<br>(sq deg)</th><th rowspan="2">Filters</th><th rowspan="2">Depth</th><th rowspan="2">Observational<br>Cadence</th><th rowspan="2">N<sub>obs</sub></th><th colspan="3">RRab</th></tr><tr><th>Total</th><th>in DES</th><th>% Found</th></tr></thead><tbody><tr><td>SDSS Stripe 82</td><td>~300</td><td>ugriz</td><td>g &le; 21</td><td>Most observed every 2 days</td><td>70–90</td><td>447</td><td>238</td><td>75</td></tr><tr><td>Gaia DR2</td><td>all sky</td><td>G<sub>BP</sub>, G, G<sub>RP</sub></td><td>G ~ 21</td><td>Uneven, follows Gaia scanning law</td><td>12–240</td><td>140,000</td><td>4609</td><td>70</td></tr><tr><td>Pan-STARRS PS1</td><td>~30000</td><td>grizY</td><td>r &le; 21.5</td><td>Two same-band obs. sep. by 25 minutes</td><td>~67</td><td>35,000</td><td>1021</td><td>79</td></tr><tr><td>Catalina Surveys</td><td>~9000</td><td>unfiltered</td><td>V &le; 19–20</td><td>Four obs. within 30 minutes</td><td>&ge;200</td><td>32,775</td><td>1463</td><td>81</td></tr><tr><td>ATLAS</td><td>~13000</td><td>c, o</td><td>r ~ 20</td><td>4 &times; per night</td><td>~200</td><td>21061</td><td>484</td><td>81</td></tr><tr><td>DES Single Epoch</td><td>~5000</td><td>grizY</td><td>g ~ 23.5</td><td>Irregular</td><td>~50*</td><td>5783</td><td>5783</td><td>...</td></tr></tbody></table>

Note. The details of DES are listed for comparison. (*): by the end of the survey (Y6).

correction to our distance estimates. The T15 RRL cover most of the spatial DES footprint, but our RRab catalog extends to fainter magnitudes than the magnitude limit probed by the T15 sample, which means we cannot assume that the metallicity distribution of this sample accurately reflects that of the entire catalog. However, if we use this subsample and assume that the stellar halo metallicity distribution function is represented by a Gaussian with a mean of [Fe/H] = -1.5 dex and standard deviation of σ = 0.3 dex (Ivezić et al. 2019) as in Sesar et al. (2017), we find a 1σ systematic uncertainty of σ<sub>[Fe/H]</sub> ≈ 0.05 mag in distance modulus. If we follow Medina et al. (2018) and estimate distance offsets for a metallicity shift of ±0.5 dex and ±1.0 dex in this subsample, we find a change in distance modulus of 0.08 mag and 0.17 mag, respectively. Again, we caution that the true distribution of metallicities in our full catalog is unknown, so these values are merely representative of the systematic uncertainty that would apply to particular stellar populations in the Milky Way halo. However, given our lack of metallicity information, we cannot quantify these systematic offsets without making such assumptions.

Another contribution to the systematic uncertainty we have not previously considered is the RRL evolution off the horizontal branch. We adopt the value σ<sub>(V)</sub> = 0.08 mag, which Vivas & Zinn (2006) estimated from RRL in globular clusters (see their Section 4 and Figure 4). Adopting the halo metallicity distribution from Sesar et al. (2017), we add both sources of uncertainty in quadrature to arrive at σ<sub>sys</sub> = (σ<sup>2</sup><sub>[Fe/H]</sub> + σ<sup>2</sup><sub>&lt;V&gt;</sub>)<sup>1/2</sup> ≈ 0.09 mag.

We verify this estimate of systematic uncertainty by comparing our estimated distance moduli for various Milky Way satellites with previously published results. For the Fornax dSph, which has a horizontal branch [Fe/H] ≈ -1.8 dex (Rizzi et al. 2007), our median distance modulus is 0.05 ± 0.04 mag closer than the values of μ = 20.72 ± 0.04 and μ = 20.72 ± 0.06 mag found by Clementini et al. (2006) and Rizzi et al. (2007), respectively. In the case of the Sculptor dSph, we find a median distance modulus 0.13 ± 0.04 mag closer than the value of μ = 19.62 ± 0.04 mag from Martínez-Vázquez et al. (2016). This larger difference is likely due to the large spread in metallicity exhibited by Sculptor's stellar populations, -2.3 ≲ [Fe/H] ≲ -1.5 dex (Martínez-Vázquez et al. 2016). Our distance estimate to the LMC, based on the RRL we found in its outskirts, is 0.12 ± 0.09 mag closer than the μ = 18.52 ± 0.09 mag found by Clementini et al. (2003). The LMC also has a spread of metallicities for HB stars, centered on [Fe/H] ≈ -1.5 dex with a dispersion of 0.4 dex (Clementini et al. 2003; Gratton et al. 2003). We expect that replacing the template's M<sub>B</sub>(P) relation with a calibrated P-L-Z or P-L-C relation in the DECam filters (K. Vivas et al. 2019, in preparation) will significantly reduce these offsets.

In summary, our distance moduli have 1σ statistical and systematic uncertainties of 0.06 and 0.09 mag, respectively. The equivalent distance uncertainties are ~2.8% (stat) and ~4.2% (sys).

# 6. Discussion

## 6.1. Detection Biases

The strength of DES lies in its wide-field coverage and depth, but the results presented here are limited by the low number of multiband observations. Figure 16 displays histograms of the total number of observations for all objects in the stellar sample (black), all objects passing variability cuts (blue), and all RRab candidates (red). The median number of total observations for each group, marked by a short colored line segment, is 10, 18, and 17, respectively. Note that most of our RRab have fewer observations than the N ~ 20 observation threshold we saw from the simulation results in Figure 14. As future DES data releases will have an increased number of observations, we expect to find more RRL and have a more robust classification of the candidates presented here.

We note that the light curves used in this analysis typically had fewer total observations than the number expected from three years of DES data. We suspect that the total number of observations for the “stellar” sample is skewed by objects near the detection limits of DES, which suffer from noisy photometry and likely have few overlapping observations across all five filters. This low number of observations is also a result of the stringent quality cuts we applied on the single-epoch photometry in Section 2.2. In future work, we aim to be more judicious in applying our photometric quality cuts so that we do not discard observations unnecessarily.

To verify that our sample is not affected by spatial fluctuations in the number of observations, we calculated the median number of total observations in each HEALPix of our Y3Q2 stellar data set. We show the median and the standard deviation of the total number of observations of light curves in each HEALPix in Figure 17. As expected, regions with the lowest number of observations fall near the edges of the survey footprint. Regions that have a median number of observations ≥25 correspond to the Science Verification region, in which 50 observations were made in the first year to demonstrate year 5 depth and the DES Supernova fields, which are observed roughly weekly (e.g., Dark Energy Survey Collaboration et al. 2016). The linear patterns of constant R.A. are a result of the survey observation strategy (Abbott et al. 2018). Beyond these patterns, the DES photometry suffers in photometric completeness in crowded stellar fields near the central regions of nearby dSph galaxies and globular clusters, thus our catalog also suffers in completeness near those regions.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Histograms of magnitudes of RRab stars from external catalogs cross-matched with our DES initial stellar catalog, as a function of the extinction-corrected weighted-average coadded DES *r* magnitude. Top left: Gaia DR2. Top right: Pan-STARRS PS1. Bottom left: Catalina Surveys DR2. Bottom right: ATLAS. Blue curves show the RRab from each catalog that were present in our sample before applying any cuts, while red curves show those that were identified as RRab in our analysis. Black curves show the distribution of DES RRab candidates and are the same in all panels. The overdensities at *r* ≈ 18.8 and *r* ≈ 21.2 correspond to the LMC outskirts and the Fornax dSph. Our catalog is deeper than the others by ~1, 1, 2, and 4.5 mag, respectively.</figcaption></figure>

Otherwise, the survey coverage is fairly uniform, and we do not expect large-scale trends in RRL detection outside of these fields of larger-than-average observation counts and dense stellar populations. We expect the addition of DES Y4–Y6 data to increase our detections of RRL considerably.

Some additional biases in our RRab sample are results of choices made to exclude non-RRab from our analysis. While we weighted our initial variability cuts by the photometric errors to make the cuts robust against spurious observations (see Section 3.2), using these error-weighted metrics biased our variable sample against RRab with smaller amplitudes located at larger distances. We also excluded some real RRab from our sample by limiting the period range to 0.44 days $\le P \le$ 0.89 days to avoid the common 1 day alias.

## 6.2. Spatial Distribution of the Candidates

The spatial distribution of the 5783 visually validated RRab candidates is shown in Figure 18. We also plot these candidates as a function of their heliocentric distance in Figure 19. In both figures, the overdensities of RRab candidates associated with (in order of decreasing heliocentric distance) the Fornax dSph, the Sculptor dSph, and the outskirts of the Large Magellanic Cloud are easily visible. We expect that the inclusion of a metallicity term into the model combined with additional epochs of DES observations in the next release of this catalog will enable further characterization of these and other substructures.

One of the largest strengths of the DES data set is its depth (see our comparisons to other wide-field surveys in Figure 13). This is extremely valuable for our understanding of the outer halo as the current census of RRL known beyond 100 kpc falls short of the thousands predicted by simulations (Sanderson et al. 2017). In this work, we identified 800 RRab candidates beyond 100 kpc (most of which have been previously discovered) and

eight RRab candidates beyond 200 kpc, all of which are new discoveries. The coadded images and light curves for the candidates beyond 200 kpc are shown in Figure 20.

The three most distant visually verified RRab candidates in our sample have heliocentric distances of ~231.6, 223.0, and 221.3 kpc. While these three stars are the most distant RRab in the Milky Way to date, they are not the most distant RRL. Medina et al. (2018) recently found two RRc with larger distances (232.9 and 261.2 kpc) using data from the HiTS Survey (Förster & Maureira 2016). Even though the candidates in our RRab sample suffer from a small number of observations and require additional follow-up for confirmation, the fact that there are so many RRab beyond 100 kpc and three RRab beyond 220 kpc provide reasonable evidence that the Milky Way stellar halo extends at least out to 220 kpc. Future DES data releases and other upcoming deep surveys such as LSST will increase the census of known RRL at this distance, enabling further characterization of the outer halo.

## 6.3. Applicability to LSST

The next-generation large ground-based LSST (Ivezić et al. 2019) is set to begin full science operations in early 2023.<fnref n="56" /> The most current LSST “Baseline Cadence” for its Wide-fast-deep Survey (WFD), which covers 18,000 deg² of the sky and comprises ~85% of its total allocated observing time, is to image each field twice 40 minutes apart once every three days in a different filter. After 10 yr of operation, each field is expected to have a median of (62, 88, 199, 201, 180, 180) visits in *ugrizy* with a single-epoch depth of (23.14, 24.47, 24.16, 23.40, 22.23, 21.57) mag. Assuming these observations are spaced uniformly over 10 yr, one can expect most light curves

<footnote n="56">⁵⁶ https://www.lsst.org/about/timeline</footnote>
<figure><img src="image_17.png" /></figure>

**Figure 14.** Recovery rates and parameter uncertainties as a function of the number of observations in the light curves and distance modulus *μ*. Colored points denote the behavior of the recovery fractions or parameter offsets for light curves with the phase-sampling flags described in Section 5.1. Uncertainties on these values were estimated using jackknife resampling. The black points in the upper-right corner of the bottom panels show the representative width of the bins in each column. Dashed gray lines show the best-fitting third-degree polynomial to the trends shown by the combined simulated data used to assign uncertainties in the RRab catalog. The coefficients for these fits are listed in Table 8. First row: fraction of simulated RRab light curves that received a classifier score $\ge$0.35. Second row: fraction of estimated periods within 1% of the true input values. Third row: standard deviation of the percent difference in period. Bottom row: standard deviation of the offset in amplitude. Note: if not visible, uncertainties are smaller than the plotting symbols.

to have ~80 multiband observations within the first year (LSST Science Collaboration et al. 2017).<fnref n="57" />

Oluseyi et al. (2012) found from their analysis of simulated LSST data that reliable RRL period estimation will require several years of operation, however, several multiband techniques have been developed since their publication. For

instance, VanderPlas & Ivezić (2015) estimated accurate periods 64% of the time on downsampled S10 light curves with ~55 observations and Sesar et al. (2017) accurately estimated periods for 85% of their PS1 training set with ~67 observations. Our simulations show that our template-fitting method is capable of estimating the correct periods to within 1% for 95% of the light curves with 20 total observations. Thus, our algorithm would be effective to identify potential RRab candidates (which would need follow-up for confirmation) within the first year of LSST operations. After the first

<footnote n="57">⁵⁷ See the most current version of the draft of the LSST Observing Strategy white paper located at https://github.com/LSSTScienceCollaborations/ObservingStrategy.</footnote>
<table><caption>**Table 8**
Coefficients for Parameter Uncertainties</caption><thead><tr><th>Value</th><th>p<sub>0</sub></th><th>p<sub>1</sub></th><th>p<sub>2</sub></th><th>p<sub>3</sub></th></tr></thead><tbody><tr><td>σ(ΔP/P)</td><td>4.8585 × 10<sup>-1</sup></td><td>-4.5912 × 10<sup>-2</sup></td><td>1.5636 × 10<sup>-3</sup></td><td>-1.7574 × 10<sup>-5</sup></td></tr><tr><td>σ(Δa)</td><td>1.5333</td><td>-1.5101 × 10<sup>-1</sup></td><td>5.3006 × 10<sup>-3</sup></td><td>-6.1034 × 10<sup>-5</sup></td></tr></tbody></table>

Note. The best-fit third-degree polynomial is of the form $\sigma(\text{Value}) = p_0 + p_1 N + p_2 N^2 + p_3 N^3$.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Difference in distance moduli between our work and Torrealba et al. (2015, T15) plotted against their photometric metallicity estimates. Stars that passed the 3σ clip are plotted as black points, while the points that were removed are plotted as gray crosses. The best-fit linear relation fit to the clipped data is plotted in red. The root mean square error (RMSE) of the fit is 0.06.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Histograms of the total number of observations across all bands for objects in the stellar sample (see Section 2.2; black), objects that passed the variability cuts (blue), and the objects identified as RRab (red). Short line segments denote the median number of observations for each group (10, 18, and 17, respectively). Note that many of the RRab candidates have very few observations and would benefit from follow-up observations to confirm their nature.</figcaption></figure>

year, the light curves in the WFD survey will be adequately sampled to use other multiband methods available in the literature.

## 7. Summary

We have presented a new physically motivated general multiband RRab template and a computationally efficient fitting procedure. We combined this method with a random forest classifier to create a powerful technique that can robustly identify these variables even when fewer than 20 observations are available. Despite the poor cadence and sampling of DES data, we detected 5783 RRab candidates, 1603 (28%) of which are previously undiscovered to the best of our knowledge. The large quantity of RRL we recovered in common with overlapping external surveys such as Gaia DR2, Pan-STARRS, Catalina Surveys, and ATLAS provides strong evidence of the effectiveness of this algorithm. Although the number of

observations is relatively uniform across the survey footprint in the DES Year 3 data, time series analyses like these will benefit immensely from the additional observations in future data releases. We make the template, these catalogs, and the light curves of the RRab candidates and the training sample available to the scientific community for future studies. Our method is especially useful for other multiband data sets that were not specifically designed for time series analysis.

The authors would like to thank the anonymous referees for insightful and constructive comments. The authors acknowledge the Texas A&M University Brazos HPC cluster that contributed to the research reported here: brazos.tamu.edu. K.M.S. would like to thank Chris Stringer and Ryan Oelkers for helpful conversations. Part of this work was completed while coauthors Z.L. and J.P.L. were research fellows in the Astrostatistics program at the Statistics and Applied Mathematical Sciences Center (SAMSI) in fall 2016. Z.L. and J.P.L. gratefully acknowledge SAMSI's support. C.N. and J. L. Myron would like to acknowledge support from the NSF grant AST-1263034, "REU Site: Astronomical Research and Instrumentation at Texas A&M University." This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

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
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Spatial map of the median number of observations across all bands in each HEALPix. As expected, the regions with the fewest observations are near the edges of the survey region. The regions with a relatively large number of observations correspond to Science Verification or supernova fields. Outside these regions, the DES coverage is relatively uniform, but suffers from a small number of time series observations. Future studies of RRL will benefit from additional years of DES data.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Map of the 5783 visually accepted RRab candidates across the DES wide-field survey footprint. The RRab are marked by dots colored by distance modulus. Large Milky Way satellite galaxies are easily distinguishable by their overdensities of RRab. The outskirts of the LMC are located near (80°, -62°), the Fornax dSph is located near (41°, -34°), and the Sculptor dSph is located near (15°, -34°).</figcaption></figure>

University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical 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
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Radial distribution of DES RRab stars. The overdensities associated with (in order of decreasing heliocentric distance) the Fornax dSph, the Sculptor dSph, and the periphery of the Large Magellanic Cloud are easily distinguishable. Note: the Sculptor and Fornax galaxies appear elongated due to uncertainties in the RRab distance moduli.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** DES coadded images and representative light curves of visually accepted candidates beyond 200 kpc, labeled with their Y3Q2 ID number. The observations and templates are colored by filter using the same convention as Figure 6. These very distant candidates have very few observations and will benefit from future DES data releases.</figcaption></figure>

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
<figure><img src="image_24.png" /><figcaption>Figure 21. Examples of simulated RRL light curves labeled with the ID number of the SDSS light curve used to generate them. The observations and templates are colored by filter using the same convention as Figure 6.</figcaption></figure>

non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

*Facility*: CTIO:Blanco (DECam).

*Software*: astropy (The Astropy Collaboration et al. 2018), easyaccess (Carrasco Kind et al. 2019), gatspy (VanderPlas & Ivezić 2015), Matplotlib (Hunter 2007), scikit-learn (Pedegrosa et al. 2012), SIMBAD (Wenger et al. 2000).

## Appendix A

### RRL Model Assumptions

The form for the RRL model is

$$m_b(t) = \mu + M_b(\omega) + R_b E(B - V) + a\gamma_b(\omega t + \phi) \quad (7)$$

where the population parameters, common to all RRL, are

$M_b(\omega)$ = absolute magnitude in band b for an RRL with frequency $\omega$

$R_b$ = total-to-selective extinction coefficient for band b

$\gamma_b$ = shape of the RRL light curve in band b.

The object specific parameters, different for each RRL, are

$\mu$ = distance modulus

$E(B - V)$ = reddening

$a$ = amplitude

$\omega$ = frequency(1/P)

$\phi$ = phase.

For one RRL, the time series photometry can be written as $\{\{t_{bi}, m_{bi}, \sigma_{bi}\}_{i=1}^{n_b}\}_{b=1}^B$, where $m_{bi}$ is the observed magnitude at time $t_{bi}$ in filter b measured with (known) uncertainty $\sigma_{bi}$. The bands are indexed 1, ..., B instead of typical letters, e.g., ugriz. The model and data are related by

$$m_{bi} = m_b(t_{bi}) + \epsilon_{bi},$$

where $\epsilon_{bi} \sim N(0, \sigma_{bi}^2)$, meaning the noise parameter $\epsilon_{bi}$ can be viewed as a random normal variable with mean = 0 and $Var(\epsilon_{bi}) = \sigma_{bi}^2$.

This model assumes all RRL share a common shape by band and that RRL are strictly singly periodic functions. These assumptions are an approximation. For example, our model does not account for the amplitude and phase modulations, which vary according to an additional period caused by the Blažko effect (Blažko 1907). Rather than construct a perfectly accurate model, the goal is to construct a model with few free parameters that provides a better approximation to RRL variation than existing methods. For example, a simple sinusoid model fit to five filters has a total of 16 free parameters (five means, five amplitudes, five phases, and one frequency) while providing only a very rough approximation to the steep rise and slow decline in brightness observed in RRL light curves. In contrast, this model provides a significantly better approximation while fitting for five free parameters (or four if light curves are corrected for extinction prior to fitting).

## Appendix B

### Determining Template Population Values

We estimated the population parameters common to all RRL ($M_b(P)$, $R_b$, and $\gamma_b$) using a combination of theory and existing data sets.

For the filter-dependent extinction coefficients $R_b$, we assumed a Galactic reddening value of $R = 3.1$ from Fitzpatrick (1999). These extinction values $R_b$ for both SDSS and DES filters are summarized in Table 9. Note that these values are only used if the templates are fit with light curves uncorrected for extinction. In general, it is better to fit with dust-corrected light curves because the model has one fewer free parameter and the uncertainty on distance is greatly reduced. In this work, we corrected the light curves for extinction prior to fitting the template (see Section 3.2).

To develop the P-L relation for this work, we determined the values of $M_b(P)$ for ugriz using version 3.2 of the BASTI synthetic horizontal branch generator,<fnref n="58" /> based on the evolutionary tracks of Pietrinferni et al. (2004, 2006). We generated synthetic absolute magnitudes for RRL spanning -0.48 ≤ log P ≤ 0.08 with a metallicity of [Fe/H] = -2.85 to use as our starting values for $M_b(P)$. Then, we shifted the template

<footnote n="58">58 Available at albione.oa-teramo.inaf.it/BASTI/WEB_TOOLS/HB_SYNT/.</footnote>
<table><caption>**Table 9**
Extinction Coefficients</caption><thead><tr><th rowspan="2">Band<br>b</th><th colspan="2">R<sub>b</sub></th></tr><tr><th>SDSS</th><th>DES</th></tr></thead><tbody><tr><td>u</td><td>4.799</td><td>N/A</td></tr><tr><td>g</td><td>3.737</td><td>3.665</td></tr><tr><td>r</td><td>2.587</td><td>2.464</td></tr><tr><td>i</td><td>1.922</td><td>1.804</td></tr><tr><td>z</td><td>1.430</td><td>1.380</td></tr><tr><td>Y</td><td>N/A</td><td>1.196</td></tr></tbody></table>

**Note.** Based on Fitzpatrick (1999). Note that these model dust coefficients differ from those listed in Schlegel et al. (1998) and Abbott et al. (2018).

<table><caption>**Table 10**
P-L Coefficients</caption><thead><tr><th>Band b</th><th>c<sub>SDSS</sub></th><th>p<sub>1,SDSS</sub></th><th>p<sub>2,SDSS</sub></th><th>c<sub>DES</sub></th><th>p<sub>1,DES</sub></th><th>p<sub>2,DES</sub></th></tr></thead><tbody><tr><td>u</td><td>1.889</td><td>-0.049</td><td>-0.319</td><td>...</td><td>...</td><td>...</td></tr><tr><td>g</td><td>0.767</td><td>0.167</td><td>-0.595</td><td>0.730</td><td>-0.020</td><td>-0.065</td></tr><tr><td>r</td><td>0.550</td><td>-0.637</td><td>-0.353</td><td>0.542</td><td>-0.739</td><td>0.997</td></tr><tr><td>i</td><td>0.505</td><td>-1.065</td><td>-0.202</td><td>0.522</td><td>-1.136</td><td>-0.057</td></tr><tr><td>z</td><td>0.510</td><td>-1.308</td><td>-0.231</td><td>0.520</td><td>-1.292</td><td>-0.535</td></tr><tr><td>Y</td><td>...</td><td>...</td><td>...</td><td>0.558</td><td>-1.392</td><td>0.657</td></tr></tbody></table>

curve in each filter to match the magnitude offsets shown in real SDSS and DES light curves. We parameterized these empirical $M_b(P)$ by using a quadratic period–absolute magnitude (P–L) relation at a fixed metallicity of [Fe/H] ≈ −1.5 (see Section 5.4 for an extended discussion of the template metallicity):

$$M_b(P) = c_{0b} + p_{1b}(\log_{10}(P) + 0.2) + p_{2b}(\log_{10}(P) + 0.2)^2, \quad (8)$$

where the c, p₁, and p₂ values for the SDSS and DES filters are shown in Table 10.

The light curve shape in filter b $γ_b$ is a function of phase that covers one pulsation period. We use RRL found by Sesar et al. (2010) to estimate $γ_b$ for the SDSS filters. We assume the same shapes for DES griz and assume the DES Y band shape is the same as the z-band shape.

To infer the shape, we first “fold” the well-sampled S82 RRL light curves from Sesar et al. (2010) into phase coordinates by taking the modulus of the Modified Julian Dates of the observations with respect to the pulsation period of each individual RRL. We “phase-align” the light curves by shifting them so that they all reach their maximum brightness at phase = 0. We smooth the light curves by removing observations with photometric errors ≥0.2 and linearly interpolating them in equally spaced phase bins. Then, we shift all of the light curves so that the curves in each filter ugriz have an average magnitude value $m_b = 0$. We sample each light curve on a grid in phase space so that a single RRL is denoted by $X_{tb}$ for t = (1, ..., T) and b = (1, ..., B), where t indexes the phase (the grid has T = 100 equally spaced phases) and b indexes the filter (total B filters).

Let $X_{itb}$ be the magnitude for the ith RRL, at phase t in band b. Let $γ_b ∈ ℝ^T$ be the template in filter b. The template matrix of all five bands is thus defined as $Γ = (γ_1, ..., γ_B) ∈ ℝ^{T×B}$. Let $a ∈ ℝ^n$ be the amplitudes for the n RRL in the SDSS S82

sample. Let $X_{i·} ∈ ℝ^{T×B}$ represent the phase-folded, shifted, and normalized photometry described in the previous paragraph for the ith RRL. To determine the Γ matrix of the template shapes, we solve the following optimization problem:

$$\min_{\Gamma, a} \sum_{i=1}^{n} \|X_{i·} - a_i \Gamma\|_{F}^{2}, \quad (9)$$

where $||a||_2 = 1$ (the Euclidean norm) for identifiability and F denotes the Frobenius norm, or the square root of the absolute squares of its elements. The resulting Γ matrix is the template shapes in each filter. To reflect the dependence of the RRL amplitudes on the filter in which they were observed, we rescale the templates so the peak-to-peak g-band amplitude is 1 and the amplitudes of the other filter shapes are fractions of the g-band amplitude. These shapes and population parameter values form the set of templates that we use for fitting in our analysis.

## Appendix C
Fitting the Model

In this section, we describe how to fit the model to the data. Directly using the inverse of the observation uncertainties as weights is known to be suboptimal when the templates are an approximation; see Long (2017). We estimate a model error term $σ_{me}$, which is then used in the least-squares fitting. To compute $σ_{me}$, we fit the template to all well-sampled SDSS RRL light curves and compute the difference between the squared residuals and the squared photometric error $σ_{bi}^2$. $σ_{me}$ is the square root of the average of these differences. The value across all of the SDSS bands ugriz is 0.0547.

The model is fit by minimizing a weighted sum of squares (“$\chi^2$ minimization”). There are at most five free parameters: $μ$, $E[B-V]$, $a$, $ω$, and $φ$. The dust can be turned off in the fitting, in which case $E(B-V)$ is set to 0. We perform a grid search across the frequency because the objective function is highly multimodal. At frequency $ω$ in the grid, we solve for the four parameters $μ$, $E(B-V)$, $a$, and $φ$ using

$$\min_{\mu, E(B-V), a, \phi} \sum_{b,i} \frac{(m_{bi} - M_b(\omega) - \mu - E(B-V)R_b - a\gamma_b(\omega t_{bi} + \phi))^2}{\sigma_{bi}^2 + \sigma_{me}^2}, \quad (10)$$

where b is the filter index and i is the epoch index. We use a block-relaxation method in which we alternate between minimizing across the ($μ$, $E(B-V)$, $a$) parameters and minimizing across the $φ$ parameter. The number of iterations can also be specified.

When minimizing across ($μ$, $E(B-V)$, $a$) at fixed $φ$, the model is linear in ($μ$, $E[B-V]$, $a$), so we find the closed-form weighted least-squares solution. Occasionally, the update will result in a negative amplitude. In this case, we do a random phase update in the next step (i.e., draw phase uniformly in [0, 1]), rather than the Gauss-Newton method described below.

When minimizing across $φ$, with fixed ($μ$, $E[B-V]$, $a$), we cannot analytically solve for $φ$. Instead we use a Gauss-Newton method. Define

$$m_{bi}^* \equiv m_{bi} - M_b(\omega) - \mu - E(B-V)R_b \quad (11)$$
and

$$ \gamma_{bi}(\phi) = a\gamma_b(\omega t_{bi} + \phi). \qquad (12) $$

Then the objective function which we seek to minimize is

$$ g(\phi) = \sum_{b,i} \frac{(m_{bi}^* - \gamma_{bi}(\phi))^2}{\sigma_{bi}^2 + \sigma_{me}^2}. \qquad (13) $$

With $\phi^{(m)}$ as our current phase estimate, the Newton update has the form

$$ \phi^{(m+1)} = \phi^{(m)} - H(g)^{-1}(\phi^{(m)})\nabla(g)(\phi^{(m)}), \qquad (14) $$

where $\nabla(g)$ and $H(g)$ are the first and second derivatives of $g$. We have

$$ \nabla(g) = \frac{\partial g}{\partial \phi} = -2 \sum_{b,i} \frac{(m_{bi}^* - \gamma_{bi}(\phi))\gamma_{bi}'(\phi)}{\sigma_{bi}^2 + \sigma_{me}^2} \qquad (15) $$

and

$$ H(g) = \frac{\partial^2 g}{\partial \phi^2} = 2 \sum_{b,i} \frac{(\gamma_{bi}'(\phi))^2 - (m_{ib}^* - \gamma_{bi}(\phi))\gamma_{bi}''(\phi)}{\sigma_{bi}^2 + \sigma_{me}^2}. \qquad (16) $$

The Gauss–Newton update approximates $H(g)$ with

$$ H^*(g) = 2 \sum_{b,i} \gamma_{bi}'(\phi)^2, \qquad (17) $$

where we substitute $H^*$ for $H$ in Equation (14), rather than using $H(g)$ in Equation (16). This is a standard approach in nonlinear regression, which avoids computing $\gamma_{bi}''$ and ensures that the second derivative is positive (see Section 14.4 in Lange 2010). We approximate $\gamma_{bi}'$ by storing numerical derivatives of the $\gamma_b$ templates. At each new $\omega$ in the grid of frequency, we obtain a warm start for the $\mu$, $E[B-V]$, $a$, and $\phi$ parameters by using estimates from the last frequency. We choose the frequency at which the RSS is minimized to be the parameters of the best-fitting template to the data.

## Appendix D
### Template Code Products

The RRab template for both SDSS and DES filters and the fitting algorithm presented in this work are available at https://github.com/longjp/rr-templates and archived in Zenodo (doi:10.5281/zenodo.1252333). The template is originally implemented in R, but can be accessed in Python via the rpy2 module as was done in this analysis. Examples of the

template usage are available in both R and Python 3, though it is also compatible with Python 2.

When using the template-fitting functions, there are two options available to the user that impact the values returned by the fits. As described in the previous section, the template model includes an optional dust term. The model already includes extinction coefficients appropriate for both the DES and SDSS filter systems and can estimate the amount of extinction affecting the light curve as one of the parameters. However, if the light curves to be fit are already corrected for dust extinction prior to fitting, this parameter can be turned off to reduce the number of estimated parameters from five to four and thus improve the quality of the fits. This option is included in the code to allow the user to choose the option most appropriate for their data.

The other option determines whether or not to use the uncertainties associated with the individual observations in the light curve when performing the fits. In this analysis, we rescaled the uncertainties and used them when fitting the template to our data. However, if one suspects that the magnitude uncertainties in the light curves are misestimated or they simply are not available, this option can be turned off and the uncertainties will not be used. We leave this option open to the user.

The repository contains examples explaining how to fit the template to both DES and SDSS light curves using all combinations of these options. To make this code accessible to a variety of users, these examples are included as R and Python Jupyter notebooks.

## Appendix E
### Data Products

To enable further work with similar data, we provide all of the RRab candidate and training light curves at https://des.ncsa.illinois.edu/releases/other/y3-rrl and an extended version of Table 6. These light curves have already had their photometric uncertainties rescaled as described in Section 3.1 and include both the dust-corrected and uncorrected magnitudes in each band. The light curves are indexed by their DES Y3Q2 QUICK_OBJECT_ID numbers, which are included as a column in the full data table.

A description of the columns in the included data table is shown in Table 11, and example sample selection criteria are shown in Table 12. All of this information is also available in the documentation at https://des.ncsa.illinois.edu/releases/other/y3-rrl.
<table><caption>**Table 11**
Description of Data Columns</caption><thead><tr><th>Column Name</th><th>Description</th></tr></thead><tbody><tr><td>QUICK_OBJECT_ID</td><td>DES Y3Q2 ID Number</td></tr><tr><td>COADD_OBJECT_ID</td><td>DES DR1 ID Number</td></tr><tr><td>R.A.</td><td>R.A. in degrees (J2000)</td></tr><tr><td>Decl.</td><td>decl. in degrees (J2000)</td></tr><tr><td>p_ab</td><td>Classifier RRab score</td></tr><tr><td>EBV</td><td>Extinction value from Schlegel et al. (1998)</td></tr><tr><td>⟨g⟩ (rizY)</td><td>Mean g(rizY) magnitude measured in light curves (not extinction corrected)</td></tr><tr><td>nobs_g(rizY)</td><td>Number of observations in DES g(rizY) in the object's light curve</td></tr><tr><td>nobs</td><td>Total number of observations in the final light curve</td></tr><tr><td>period_0(1,2)</td><td>Period of the first (second, third) best-fitting template (days)</td></tr><tr><td>sigma_dp_p</td><td>Uncertainty in the ΔP/P of the best-fit period</td></tr><tr><td>amp_0(1,2)</td><td>Amplitude of the first (second, third) best-fitting template (mag)</td></tr><tr><td>sigma_da</td><td>Uncertainty in the Δa of the best-fit amplitude</td></tr><tr><td>mu_0(1,2)</td><td>Distance modulus of the first (second, third) best-fitting template (mag)</td></tr><tr><td>phase_0(1,2)</td><td>Phase offset of the first (second, third) best-fitting template</td></tr><tr><td>rss_0(1,2)</td><td>Residual sum of squares (RSS) of the first (second, third) best-fit template</td></tr><tr><td>chi2_g(rizY)</td><td>Reduced chi-squared of the light curve from a constant value in DES g(rizY)</td></tr><tr><td>sig_g(rizY)</td><td>"Significance" of the light curve in DES g(rizY)</td></tr><tr><td>rss_dof_0(1,2)</td><td>RSS per degree of freedom of the first (second, third) best-fitting template</td></tr><tr><td>lchi_med</td><td>Median log(reduced chi-squared) across DES grizY</td></tr><tr><td>rss_lchi_med</td><td>(RSS/dof)/median log(reduced chi-squared) across DES grizY</td></tr><tr><td>amp_rss_0(1,2)</td><td>Amplitude/(RSS/dof) for the first (second, third) best-fitting template</td></tr><tr><td>f_dist1_0(1,2)</td><td>Distance of the first (second, third) best-fit amplitude/period from the Sesar et al. (2010) Oosterhoff I relation</td></tr><tr><td>f_dist2_0(1,2)</td><td>Distance of the first (second, third) best-fit amplitude/period from the Sesar et al. (2010) division between the Oosterhoff I and II groups</td></tr><tr><td>kappa_0(1,2)</td><td>von Mises-Fisher concentration parameter of phases for the first (second, third) best-fitting template</td></tr><tr><td>flag_minmax</td><td>Light-curve sampling flag, 0:≥2 obs at max and min, 1:,&lt;2 obs at min, 1:,&lt;2 obs at max, 3:,&lt;2 obs at max and min</td></tr><tr><td>GaiaDR2_ID</td><td>Gaia DR2 source_id for cross-matched RRab in Clementini et al. (2018)</td></tr><tr><td>PS1_RRab</td><td>Present in Sesar et al. (2017) RRab catalog? 0 = no, 1 = yes</td></tr><tr><td>CSDR2_ID</td><td>ID number from associated Catalina Surveys DR2 RRab catalogs</td></tr><tr><td>ATLAS_ID</td><td>ID from Heinze et al. (2018) ATLAS RRab catalog</td></tr><tr><td>SDSS_ID</td><td>SDSS DR7 ID</td></tr><tr><td>SDSS_class</td><td>Classification from S82 studies from Sesar et al. (2010) or Ivezić et al. (2007)</td></tr><tr><td>Simbad_ID</td><td>Object identifier in SIMBAD database (Wenger et al. 2000)</td></tr><tr><td>Simbad_class</td><td>Object classification in SIMBAD database (Wenger et al. 2000)</td></tr><tr><td>OGLE_ID</td><td>OGLE ID from LMC and SMC RRAb catalogs from Soszyński et al. (2016)</td></tr><tr><td>CEMV_ID</td><td>ID from Martínez-Vázquez et al. (2016)</td></tr><tr><td>GCVS5_ID</td><td>Cross-matched ID with the General Catalogue of Variable Stars v5.1 (Ram et al. 2017)</td></tr><tr><td>GCVS5_Type</td><td>Variable type for cross-matched object as listed in the General Catalogue of Variable Stars v5.1 (Ram et al. 2017)</td></tr><tr><td>hpx32</td><td>Object's Healpix (nside = 32) used to subdivide the light curves</td></tr><tr><td>comments</td><td>Comments from first author during visual validation followed by a reason code</td></tr><tr><td>train</td><td>Identifies objects used to train the random forest classifier. 0 = no, 1 = yes</td></tr><tr><td>rrab</td><td>Identified as a high-confidence RRab in this study. 0 = no, 1 = yes</td></tr><tr><td>filepath</td><td>Filepath to the object's light curve</td></tr></tbody></table>

Note. All of these features are included in the table included in the data products available at https://des.ncsa.illinois.edu/releases/other/y3-rrl.

For the "comments" column, the possible comments and their meanings are as follows:

confirmed: convincing RRab

no_?: rejected followed by a reason code

maybe_?: ambiguous candidate followed by a reason code

missing_image: missing coadd image in SkyViewer

Reason codes:

galaxy—DES image showed galaxy

bad_fit—poor template fit

n_points—few observations and/or poor phase coverage

crowding—object is close to another object

misc—miscellaneous reasons
<table><caption>**Table 12**
Example Selection Criteria</caption><thead><tr><th>To select:</th><th>Choose:</th></tr></thead><tbody><tr><td>Objects that were identified as RRab by the classifier</td><td>p_ab ≥ 0.35</td></tr><tr><td>Objects found in the Clementini et al. (2018) Gaia DR2 RRab catalog</td><td>GaiaDR2_ID > 0</td></tr><tr><td>Objects found in S82</td><td>SDSS_class != 0</td></tr><tr><td>Objects found in the Sesar et al. (2017) Pan-STARRS catalog</td><td>PS1_RRab > 0</td></tr><tr><td>Objects found in the Catalina Surveys RRab catalog</td><td>CSDR2_ID != 0</td></tr><tr><td>Objects found in the ATLAS RRab catalog</td><td>ATLAS_ID != 0</td></tr></tbody></table>

## Appendix F
Simulating DES RRab Light Curves

Because DES is deeper than most overlapping surveys and much deeper than our training set in S82, we use simulated light curves to estimate the recovery rates at fainter magnitudes. We create the simulated light curves as follows:

(1) We generate light curve shapes from all 379 smoothed RRab light curve templates from Sesar et al. (2010) with gatspy's RRLyraeGenerated function (VanderPlas & Ivezić 2015). These generated template light curves already include the measured period and amplitude from their real observed counterpart light curves.

(2) We shift the light curves into the DES filter system using the DES-SDSS filter transformation relations:<fnref n="59" />

$$g_{DES} = g_{SDSS} + 0.001 - 0.075(g - r)_{SDSS}$$

$$\text{rms} = 0.021 \text{ per star}$$

$$r_{DES} = r_{SDSS} - 0.009 - 0.069(g - r)_{SDSS}$$

$$\text{rms} = 0.021 \text{ per star}$$

$$i_{DES} = i_{SDSS} + 0.014 - 0.214(i - z)_{SDSS} - 0.096(i - z)^2_{SDSS}$$

$$\text{rms} = 0.023 \text{ per star}$$

$$Y_{DES} = z_{SDSS} + 0.045 - 0.306(i - z)_{SDSS}$$

$$\text{rms} = 0.030 \text{ per star}$$

$$z_{DES} = z_{SDSS} + 0.022 - 0.068(i - z)_{SDSS}$$

$$\text{rms} = 0.025 \text{ per star}$$

(3) We use the distance estimates from Sesar et al. (2010) in the distance modulus equation to shift the light curves to a distance of 10 pc so that the light curves reflect the absolute magnitudes of the stars.

(4) We then sample the photometric measurements in the light curve at phases corresponding to the real DES cadence. The cadence is randomly selected from 1808 distinct fields in the DES wide-field footprint with unique observation times. This results in a light curve that is sampled in the same manner as the light curve for a real object somewhere in the DES footprint.

(5) To simulate the effects of distance on our recovery, we shift the downsampled light curve magnitudes to the apparent magnitudes they would have at a specified distance within the range that DES could detect in the single-epoch images. Once the magnitudes are shifted, we remove any magnitudes that are fainter than the median magnitude depth for that filter in the DES single-epoch images: $g \sim 23.5$, $r \sim 23.3$, $i \sim 22.8$, $z \sim 22.1$,

$Y \sim 20.7$ (Abbott et al. 2018). Thus, the light curves reflect the magnitude limits of each band at fainter magnitudes.

(6) Last, we assign a photometric uncertainty to each observation in the light curves. Following a method similar to Medina et al. (2018), we calculate the standard deviation of error-rescaled light curves in the survey region as a function of their mean magnitudes in each band as shown in Figure 2. We apply a shift in magnitude to the simulated observations using a Gaussian distribution with the appropriate standard deviation for that magnitude and band.

Following this procedure, we created 5685 simulated RRab light curves. A sample of these simulated light curves is shown in Figure 21. Because the photometric uncertainties were sampled from the rescaled uncertainties we applied to the real data, there was no need to rescale the errors using the method described in Section 3.1. Because we did not also simulate nonvariable light curves to analyze alongside the simulated RRab light curves, rescaling the errors using the same procedure would have removed real variable objects. We fit the template to all 4751 simulated light curves that passed the initial variability cuts and had at least five observations, and used the results to determine our detection efficiency. Results of this analysis are detailed in Section 5.3.

## ORCID iDs

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

J. P. Long https://orcid.org/0000-0001-5853-5938

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

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

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

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

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

E. Morganson https://orcid.org/0000-0001-7180-109X

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

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

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

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

E. Rykoff https://orcid.org/0000-0001-9376-3135

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

L. Strigari https://orcid.org/0000-0001-5672-6079

F. Paz-Chinchón https://orcid.org/0000-0003-1339-2683

D. Hatt https://orcid.org/0000-0003-2767-2379

A. H. Riley https://orcid.org/0000-0001-5805-5766

E. O. Nadler https://orcid.org/0000-0002-1182-3825

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

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

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

C. B. D'Andrea https://orcid.org/0000-0002-8198-0332

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

D. Gruen https://orcid.org/0000-0003-3270-7644

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

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

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

P. Melchior https://orcid.org/0000-0002-8873-5065

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

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

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

M. Smith https://orcid.org/0000-0002-3532-6970

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

<footnote n="59">⁵⁹ http://www.ctio.noao.edu/noao/node/5828#transformations</footnote>
F. Sobreira https://orcid.org/0000-0002-7822-0658

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

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

B. Yanny https://orcid.org/0000-0002-9541-2678

References

Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18

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

Baker, M., & Willman, B. 2015, AJ, 150, 160

Balbinot, E., Yanny, B., Li, T. S., et al. 2016, ApJ, 820, 58

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

Belokurov, V., Evans, N. W., Bell, E. F., et al. 2007, ApJL, 657, L89

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2008, ApJL, 686, L83

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2009, MNRAS, 397, 1748

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2010, ApJL, 712, L103

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

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006b, ApJL, 647, L111

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007, ApJ, 654, 897

Bergemann, M., Sesar, B., Cohen, J. G., et al. 2018, Natur, 555, 334

Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2014, MNRAS, 443, L84

Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2016, MNRAS, 463, 1759

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

Boettcher, E., Willman, B., Fadely, R., et al. 2013, AJ, 146, 94

Breiman, L. 2001, Machine Learning, 45, 5

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

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

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

Carrasco Kind, M., Drlica-Wagner, A., Koziol, A. M. G., & Petravick, D. 2019, JOSS, 4, 1022

Catelan, M. 2009, Ap&SS, 320, 261

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

Catelan, M., & Smith, H. A. 2015, Pulsating Stars (New York: Wiley-VCH)

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Clementini, G., Gratton, R., Bragaglia, A., et al. 2003, AJ, 125, 1309

Clementini, G., Greco, C., Held, E. V., et al. 2006, MmSAI, 77, 249

Clementini, G., Ripepi, V., Molinaro, R., et al. 2018, arXiv:1805.02079

Coupon, J., Kilbinger, M., McCracken, H. J., et al. 2012, A&A, 542, A5

Dark Energy Survey Collaboration, Abbott, T., Abdalla, F. B., et al. 2016, MNRAS, 460, 1270

Deason, A. J., Belokurov, V., Hamren, K. M., et al. 2014, MNRAS, 444, 3975

Desai, S., Armstrong, R., Mohr, J. J., et al. 2012, ApJ, 757, 89

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

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., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870

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

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

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

Fisher, R. 1953, RSPSA, 217, 295

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

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

Förster, F., Maureira, J. C., & San Martín, J. 2016, ApJ, 832, 155

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gilbank, D. G., Gladders, M. D., Yee, H. K. C., & Hsieh, B. C. 2011, AJ, 141, 94

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

Gratton, R. G., Bragaglia, A., Carretta, E., et al. 2003, A&A, 408, 529

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

Grillmair, C. J., & Carlin, J. L. 2016, ASSL, 420, 87

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

Heinze, A. N., Tonry, J. L., Denneau, L., et al. 2018, AJ, 156, 241

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

High, F. W., Stubbs, C. W., Rest, A., Stalder, B., & Challis, P. 2009, AJ, 138, 110

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

Huijse, P., Estévez, P. A., Förster, F., et al. 2018, APJS, 236, 12

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

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

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Ivezić, Ž., Lupton, R., Schlegel, D., et al. 2004a, in ASP Conf. Ser. 327, Satellites and Tidal Streams, ed. F. Prada, D. Martinez Delgado, & T. J. Mahoney (San Fransisco, CA: ASP), 104

Ivezić, Ž., Lupton, R. H., Schlegel, D., et al. 2004b, AN, 325, 583

Ivezić, Ž., Sesar, B., Jurić, M., et al. 2008, ApJ, 684, 287

Ivezić, Ž., Smith, J. A., Miknaitis, G., et al. 2007, AJ, 134, 973

Jupp, P. E., & Mardia, K. V. 1989, AnSta, 7, 3

Kaluzny, J., Stanek, K. Z., Krockenberger, M., et al. 1998, AJ, 115, 1016

Kelly, P. L., von der Linden, A., Applegate, D. E., et al. 2014, MNRAS, 439, 28

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

King, D. S., & Cox, J. P. 1968, PASP, 80, 365

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

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

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

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

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

Lange, K. 2010, Numerical Analysis for Statisticians (Berlin: Springer)

Li, T. S., Balbinot, E., Mondrik, N., et al. 2016, ApJ, 817, 135

Long, J. P. 2017, EJS, 11, 1

LSST Science Collaboration, Marshall, P., Anguita, T., et al. 2017, arXiv:1708. 04058

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

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

MacDonald, E. C., Allen, P., Dalton, G., et al. 2004, MNRAS, 352, 1255

Majewski, S. R., Ostheimer, J. C., Rocha-Pinto, H. J., et al. 2004, ApJ, 615, 738

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., Stetson, P. B., Monelli, M., et al. 2016, MNRAS, 462, 4349

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

Mateu, C., Vivas, A. K., Downes, J. J., et al. 2012, MNRAS, 427, 3374

Mateu, C., Vivas, A. K., Zinn, R., Miller, L. R., & Abad, C. 2009, AJ, 137, 4412

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

Medina, G. E., Muñoz, R. R., Vivas, A. K., et al. 2017, ApJ, 845, L10

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

Mondrik, N., Long, J. P., & Marshall, J. L. 2015, ApJL, 811, L34

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

Neilsen, E., & Annis, J. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Fransisco, CA), 77

Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245

Odenkirchen, M., Grebel, E. K., Rockosi, C. M., et al. 2001, ApJL, 548, L165

Oluseyi, H. M., Becker, A. C., Culliton, C., et al. 2012, AJ, 144, 9

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

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

Pickering, E. C., Colson, H. R., Fleming, W. P., & Wells, L. D. 1901, ApJ, 13, 226P

Pieres, A., Santiago, B. X., Drlica-Wagner, A., et al. 2017, MNRAS, 468, 1349
Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168

Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2006, ApJ, 642, 797
Prudil, Z., Grebel, E. K., Dékány, I., & Smolec, R. 2018, MNRAS, 480, 669

Richards, J. W., Starr, D. L., Butler, N. R., et al. 2011, ApJ, 733, 10

Rocha-Pinto, H. J., Majewski, S. R., Skrutskie, M. F., Crane, J. D., & Patterson, R. J. 2004, ApJ, 615, 732

Saha, A., & Vivas, A. K. 2017, AJ, 154, 231

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

Samus', N. N., Kazarovets, E. V., Durlevich, O. V., Kireeva, N. N., & Pastukhova, E. N. 2017, ARep, 61, 80

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

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

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

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

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

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

Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, AJ, 134, 2236

Sharma, S., Johnston, K. V., Majewski, S. R., et al. 2010, ApJ, 722, 750
Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114

Soszyński, Udalski, A., Szymański, M. K., et al. 2016, ApJ, 862, 114

Sra, S. 2016, arXiv:1605.00316

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

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

Süveges, M., Sesar, B., Váradi, M., et al. 2012, MNRAS, 424, 2528

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

The Dark Energy Survey Collaboration 2005, arXiv:astro-ph/0510346

Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, PASP, 130, 064505

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

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2018, arXiv:1811.04082

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

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

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

van Albada, T. S., & Baker, N. 1973, ApJ, 185, 477

VanderPlas, J. T., & Ivezić, Ž. 2015, ApJ, 812, 18

Vivas, A. K., Olsen, K., Blum, R., et al. 2016, AJ, 151, 118

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

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, ApJ, 554, L33

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

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

Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9

Willman, B. 2010, AdAst, 2010, 285454

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

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

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579

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

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