The Next Generation Virgo Cluster Survey. XXXVII.
Distant RR Lyrae Stars and the Milky Way Stellar Halo out to 300 kpc

The Next Generation Virgo Cluster Survey (NGVS) is a deep optical imaging survey of the Virgo cluster of galaxies, carried out with the 1 deg² MegaCam instrument on the 3.6-m Canada-France-Hawaii Telescope (CFHT) in the $u^{*}$, $g^{\prime}$, $i^{\prime}$, and $z^{\prime}$ bands (with limited additional coverage in the $r^{\prime}$ band—not used in this work). The survey consists of 117 pointings (not including 4 background pointings), with slight overlap between adjacent pointings. The resulting NGVS footprint covers a contiguous area of 104 deg² of the nearby Virgo cluster, out to the virial radii for both the Virgo A and B subclusters. The NGVS project motivations, strategy, and observational program are discussed by Ferrarese et al. (2012). The main goal of the NGVS is to study the galaxy population of the Virgo cluster. The observing strategy adopted, however, was atypical of deep, extragalactic optical surveys, and involved imaging multiple fields in sequence before returning to the same field (see § 2 of Ferrarese et al., 2012, for a full description). The purpose of this strategy was to improve sky subtraction on large scales, and it naturally resulted in an observing cadence where individual observations of a single field were spaced by a minimum of $\sim\!1$ hr and by as much as a few yr. This range of temporal spacing makes the NGVS data well-suited to finding both short- and long-term variables.

The NGVS images reach point source depths of [$u^{*}$ ($S/N=5$), $g^{\prime}$ ($S/N=10$), $i^{\prime}$ ($S/N=10$), $z^{\prime}$ ($S/N=5$)]$\;\approx\;$[26.3, 25.9, 25.1, 24.8] mag for the stacked images, and [24.5, 24.4, 23.6, 23.1] mag for single exposures. This paper is based on all NGVS single exposure images obtained during the period March 2009 through May 2013. The total number of epochs across the four bands is roughly the same ($\sim\!38$) for each MegaCam pointing (tile in the survey mosaic), with the exception of small overlap regions between adjacent pointings which have twice as many epochs or more. The single-exposure integration times [and approximate number of exposures per band] are 582 s [$\times 13$] in $u^{*}$, 634 s [$\times 7$] in $g^{\prime}$, 411 s [$\times 8$] in $i^{\prime}$, and 550 s [$\times 10$] in $z^{\prime}$. The sky was clear during most exposures, with a representative seeing FWHM $\sim 0\farcs 8$, although the image quality in $i^{\prime}$ was always better than $0\farcs 6$. The basic parameters of the NGVS dataset are summarized in Table 1. The NGVS dataset represents the deepest uniform imaging available across the entire Virgo cluster, and is likely to remain so for some time.

Figure 1 shows how the NGVS compares to recent and future surveys of RR Lyrae stars in terms of survey area, number of epochs, and depth. The NGVS occupies an interesting part of this survey parameter space: its area is comparable to Sloan Digital Sky Survey (SDSS) Stripe 82 and HiTS, its number of epochs is comparable to PS1 and DES Y6, and its single-epoch depth is $\sim\!1$ mag fainter than DES. Only with Rubin/LSST will there be a survey of the distant MW halo with better single epoch depth than NGVS and with a substantially larger area.

The pre-processed individual NGVS images were stacked for each of the four bands, and the stacked images in the four bands were re-gridded onto a common astrometric system. These steps were carried out using the MegaPipe software (see Gwyn, 2008; Ferrarese et al., 2012). The SExtractor software was used in “dual image mode” to detect sources in the $g^{\prime}$ band and carry out forced photometry in all four bands (of the four, $g^{\prime}$ is the one in which faint RR Lyrae have the highest signal-to-noise ratio in the stacked NGVS images). Aperture photometry was carried out for all detected objects on the stacked NGVS $u^{*}$, $g^{\prime}$, $i^{\prime}$, and $z^{\prime}$ images using the following set of aperture diameters: 3, 4, 5, 6, 7, 8, and 16 pixels, and a local sky annulus of $(r_{\text{inner}},~{}r_{\text{outer}})=(16,~{}20)$ pixels. where each pixel corresponds to 0$\farcs$187.

