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

DES (DES Collaboration 2005; DES Collaboration et al. 2016) was a six-year optical/near-infrared imaging survey covering ∼5000 deg² of the southern Galactic cap using the Dark Energy Camera (DECam; Flaugher et al. 2015) mounted at the prime focus of the 4 m Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO). Observations were completed in 2019 January. DES obtained ∼10 × 90 s exposures in five broadband filters, grizY (Neilsen et al. 2019). ⁵¹ DES observed with gri in dark time, iz in gray time, and zY in bright time, with each field of the footprint being observed two to three times per year (Diehl et al. 2016, 2019; Neilsen et al. 2019). The 5σ point-source depth of the DES exposures is estimated to be grizY ∼ (24.3, 24.1, 23.5, 22.9, 21.5) (Morganson et al. 2018).

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

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

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

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

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

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

$\begin{eqnarray}&&{\tilde{\chi }}_{\nu ,b}^{2}=\displaystyle \frac{1}{{N}_{b}-1}\displaystyle \sum _{1}^{{N}_{b}}\displaystyle \frac{{\left({m}_{i,b}-\mathrm{med}({m}_{b})\right)}^{2}}{{\sigma }_{i,b}^{2}},\end{eqnarray} \tag{ 1 }$

where m_i,b is the observed PSF magnitude of an object in observation i, σ_i,b is the reported magnitude uncertainty on that observation, and N_b is the total number of observations of that object in that band. We fit a quadratic function of the form

$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}({\tilde{\chi }}_{\nu ,b}^{2}) & = & \ {c}_{0,b}\\ & & +{c}_{1,b}(\mathrm{med}({m}_{b})-20)\\ & & +{c}_{2,b}{\left(\mathrm{med}({m}_{b})-20\right)}^{2}.\end{array}\end{eqnarray} \tag{ 2 }$

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

Zoom In Zoom Out Reset image size

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

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

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

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

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

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

Zoom In Zoom Out Reset image size

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

Since numerous color, magnitude, and variability measurements are calculated (many of which are correlated), we use the random forest algorithm to "learn" the optimal boundaries in feature space to separate RRab from non-RRab. Random forests use a collection of decision trees to predict an object's type. Each decision tree is trained on a random subset of the training population by repeatedly subdividing the sample based on feature values until a user-defined maximum depth or other specified stopping condition is reached. The specific features that are used in each split are chosen randomly and can influence the characteristics of the population that a tree learns to detect. A random forest classifier takes advantage of the learning differences between individual trees by averaging the predictions from a large number of trees to produce an aggregate score (Amit & Geman 1997; Breiman 2001). Random forests are a good choice for this task because they are largely insensitive to uninformative features, so the inclusion of a feature that does not separate the object types well will not harm the overall results. For our classifiers, we use the RandomForest implementation in scikit-learn (Pedregosa et al. 2012). To reduce our sample of light curves down to a manageable size for template fitting, we divide our pretemplate selection criteria into two phases: (1) remove nonvariable sources, and (2) remove common variable contaminants (i.e., QSOs).

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

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

$\begin{eqnarray}&&{F}_{\beta }=\displaystyle \frac{{\left(1+\beta \right)}^{2}\cdot \mathrm{TP}}{(1+{\beta }^{2})\cdot \mathrm{TP}+{\beta }^{2}\cdot \mathrm{FN}+\mathrm{FP}},\end{eqnarray} \tag{ 3 }$

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

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

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

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

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

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

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

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

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

$\begin{eqnarray}\begin{array}{rcl}\mathrm{MCC} & = & \displaystyle \frac{\mathrm{TP}\times \mathrm{TN}-\mathrm{FP}\times \mathrm{FN}}{\sqrt{(\mathrm{TP}+\mathrm{FP})(\mathrm{TP}+\mathrm{FN})(\mathrm{TN}+\mathrm{FP})(\mathrm{TN}+\mathrm{FN})}},\end{array}\end{eqnarray} \tag{ 4 }$

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

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

Zoom In Zoom Out Reset image size

Due to the limited depth and area of external RRL catalogs, we are required to use simulated RRL to assess the completeness of our selection for faint RRL (D ≳ 100 kpc). To assess the recovery of our template-fitting and classification process and to assess how frequently RRc are misclassified as RRab, we follow the procedure described in Section 3.2 to simulate an independent set of light curves for ∼3100 RRab and ∼900 RRc. We subject these simulated light curves to the same minimum number of observations, color cuts, variability classifier, template fitting, and final classifier that are applied to the real data. The recovery rate and precision of the parameter estimates are shown in Figure 6.

Zoom In Zoom Out Reset image size

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

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