A three-step calibration/correction process was applied in order to obtain optimal point-source photometry, while accounting for the fact that the point spread function (PSF) varies from field to field across the NGVS footprint and from band to band. First, a sample of bright (but unsaturated in NGVS) stars were selected in each field and their 16-pixel (3$\farcs$0) aperture instrumental magnitudes were photometrically calibrated to their SDSS PSF magnitudes. Second, a map of aperture correction as a function of sky position was created based on the curves of growth of a larger sample of bright stars (extending to fainter magnitudes than the SDSS calibration stars) in each of the four bands. Third, the aperture magnitudes (diameters of 3, 4, 5, 6, 7, and 8 pixels) of all objects were corrected to 16-pixel aperture total magnitudes based on the aperture correction map. The 16-pixel aperture contains practically all of the flux of a point source, given the sub-arcsecond to arcsecond (FWHM) seeing conditions under which the NGVS data were obtained (Table 1). In this way, we have compiled a homogeneous point-source photometry database that is calibrated against SDSS. To summarize, $m_{b,N}$ for each point source represents its total magnitude in band $b$, aperture-corrected from its $N$-pixel aperture magnitude on the stacked image. In this study, we use the aperture-corrected 4-pixel ($0\farcs 75$) diameter total apparent magnitudes in the four bands, $u^{*}_{4}$, $g^{\prime}_{4}$, $i^{\prime}_{4}$, and $z^{\prime}_{4}$ ; they are hereafter referred to simply as $u^{*}$, $g^{\prime}$, $i^{\prime}$, and $z^{\prime}$ magnitudes, respectively.¹¹1The stacked-image photometry used through most of this paper is the flux-averaged apparent brightness over the NGVS cadence. This cadence is sparse and varies from band to band for any given source. For our NGVS RR Lyrae candidates, we compute the complete flux-averaged apparent magnitude in each band $b$, $\langle{m_{b}}\rangle$, by integrating over the full light curve based on the best-fit pulsational parameters (see §§ 3–5).

At the faint apparent magnitudes that this study explores, any survey of stars must contend with contamination by background galaxies and quasi-stellar objects (QSOs). In order to distinguish between point sources and extended sources, we adopt a fuzziness index: $\Delta i^{\prime}=i^{\prime}_{4}-i^{\prime}_{8}$. An object is considered to be a point source if $-0.05\leq\Delta i^{\prime}\leq+0.05$. With the above definition of $m_{b,N}$, we expect the distribution of $m_{b,N_{1}}-m_{b,N_{2}}$ for point sources to concentrate around 0. Of the four NGVS bands, we chose to use $i^{\prime}$ magnitudes because: (1) the $i^{\prime}$ images were obtained under the best seeing conditions ($\rm FWHM\sim 0$$\farcs$6); and (2) the field-to-field variance in the seeing FWHM was smallest in the $i^{\prime}$ band. Our point source selection criteria are similar to those used in previous NGVS studies (Durrell et al., 2014; Liu et al., 2015). The distribution of $\Delta i^{\prime}$ versus $g^{\prime}$ magnitude shown in Figure 2 displays a clear vertical locus of point sources (mostly stars, along with some QSOs) around $\Delta i^{\prime}=0$. The rest of this paper is limited to NGVS point sources as defined by the above $\Delta i^{\prime}$ criteria.

To search for RR Lyrae in the NGVS database, we first select all point sources in the region of $(u^{*}-g^{\prime})$ versus $(g^{\prime}-i^{\prime})$ and $(g^{\prime}-i^{\prime})$ versus $(i^{\prime}-z^{\prime})$ color-color space that is occupied by known RR Lyrae—i.e., those that are located to the lower left of the blue dashed lines in both color-color diagrams (Figure 3). Since these colors are based on stacked-image photometry, they can be quite different from the true colors for variable objects, because the NGVS cadence is sparse and different from band to band. In other words, it is not surprising that known RR Lyrae display a relatively large scatter in our color-color diagrams. The dashed lines we use take this scatter into account.

Of the 94 RR Lyrae discovered by Sesar et al. (2017) in the NGVS footprint using PS1 survey data: 84 are confirmed as RR Lyrae in this work (see § 3); 7 are bright stars ($g^{\prime}\sim 18.0$–18.5) for which a portion of the NGVS time-series photometry in one or more of the bands suffers from saturation; and 3 are fainter objects ($g\sim 20.5$) that are classified as QSOs or non-variable in our work. Our final color-color selection criteria, shown as blue dashed lines in Figure 3, were designed to include all 84 of these reconfirmed PS1 RR Lyrae whose photometry is free of saturation in NGVS.

We also defined regions around the stellar locus in the two color-color diagrams, as indicated by the pink shaded regions in Figure 3, to select bright ($18<g^{\prime}<21$), non-variable stars. This yields 54,043 point sources ($\sim 520$ per NGVS pointing) that are used as photometric reference stars for the calibration of our time-series photometry, as described in § 2.4 below. Of these, 402 sources (0.7%) turned out to be variable candidates (see § 2.5); the variable fraction is so low that it has a negligible effect on the calibration of the NGVS time-series photometry.

Our study is the first to explore time-domain photometry in the NGVS database. As a result, we had to develop specific data analysis methods to account for exposure-to-exposure temporal variations in the seeing FWHM in each band and residuals in the atmospheric transparency correction.

To estimate the systematic uncertainty of our photometry, and to ensure our measurements across different epochs are self-consistent, we did an error-rescaling by calibrating the reduced chi-squared statistic, $\chi_{\nu,\,b}^{2}$. This was calculated using the median magnitude $\overline{m_{b}}$ in a given band $b$ for each light curve:

in which $N_{b}$ is the number of epochs of star $\nu$ in band $b$, $m_{i,\,b}$ is the $i^{\rm th}$ measurement in that band, and $\sigma_{i,\,b}$ is the photometric uncertainty which is a combination of random (Poisson) error and systematic error:

We assume that the systematic error $\sigma_{\rm sys}$ should consist of two terms: a constant term $\sigma_{0,\,b}$ that accounts for bright nearby objects and/or bad sky subtraction, and a linear term $\sigma_{1,\,b}\times 10^{\frac{\langle m_{\nu,\,b}\rangle}{2.5}}$ that scales with the apparent brightness of the star and accounts for effects like CCD non-homogeneity. For non-variable sources, which are the majority of our sample, it is reasonable to expect that $\chi_{\nu,\,b}^{2}\sim 1$. Therefore, the parameter values $(\sigma_{0,\,b},\,\sigma_{1,\,b})$ that quantify the systematic error in each band $b$ are optimized to ensure that the $\log\,\chi_{\nu,\,b}^{2}$ versus $m_{b}$ trend line is horizontal and centered on zero, as shown in the lower panel of Figure 4. This parameter optimization involves minimizing the sum of sigma-clipped $\left|\log{\chi_{\nu,\,b}^{2}}\right|$ for all NGVS stars in each band $b$. The photometric uncertainties reach $\sim 0.05$ mag at $g^{\prime}\sim 24$ (upper panel of Figure 4). Compared to the two predecessor surveys, our NGVS single-epoch photometry is 1.7 mag deeper than DES and 2.3 mag deeper than HiTS at an error level of 0.05 mag.

After the time-series photometry is calibrated, we select variable candidates that satisfy all three of the following criteria: (1) $\sigma_{\nu,\,b}\geq\sigma_{\text{med},\,b}$; (2) $\ln{\chi^{2}_{\nu,\,b}}\geq 2$; and (3) $N_{\nu,\,b}\geq 3$, where the meanings of subscripts are the same as in § 2.4. In other words, if an object $\nu$ has elevated photometric RMS and $\chi^{2}$ in all four bands relative to the corresponding median values for objects of comparable apparent magnitude, and has at least three photometric measurements in each band, it is considered a variable candidate. In total, 1685 variable candidates passed the variability cuts and are used as the input for our template fitting algorithm.

To identify RR Lyrae and derive robust estimates of their light curve parameters like period and amplitude, we performed template fitting based on the empirical RR Lyrae light curves generated from high cadence observations of 483 RR Lyrae in the SDSS Stripe 82 (Sesar et al., 2010). Of these 483 templates, 379 are of type RRab and 104 are of type RRc. Our light curve model has the form:

in which $a_{b}$, $m_{0,\,b}$ are free parameters (amplitude and magnitude at peak brightness, respectively), $T_{b,\,k}$ is the $k$th normalized light curve in band $b=\left\{u^{*},\,g^{\prime},\,i^{\prime},\,z^{\prime}\right\}$, and phase $\phi$ is determined by:

Note the $\phi_{0}$ in our work was calculated by using the MJD system for the time parameter $t$, and the value of $\phi$ was restricted to $0\leq\phi\leq 1$. To get the best-fit values of the free parameters $\left\{a_{b},\;m_{0,\,b},\;P,\;\phi_{0}\right\}$, we minimize the following $\chi^{2}$ value:

in which $t_{b,\,n}$, $m_{b,\,n}$ and $\sigma_{b,n}$ are the time, apparent magnitude, and photometric uncertainty of the $n$th observation in band $b=\left\{u^{*},\,g^{\prime},\,i^{\prime},\,z^{\prime}\right\}$. In the end, a periodogram score, $s_{\text{fitting}}=1-\frac{\chi^{2}}{\chi_{0}^{2}}$, is estimated for the best-fit period to quantify the goodness of fitting. This fitting process is very similar to that which is described in Sesar et al. (2017), but the difference is that we applied the normalized light curve model and did not fix the ratio of the light curve amplitudes and colors between the $u^{*}$ band and the other three bands. The periodogram is calculated with a brute-force search in period space from 0.2 to 0.9 d, with a step size of 0.5 s. For each given period, the fitting of pulsational parameters was done with the Python routine scipy.optimization.minimize. Here is our rationale for our adopted template-fitting scheme:

1. 1.

   We decided not to apply the color and amplitude restrictions from the SDSS RR Lyrae templates to our fitting, because the transformations from the SDSS to the CFHT MegaCam photometric system is not well-calibrated for horizontal branch stars like RR Lyrae, especially at the blue end.

2. 2.

   The grid interval used in our initial fitting is estimated as

   $$\Delta t\sim f\frac{P^{2}}{T_{\rm obs}}$$

   (7)

   where $T_{\rm obs}$ is the whole timespan of the NGVS observation, $P$ is an estimated typical period of RR Lyrae and $f$ is an arbitrarily chosen fudge factor to define the desired phase accuracy. The above relationship is derived based on a series of heuristic assumptions. We assume the maximum number of pulsation cycles of RR Lyrae to be $\sim\frac{T_{\rm obs}}{P}$, and the single-cycle pulsation phase shift resulting from the period search grid size to be $\frac{\Delta t}{P}$. The cumulative phase error would therefore be $\frac{T_{\rm obs}\Delta t}{P^{2}}$, and with a phase accuracy level $f$ selected (we used $f=0.05$ in this work), we can derive the above equation to estimate the necessary grid interval for the period search.