The variable star content of the Magellanic Clouds (MCs) has been extensively studied, with the OGLE-IV catalog of Soszyński et al. (2019) containing 47,828 RRL in both MCs. Although the DES footprint only covers the outskirts of the MCs, a comparison with the OGLE catalog still yields a relatively large number of objects in common (792) and allows us to assess the recovery of light-curve properties of the DES Y6 RRab candidates. The OGLE catalog is ideal for this comparison because its light curves have <100 epochs (and consequently, period estimates and classifications are very robust) and it contains many types of variables (allowing us to investigate the contamination rate of our catalog).

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

Zoom In Zoom Out Reset image size

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

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

We search for new RRab in the MCs using the DES Y6 coverage beyond the OGLE footprint. We start by selecting a range of distance moduli for each MC based on stars in common with OGLE. For the LMC, the distance modulus distribution showed a clear peak at 17.9 < μ < 18.8, containing 608 RRL. The SMC overlap with OGLE is significantly smaller, but there is still a clear peak in the distance modulus distribution at 18.1 < μ < 19.2, with 27 RRL matching OGLE in that distance range. Matched OGLE-DES stars are shown as black diamonds in Figure 8. We then select DES Y6 RRab candidates that were not in OGLE but have distance moduli consistent with the selection defined above and have an angular separation of less than 25° and 15° from the centers of the LMC and SMC, respectively. These stars are potential new members of the MCs. Most of them are located outside the OGLE footprint, although there are also some stars that may have been missed by OGLE.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

The old stellar population of the LMC extends to large angular separations, as shown in previous studies. Saha et al. (2010) traced the LMC population with main-sequence turnoff stars out to 16° from its center. Using a similar technique, Nidever et al. (2019) further extended the detection to 21°. In addition, Belokurov & Koposov (2016) used publicly available DES images to find a lumpy distribution of blue horizontal branch (BHB) stars extending out to 30° and possibly up to 50°. We note in passing that some of our candidates match the location and distance of their "S1" substructure. Our strong detection of RRL in the LMC periphery confirms these previous results and strengthens the case for an extended, old (≳10 Gyr) stellar population in the LMC.

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

An interesting feature associated with the SMC in this part of the sky is the so-called Small Magellanic Cloud Northern Over-Density (SMCNOD; Pieres et al. 2017), whose area is indicated by a blue ellipse in the left panel of Figure 8. Pieres et al. (2017) estimated that the SMCNOD is primarily composed of an intermediate-age population of ∼6 Gyr, which is too young to contain RRL. Indeed, the variable star content in the SMCNOD was examined by Prudil et al. (2018) with an earlier release of the OGLE catalog (Soszyński et al. 2017). Although they found eight OGLE RRab inside the SMCNOD area, they concluded their density was compatible with the expected SMC background at that distance. In contrast, we find almost twice as many (15 RRab candidates) in the same region.

In addition to the MCs, two other prominent overdensities of RRab candidates are within the DES footprint: the classical dwarf spheroidal (dSph) galaxies Sculptor (α, D = 1502, 86 kpc) and Fornax (α, D = 3996, 147 kpc). Both can be clearly seen in Figure 4.

Martínez-Vázquez et al. (2016a, hereafter MV16) presented the most complete and extensive study of the variable star population in Sculptor, reporting 536 RRL, of which 289 were RRab. We search for potential Sculptor members by selecting all RRL in our catalog out to 1.5× its tidal radius of 691 (Munoz et al. 2018) and with distance moduli within 19.67 ± 0.75. We find 116 RRL within these limits, all but four matching the catalog of MV16. The lower-than-average completeness in this region (39% at μ ∼ 19.5) is due to crowding effects near the center of the galaxy.

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

Three of the four new candidate members associated with this galaxy lie outside the footprint explored by MV16. One is located at 93', significantly beyond the tidal radius of Sculptor but with a compatible Gaia DR2 proper motion. The mean distance modulus of our Sculptor RRL is 〈μ〉 = 19.49 ± 0.09 mag, while Martínez-Vázquez et al. (2015) found a value of 19.62 ± 0.04 mag. Although small, the difference can be attributed to the fact that our analysis assumes all RRL have the mean metallicity of the Milky Way halo, while the old stellar population in Sculptor has a wide range of metallicities extending down to [Fe/H] ∼ −2.4 (Martínez-Vázquez et al. 2016b).

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

Surprisingly, 51 of our Fornax RRL candidates are located beyond its tidal radius. These stars are uniformly distributed around the galaxy (Figure 10), reaching out to 114' from its center (note that our search radius only extended to 1162). In a recent study based on data from DES Y3, Wang et al. (2019) concluded that no significant extratidal disturbances are observed down to a surface brightness limit of ∼32.1 mag arcsec⁻². Fornax is known to have several bursts of star formation and is dominated by a population of age ∼5 Gyr (de Boer et al. 2012; Rusakov et al. 2021). While this dominant population is too young to have produced the observed RRL, an old stellar population (>10 Gyr) is also present in Fornax. This older stellar population has been found to be more spatially extended than the younger populations (Wang et al. 2019). The detection of RRL outside the tidal radius suggests a very low surface brightness, old extratidal population.