We understand that the periodogram score defined above is more about the goodness of the fitting and therefore an imperfect indicator of whether an object is an RR Lyrae. We also understand that the chi-square minimization algorithm we have used does not guarantee that the global minimum has been found, mainly because of the extra degrees of freedom we allowed for the $u^{*}$ band data. The template fitting scheme we have used, however, is a practical choice given our available computational resources, the number of the NGVS photometry epochs, and our multi-band light curve models. Our initial fitting is at least good enough to track the peaks where the global optimal solution resides, and based on the definition of $\Delta t$ shown above, it is reasonable to expect that the periodogram should be mostly smooth within the $\Delta t$ scale.

After the initial light curve fitting process was completed, our sample of variable candidates was cross-matched against the known QSOs and AGNs classified with SDSS and XMM-Newton data (Zhang et al., 2021). We found 131 matches with spectroscopically confirmed galaxies in SDSS, photometric redshift based galaxies in SDSS, and X-ray luminous QSOs and active galactic nuclei (AGNs) in XMM-Newton. Presumably, the variability we are detecting is related to nuclear activity in these background galaxies. It is reassuring that, for these objects, our light curve fitting algorithm returns low fitting scores: $s_{\rm fitting}<0.3$.

We visually vetted the light curves of the remaining 1554 objects, and classified them into four groups; (a) 366 probable RR Lyrae candidates, with a relatively high fitting score in the initial round ($s_{\text{fitting}}>0.5$), good phase coverage, and obvious short-term variability in its unfolded light curve; (b) 71 marginal variable sources; (c) 615 suspected QSO/AGN contaminants with a low fitting score ($s_{\text{fitting}}<0.3$), that display clear long-term variability through multiple observing seasons with similar photometric trends across the different bands; and (d) 502 sources that appear to be non-variable that are likely to be affected by systematic errors (e.g., saturation, cosmic rays, charge bleeds from neighboring bright stars, effect of variable seeing on close neighbors, detector artifacts, and other possible calibration errors in our time-series photometry). Representative examples of these four categories are shown in Figure 7. All objects classified as type (a) and (b) are further analyzed with a more detailed light curve fitting process.

We perform a zoom-in fitting to all variable candidates around the 10 periodogram peaks shown in their initial fitting results, in order to avoid missing the global minimum of the periodogram as the ideal $\Delta t$ interval in Equation 7 may vary for RR Lyrae with different periods and observation baselines. For each peak $P_{\text{peak}}$ in the initial periodogram, we search its nearby $\pm 0.5$ s range with a zoom-in fitting grid interval $\Delta t_{\text{zoom}}=0.02$ s. The period value that corresponds to the highest fitting score in the zoom-in fitting process is then recorded as the best-fit period $P_{\text{best}}$. In most cases, the $P_{\text{best}}$ value occurs around the tallest peak in the initial fitting, as shown in the upper panels of Figure 5. However, we also found cases where $P_{\text{best}}$ occurs around suboptimal peaks in the initial fitting, as shown in the lower panels of Figure 5, suggesting that we were slightly overestimating the ideal $\Delta t$ interval in the initial fitting. The best-fit period $P_{\text{best}}$ and its corresponding fitting parameters (best-fit templates, pulsational parameters, and fitting score) were then used for further vetting process. Each variable candidate is assigned the type (RRab or RRc) of the best-fit SDSS Stripe 82 RR Lrae template.

The result of applying this fitting procedure to the unfolded light curves of the 84 PS1 RR Lyrae candidates is illustrated in Figure 6. Our multi-band template fitting method accurately measures periods for $91\%$ of RR Lyrae stars ($97\%$ of RRab and $70\%$ of RRc stars). The period is recovered to within 1 sec for $73\%$ of RR Lyrae stars. For all fittings with score $s_{\text{fitting}}>0.93$, the period differences are within $3$ min. Note that if the period fitting returns a discrepant value, this can predominately be attributed to one-day beat frequency aliasing.²²2The beat frequency is: $f_{\text{beat}}=f_{\text{true}}+Nf_{\text{sample}}\;(N=\pm 1,\pm 2,...)$. For ground based surveys like NGVS, $f_{\text{sample}}\approx 1$ d^-1. The one-day beat period (in d) is therefore given by: $P_{\text{beat}}=\frac{P_{\text{true}}}{1+NP_{\text{true}}}$, as shown by the curved dashed lines in the upper panel of Figure 6 (for $N=\pm 1$).

The differences in the fitted pulsational parameters are remarkably small considering that we are using completely different observational data, as well as sparsely sampled observations compared with the PS1 study.

We also note that we recover the halo RR Lyrae star discovered by Ciardullo et al. (1989) near M49 (NGC 4472), which was the most distant MW halo star observed at that time. The period that we fit is within 15 s (0.03%) of that measured by the original discoverers.

We visually inspected the stacked and single exposure $g^{\prime}$ band images of all variable candidates with $s_{\text{fitting}}>0.9$ to guard against a couple of factors could lead to false positives in our RR Lyrae search. First, the NGVS footprint contains a large number of star-forming Virgo cluster galaxies whose photometrically variable supergiants could masquerade as MW halo RR Lyrae (e.g., the study of NGC 4535 by Spetsieri et al., 2018). Second, the combination of variable seeing and sky subtraction errors for point sources that are located in complicated backgrounds (e.g., within the dust lanes/spiral arms of star-forming galaxies in the Virgo cluster) can compromise the fidelity of the stacked-image and time-series photometry. We found 8 objects with complicated backgrounds and noticed that their scores are in the range: $0.9<s_{\text{fitting}}<0.93$. This led us to use an RR Lyrae selection criterion of $s_{\text{fitting}}>0.93$.

We searched for our faintest NGVS RR Lyrae candidates ($g^{\prime}>21.5$) in the SMOKA archive to check if they had been observed with the Subaru 8-m telescope, and found 11 matched stars. The single exposure Hyper Suprime-Cam HSC-$g$ band images of these objects, obtained between 2014 and 2018, show that none of them are transients. A similar determination was made for 2 additional NGVS RR Lyrae candidates that were matched against Hubble Space Telescope Advanced Camera for Surveys images in the MAST archive.

In total, we have found 180 RR Lyrae candidates in the NGVS field, of which there are 139 RRab and 41 RRc. The light curve fittings of all RR Lyrae candidates passed the threshold periodogram score value of $s_{\text{fitting}}>0.93$ and were checked visually. The best-fit pulsation parameters and other information are available in a machine-readable table, with the column definition and the first few table entries shown in Table 2. For each NGVS RR Lyrae candidate, we integrate over its best-fit light curve—based on its best-fit pulsational parameters: period, amplitude, phase, and light curve shape—to calculate its flux-averaged apparent magnitude $\langle{m_{b}}\rangle$ in each band (where $b=u^{*}$, $g^{\prime}$, $i^{\prime}$, and $z^{\prime}$). The extinction was corrected using the dust map in Chiang (2023) and the $u^{*}g^{\prime}i^{\prime}z^{\prime}$ extinction coefficients calibrated by Muñoz et al. (2014).

All of our RR Lyrae candidates are projected more than 5 $r_{\text{eff}}$ from the center of the closest (in projection) bright Virgo cluster galaxy, with two exceptions: (1) one is projected within 1 $r_{\text{eff}}$ of the center of the bright and smooth giant elliptical galaxy M49 in the Virgo cluster, and is an RR Lyrae that was independently discovered by Ciardullo et al. (1989); and (2) the other is projected 4.8 $r_{\text{eff}}$ from the center of the early-type spiral galaxy NGC 4492 in the Virgo cluster, and is an RRab that was independently discovered by PS1 (Sesar et al., 2017).

Although all of our RR Lyrae candidates have been visually inspected, it would nevertheless be useful to obtain additional photometry to confirm the RR Lyrae classification of those candidates that are poorly sampled in one or more bands. In Figure 8, we show a sky map of our NGVS RR Lyrae candidates.

We directly adopted the $i^{\prime}$-band $P$-$L$ relation in Sesar et al. (2017), which was calibrated with PS1 RR Lyrae data, since there are no known MW star clusters or satellite dwarf galaxies in the NGVS footprint for an independent calibration. We calculated the distances of our RR Lyrae candidates with their best-fit $i^{\prime}$-band flux-averaged magnitude and the $P$-$L$ relation in Sesar et al. (2017), where the absolute $i^{\prime}$-band magnitude is estimated as:

and for RRc stars their periods are fundamentalized before calculating their absolute magnitudes (Catelan, 2009):

The uncertainty of the distance estimation calculated by the $i^{\prime}$-band $P$-$L$ relation in Sesar et al. (2017) was reported to be $\sigma_{M_{i^{\prime}}}=0.06~{}(\text{random})\pm 0.03~{}(\text{systematic})$. Here, we validate that the uncertainty of the distance estimation in our study can be constrained from the PS1 result. First, the differences between the flux-averaged $i$-band magnitudes in our study and in PS1 for overlapping RR Lyrae stars are at a remarkably low level ($<0.01$ mag for non-aliasing cases). Second, the error induced by the “fundamentalization” of RRc star periods (by Catelan 2009) is smaller than 0.002 mag. Third, the $P$-$L$ relation has only a weak dependence on metallicity ($Z$). Sesar et al. claimed that their period-luminosity-metallicity ($P$-$L$-$Z$) relation has a scatter of $\sigma_{\rm DM}=0.09$ with a reference metallicity of $\rm[Fe/H]=-1.5$. While faint stars in the Virgo direction reside in the outer halo and could be metal-poor, the dependence of absolute magnitude on metallicity is weak. Even for the extreme metal-poor case, the resulting error in the distance modulus is $\sigma\sim 0.1$ mag, corresponding to a $5\%$ error in distance even at 300 kpc.