Zoom In Zoom Out Reset image size

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Large-scale anisotropies in the halo RRL distribution have been claimed by Iorio et al. (2018) using a combined catalog of Gaia+2MASS RRL. Boubert et al. (2019) provided supporting evidence for this structure using observations from the Catalina Surveys (e.g., Torrealba et al. 2015; Drake et al. 2017) and the sample of RR Lyrae variables identified in PS1 (Hernitschek et al. 2016). While this structure is significantly closer than the excess observed here (Galactocentric distance of ∼20 kpc), Boubert et al. (2019) claim a Magellanic origin for this overdensity, which could extend to greater distances where the infall track of the Magellanic Clouds crosses the DES footprint. Further investigation of these distant candidates is necessary to better understand possible anisotropies in the RRL distribution at distances ≳80 kpc.

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

$\begin{eqnarray}\rho (x,y,z)={\rho }_{0}\left\{\begin{array}{cc}{r}_{e}^{{n}_{1}}, & {r}_{e}\lt {R}_{0}\\ {r}_{e}^{{n}_{2}}\cdot {R}_{0}^{{n}_{1}-{n}_{2}}, & {r}_{e}\geqslant {R}_{0},\end{array}\right.\end{eqnarray} \tag{ 5 }$

where ρ₀ is the density normalization, n₁ is the inner power-law index, n₂ is the outer power-law index, and R₀ is the break radius (in elliptical Galactocentric coordinates). We perform a binned Poisson maximum-likelihood fit for the observed counts in each bin, k, as a function of the model parameters, θ = {ρ₀, R₀, n₁, n₂}. Our likelihood analysis is described in more detail in Appendix C. Within each bin, we correct the observed number of RRL for the detection completeness of our catalog estimated from simulated RRab described in Section 5.1. We fit the model parameters using a Markov Chain Monte Carlo ensemble sampler (emcee; Foreman-Mackey et al. 2013), and we report the median, 16th, and 84th percentiles of the marginalized posterior distributions of each parameter in Table 2. We perform this analysis for both the full DES Y6 RRab candidate sample and the high-confidence RRab associated with the Oosterhoff I sequence (Figure 12).

Table 2. Broken-power-law Halo Density Parameters

Sample	ρ0	R0	n1	n2
	(kpc−3)	(kpc)
Full	${990}_{-230}^{+300}$	${32.1}_{-0.9}^{+1.1}$	$-{2.54}_{-0.09}^{+0.09}$	$-{5.42}_{-0.14}^{+0.13}$
Oo-I	${240}_{-80}^{+110}$	${31.3}_{-1.0}^{+1.3}$	$-{2.32}_{-0.13}^{+0.14}$	$-{5.59}_{-0.19}^{+0.17}$

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

Download table as:  ASCIITypeset image

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

Zoom In Zoom Out Reset image size

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

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

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

We calculate the significance of a putative satellite at each search location (i.e., at the location of each RRab candidate in our catalog) as the Poisson probability of detecting k or more RRL given an expectation of λ:

$\begin{eqnarray}&&P(k| \lambda )=\displaystyle \frac{{\lambda }^{k}{e}^{-\lambda }}{k!}.\end{eqnarray} \tag{ 6 }$

To select candidate satellites, we apply a significance threshold of p < 3 × 10⁻⁷, which corresponds to a one-sided Gaussian significance of 5σ. We also require that satellite candidates consist of at least three RRab candidates.

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

$\begin{eqnarray}&&\mathrm{log}({N}_{\mathrm{RRL}})=-0.29\times {M}_{V}-0.80.\end{eqnarray} \tag{ 7 }$

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

Zoom In Zoom Out Reset image size

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

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

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

A second overdensity (SubId = 2) consists of three RRab candidates located at (α,δ) ∼ 349.8, −55.6 at a distance of 59.5 kpc. This overdensity is in the same region of the sky and at roughly the same distance as the Tucana II dwarf galaxy, D = 57.5 kpc; however, it is separated from Tucana II by ∼5 deg on the sky. Interestingly, the Gaia DR2 proper motions of the RRab candidates in this overdensity are similar to the proper motions of confirmed members of Tucana II, albeit with large uncertainties (Figure 15). It has been suggested that the diffuse structure of Tucana II could be an indication of tidal disruption (e.g., Bechtol et al. 2015). Belokurov & Koposov (2016) showed evidence for extended stellar structure around Tucana II using BHB stars (i.e., their "S2a" cloud at μ ∼ 18.8), and Chiti et al. (2020) used Gaia proper motions combined with photometric metallicities to identify probable members >1 deg from Tucana II. Further observations are required to confirm an association between this overdensity of RRab candidates and the Tucana II dwarf galaxy.

Zoom In Zoom Out Reset image size

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