We estimate the completeness of our NGVS RR Lyrae sample by using a large synthetic RR Lyrae data set that mimics the cadence and photometric errors of the NGVS. In order to generate the light curve of a synthetic RR Lyrae with a flux-averaged magnitude $\langle g^{\prime}_{\text{synth}}\rangle$, the set of $ugiz$ light curve templates associated with a random ‘parent star’ with flux-averaged magnitude $\langle g_{\text{parent}}\rangle$ is selected from the Sesar et al. (2010) catalog of 483 SDSS RR Lyrae. The pulsational parameters of the parent star are adopted: period and $ugiz$ amplitudes. The difference $\Delta\langle g^{\prime}_{\text{synth}}\rangle=\langle g^{\prime}_{\text{synth}}\rangle-\langle g_{\text{parent}}\rangle$, which is equal to the difference in distance moduli between the synthetic and parent RR Lyrae, is then added to the light curves of the parent star in all four bands, and a random initial phase is assigned. The mock NGVS observation of this synthetic RR Lyrae is simulated by randomly selecting a point source in our NGVS time-series catalog, importing its epochs in the four bands, and calculating the folded phases and magnitudes of each simulated NGVS observation based on the light curve of the synthetic star. Next, this time-series data set is converted from the SDSS $ugiz$ system to the NGVS $u^{*}g^{\prime}i^{\prime}z^{\prime}$ system by applying the inverse of the photometric calibration procedure mentioned in § 2.2 (for details, see Gwyn, 2008). Finally, we incorporate realistic photometric errors in order to create a synthetic RR Lyrae time-series data set.

In total, we generated 7200 synthetic RR Lyrae (4800 RRab and 2400 RRc) time-series data sets that sample the full range of periods, amplitudes, and light curve shapes of the 483 SDSS Stripe 82 RR Lyrae, but shifted to the apparent magnitude range $20.5\leq\langle g^{\prime}\rangle\leq 24.5$. This mock RR Lyrae sample was analysed using the same light curve template fitting process described in §§ 3.1–3.4. A synthetic RR Lyrae is considered recovered if it satisfies the multi band variability criteria and fitting score threshold, and the best-fit period is within $\pm 5$ min of the period of the parent star.

The resulting completeness results are shown in Figure 11. For RRab stars, the completeness is in the range 85%–90% at bright magnitudes, and drops from 88% to 58% over the range $\langle g^{\prime}\rangle=23.4$–24.45. For RRc stars, the completeness is in the range 80%–85% at bright magnitudes, and drops from 80% to 36% over the range $\langle g^{\prime}\rangle=22.8$–24.4. The smaller pulsation amplitudes, shorter periods, and more symmetric light curve shapes of RRc relative to RRab makes it harder for the light curve fitting procedure to recover them. The apparent magnitude corresponding to the 85% recovery rate for RRab is 1.8 mag deeper for NGVS than HiTS (Medina et al., 2018) and 1.3 mag deeper for NGVS than DES (Stringer et al., 2021). The NGVS RRab completeness curve shown in Figure 11 is used to correct the raw radial density profile of MW halo RR Lyrae, as described in § 4.2. The RRab completeness level is assumed to be constant at 90% (average of the first three bins in Figure 11) from $\langle g^{\prime}\rangle=19.0$–20.5, a range over which there is no saturation and the photometric accuracy is high. It is worth mentioning that these completeness estimates do not take into account the areal completeness of the NGVS survey. However, given that the 5 $r_{\text{eff}}$ regions of brightest 350 galaxies combined cover only 1.1% of the NGVS survey footprint (shown in Figure 8), the effect of the Virgo cluster background on our completeness estimates is negligible.

In Figure 9, we demonstrate the robustness of our NGVS RR Lyrae sample by showing the folded and unfolded light curves for the 6 most distant (faintest) candidates. The quality of the template fits for these 6 distant RR Lyrae is comparable to that of brighter NGVS RR Lyrae (Figure 10) and better than that of comparably distant/faint HiTS and DES RR Lyrae candidates (Medina et al., 2018; Stringer et al., 2021, e.g., see their Figure 3). The properties of the full sample of 180 NGVS RR Lyrae are presented in Table 2.

Our NGVS sample of 180 RR Lyrae candidates roughly doubles the number of known RR Lyrae in this region of the sky: PS1 reported 94 RR Lyrae, whereas the additional RR Lyrae we have found in the NGVS database are, for the most part, fainter than the detection limit of PS1 (§ 2.2). The primary advantage of our NGVS data set over HiTS (Medina et al., 2018) and DES (Stringer et al., 2021) is significantly greater single-epoch photometric precision/depth (Figure 4). As a result, we expect that our sample of NGVS RR Lyrae candidates represents a more robust detection and reliable characterization of the most distant known RR Lyrae in the MW halo. Our NGVS sample of RR Lyrae candidates includes 39 with $d_{\rm helio}>100$ kpc and 7 with $d_{\rm helio}>200$ kpc. We show the full sample of RR Lyrae stars as a function of distance and RA in Figure 12; the six most distant objects whose light curves are shown in Figure 9 are marked by the bold black star symbols.

We construct the MW stellar halo number density radial profile $\rho({R_{\rm GC}})$ in this direction of the sky based on our NGVS RR Lyrae sample, where $R_{\rm GC}$ is the distance to the center of the Galaxy. For this number density calculation, we consider only the 140 RRab candidates identified in our study, spread over an area of 104 deg² in the direction of the Virgo cluster. The Galactocentric radius $R_{\rm GC}$ is calculated as follows:

where $d_{\rm helio}$ is the heliocentric distance, $l$ and $b$ are the Galactic latitude and longitude of each star, respectively, and $R_{\odot}$ is the distance from the Sun to the Galactic center. We adopt an $R_{\odot}$ value of 7.9 kpc (VERA Collaboration et al., 2020). We fit the radial number density profile of our RRab sample with a single power-law model: $\rho=\rho_{0}(R_{\rm GC}/R_{0})^{n}$ for a spherical halo, in which $\rho_{0}/R_{0}^{-n}$ is the normalization constant. The bins are evenly spaced on a log scale as shown in Figure 13, and the fitting was weighted with Poisson errors of the density values in each bin. We restrict the power-law fit to $R_{\rm GC}>45$ kpc because: (1) as these distances, the NGVS sample is unaffected by saturation; and (2) this avoids the well-known Sagittarius stream and Virgo Overdensity substructure along this line of sight (e.g., Donlon et al., 2019). Various observational studies have reported the existence of a break in the halo density radial profile at $R_{\rm GC}\sim 20$–35 kpc (e.g., Zinn et al., 2014; Xue et al., 2015) and have therefore adopted a broken power-law model to mark the transition between the inner (in-situ) halo and outer (ex-situ) halo. For this work, however, we only probe the outer halo beyond 45 kpc in the Virgo direction, so use only a single power-law model.

We detect the existence of RR Lyrae stars out to $R_{\rm GC}\sim 300$ kpc without any clear evidence of a break in the power-law density profile out to this radius. This corroborates similar results in the DES Y6 RR Lyrae catalog (Stringer et al., 2021), and the most distant Mira stars recently presented by Nikzat et al. (2022) with VVV survey data. Also, the Galactic splashback radius, namely the edge of the MW halo, is predicted to be $\sim 0.8r_{200\rm m}=290\pm 61$ kpc by Deason et al. (2020), which corresponds to the Galactocentric radii of the most distant RR Lyrae presented in our work.

The slope of the number density radial profile of the Galactic outer stellar halo has been studied using different stellar tracers, including RR Lyrae, K giants, blue horizontal branch/blue straggler stars, and A stars (see Hernitschek et al., 2018; Medina et al., 2018; Stringer et al., 2021). The derived density slopes of the outer halo lie in a wide range, between $n=-3.8\pm 0.1$ to $n=-5.4\pm 0.1$. Our result, $n=-4.09\pm 0.10$, is consistent with the $n=-4.17^{+0.18}_{-0.20}$ slope obtained by Medina et al. (2018) using RR Lyrae stars, and is within the range of other slope estimates. It should be emphasized that our sample of NGVS RR Lyrae and the DES Y6 RR Lyrae are currently the only two probes of the outer stellar halo out to $R_{\rm GC}\sim 300$ kpc. However, Stringer et al. (2021) opted to only fit the DES RR Lyrae stellar density radial profile in the distance range $30<d_{\rm helio}<100$ kpc, and obtained a steep outer halo slope of $n=-5.42\pm 0.13$. The density profile of their RR Lyrae candidates is significantly shallower than this power-law slope beyond 100 kpc (see their Figure 11). Their choice to limit the power-law fit to $d_{\rm helio}<100$ kpc was motivated by two factors: (1) they believe that there is unaccounted for QSO contamination in their sample of RR Lyrae candidates at large distances; and (2) any anisotropy in the MW halo caused by the infall of the Magellanic Clouds may invalidate the assumption of spherical symmetry. Our sample of NGVS RR Lyrae candidates, with its higher precision single-epoch photometry, multi-year time baseline, and visual vetting to guard against long period variables, is likely to be cleaner of QSOs, even at these large distances.

We compare our measured MW halo density profile slope with predictions from galaxy formation simulations. Pillepich et al. (2014) fit the logarithmic slope of the spherically-averaged stellar density profile out to the virial radius for $\sim 5000$ MW analogs in the Illustris simulation. Overall, they measured slopes that range from $-5.5<n<-3.5$ with a mean value of $n\sim-4.5$ by combining the results from light ($6\times 10^{11}M_{\odot}<M_{\rm halo}<9\times 10^{11}M_{\odot}$) and massive ($9\times 10^{11}M_{\odot}<M_{\rm halo}<2\times 10^{12}M_{\odot}$) dark matter halos, with more massive halos exhibiting flatter slopes. They also showed that halos that formed more recently (i.e., had a recent merger event) or accreted a larger fraction of their stellar components from satellite galaxies, exhibit flatter stellar halo slopes. Our NGVS slope result, $n=-4.09\pm 0.1$, which is slightly flatter than the average simulation slope, may suggest a more massive Galactic halo mass, or a relatively active recent accretion history of the MW’s outer halo.

A major shortcoming of our study is the small sky coverage of NGVS (104 deg²) relative to DES ($\sim 5,000$ deg²; Stringer et al., 2021) and PS1 ($\sim 30,000$ deg²; Sesar et al., 2017). The smaller the field of view, the greater the likelihood that the measured MW halo density profile could be biased by underlying substructure. With this caveat in mind, our NGVS MW halo RR Lyrae sample is consistent with the Sanderson et al. (2017) model predictions for MW analogs: $\sim 10$ field RR Lyrae at $\sim 300$ kpc over the whole sky, while we found one such star in 104 deg². Their Figure 7 shows significant anisotropy in the distribution of distant RR Lyrae, and their Figure 5 shows substantial variation in the number of such stars across their different halo simulations.

Oosterhoff (1939) discovered the bimodality of RRab stars (Oo I and Oo II) from different globular clusters (GCs) in period versus amplitude space (the “Bailey diagram”), and subsequent studies refined this empirical scenario by demonstrating that Oo I GCs are more metal-rich than Oo II GCs (e.g., Kinman, 1959). Moreover, the distribution of RR Lyrae on the Bailey diagram is independent of the uncertainties in their distance and reddening. The above factors make the Bailey diagram useful in investigating the formation environment of RR Lyrae. Metallicity differences among RR Lyrae stars from different progenitor dwarf galaxies should be imprinted in the Oosterhoff dichotomy in the Bailey diagram. A recent study by Fabrizio et al. (2019) compiled a catalog of all spectroscopically confirmed field RR Lyrae stars and analyzed their distribution in the Bailey diagram. Unlike RRab stars in GCs, field RRab stars have a continuous (unimodal, rather than bimodal) distribution in the Bailey diagram. While the period, amplitude, and metallicity of field RRab stars do not follow simple linear correlations, RRab stars trend towards being more metal-rich when one moves from long to short periods at fixed amplitude (see Figure 13 of Fabrizio et al., 2019).

In Figure 14, we show the distribution of our NGVS RRab stars in the Bailey diagram: amplitude $A_{g^{\prime}}$ (mag) versus the logarithm of the period log${}_{10}P$ (d). In Figure 15, we plot each RRab stars’ horizontal distance from the Oo I locus in the Bailey diagram $\Delta{\rm log}_{10}P$ (d), a rough proxy for metallicity, versus Galactocentric distance $R_{GC}$ (kpc). The analytic relation for the Oo I locus is a quadratic line fit by Sesar et al. (2010):

in which $A_{g^{\prime}}$ is the RRab’s $g^{\prime}$-band amplitude (mag) and $P$ is the pulsation period (d). The Oo II locus is fitted by offsetting the Oo I locus by +0.03 in the $\log{P}$ direction. This proxy for the Oo II locus is the same as the definition in § 4.6 of Sesar et al. (2010). The Oo II line is used to characterize the long-period field RRab subset that do not form a tight secondary sequence. We use the definition: $\Delta\log{P}=\log{P}-\log{P^{\prime}}$, where $P^{\prime}$ is the period interpolated from the Oo I locus at the same $A_{g^{\prime}}$.

Figure 15 indicates some degree of clustering around the Oo I locus for stars within $R_{\rm GC}\sim 100$ kpc, but the distribution appears to gradually shift to a broader distribution around the Oo II locus with increasing $R_{\rm GC}$. This shift in Figure 15, or equivalently the systematic shift from the shorter to longer periods in Figure 14, suggests a change from the metal-rich to the metal-poor regime for the RR Lyrae stars with increasing Galactocentric distance. Our estimates of the amplitude and period of RR Lyrae stars may be affected by Blazhko modulations that are undetectable given the sparse sampling of the NGVS RR Lyrae light curves, especially for distant stars with relatively low-accuracy photometry. This is a possible explanation for the broad distribution of our distant RR Lyrae around the Oo II sequence. On the other hand, Fabrizio et al. (2019), show that the distribution of field RR Lyrae in the Bailey diagram is intrinsically broad and continuous. We believe the variance of the periods of our distant stars reflect the metallicity variance of their progenitor host environments, which are mostly metal-poor for stars beyond 150 kpc. Also, the mean period of our most distant ($R_{\rm GC}>150$ kpc) RRab sample ($P_{\rm mean}=0.658\pm 0.03$ d) is broadly consistent with the recent census results of Gaia RR Lyrae stars in the nearby ultra faint dwarf (UFD) satellites ($P_{\rm mean}\sim 0.667$ d; see Vivas et al., 2020). This, combined with the gradient we see towards longer periods with increasing Galactocentric radius (Figure 15), suggests that metal-poor UFD satellites could be a main contributor to the stellar component of the outermost Galactic halo.