Taking the pulse of the outer Milky Way with HOWVAST: an RR Lyrae density proﬁle out to >200 kpc

In order to constrain the evolutionary history of the Milky Way, we hunt for faint RR Lyrae stars (RRLs) using Dark Energy Camera data from the High cadence Transient Survey (HiTS) and the Halo Outskirts With Variable Stars (HOWVAST) survey. We report the detection of ∼ 500 RRLs, including previously identiﬁed stars and ∼ 90 RRLs not yet reported. We identify 9 new RRLs beyond 100kpc from the Sun, most of which are classiﬁed as fundamental-mode pulsators. The periods and amplitudes of the distant RRLs do not place them in either one of the two classical Oosterhoﬀ groups, but in the Oosterhoﬀ intermediate region. We detect two groups of clumped distant RRLs with similar distances and equatorial coordinates, which we interpret as an indication of their association with undiscovered bound or unbound satellites. We study the halo density proﬁle using spheroidal and ellipsoidal ( 𝑞 = 0 . 7) models, following a Markov chain Monte Carlo methodology. For a spheroidal halo, our derived radial proﬁle is consistent with a broken power-law with a break at 18 . 1 + 2 . 1 − 1 . 1 kpc separating the inner and the outer halo, and an outer slope of − 4 . 47 + 0 . 11 − 0 . 18 . For an ellipsoidal halo, the break is located at 24 . 3 + 2 . 6 − 3 . 2 kpc and the outer slope is − 4 . 57 + 0 . 17 − 0 . 25 . The break in the density proﬁle is a feature visible in diﬀerent directions of the halo. The similarity of these radial distributions with previous values reported in the literature seems to depend on the regions of the sky surveyed (direction and total area) and halo tracer used. Our ﬁndings are compatible with simulations and observations that predict that the outer regions of Milky Way-like galaxies are mainly composed of accreted material.


INTRODUCTION
In the currently favoured cosmological framework, the Λ cold dark matter (Λ-CDM) model, galaxies assemble hierarchically through the accretion of smaller systems.The Milky Way (MW) and similar massive disc galaxies likely experienced numerous mergers in their early history as part of their hierarchical formation (see, e.g., Press & Schechter 1974;Blumenthal et al. 1984;Bullock & Johnston 2005;Montalbán et al. 2021).The stellar halos of these galaxies provide key information to help reconstruct their formation conditions.For the MW, in particular, compelling evidence for past and ongoing accretion events have been identified in present-day inner and outer halo stellar populations, unveiling details of gravitational interactions with massive satellites such as the Sagittarius stream (e.g., ★ E-mail: gustavo.medina@astro.utoronto.caIbata, Gilmore, & Irwin 1994;Majewski et al. 2003;Vivas & Zinn 2006), Gaia-Sausage-Enceladus (GSE; e.g., Belokurov et al. 2018b;Helmi et al. 2018;Haywood et al. 2018), and the infall of the Magellanic Clouds (e.g., Mathewson, Cleary, & Murray 1974;Besla et al. 2007;Zaritsky et al. 2020;Erkal et al. 2021).
The accretion history of a particular halo is also imprinted in the shape of its stellar density profile (e.g., Bullock & Johnston 2005;Cooper et al. 2013), as the mass distribution is sensitive to properties such as the halo formation time, the amount of stellar mass accreted, and how long ago the last mergers took place (see e.g.Pillepich et al. 2014).The slope of the number density profile of outer halo stars, in particular, has been shown to be a parameter of cosmological significance, closely related to the accretion his-cially close to the 'edge' of the MW (292±61 kpc, when defined as the point at which virialized material has completed at least two pericentric passages; Deason et al. 2020).Among the commonly used stellar tracers are the RR Lyrae stars (RRLs; see e.g.Drake et al. 2013b;Medina et al. 2018;Stringer et al. 2021;Huang & Koposov 2022).RR Lyrae variables are old (ages > 10 Gyr) and metal-poor ([Fe/H] typically < −0.5) stars of paramount importance for Galactic studies given their status of precise distance indicators (∼5 per cent from period-luminosity relations; Catelan & Smith 2015, and references therein), and their easy identification in time-domain surveys (given their intrinsically high luminosity and light curve shapes).The pulsation periods of RRLs typically range from 0.2 to 0.9 d (e.g., Smith 1995) and they are mainly classified into three subtypes, according to the nature of the pulsations, their periods, amplitudes, and light curve shapes.RRab stars are fundamental-mode pulsators with sawtooth light curves and typically longer periods, whereas RRc stars pulsate in the first overtone and have sinusoidal light curves with shorter periods.RRd stars are a less common subtype and pulsate in the fundamental mode and the first overtone simultaneously, with the first overtone being the principal mode of pulsation.
Given that the distances of RRLs are known with great precision, their spatial distribution can be derived and thus they can be used to study the radial density profile of the Galaxy (Wetterer & McGraw 1996;Vivas & Zinn 2006;Cohen et al. 2017;Iorio et al. 2018).This also makes them excellent tracers of the outermost limits of our Galaxy, as well-characterized stars at such large distances (beyond 100 kpc) are scarce (Sesar et al. 2017;Medina et al. 2018;Stringer et al. 2021), thus playing a key role in the estimation of the MW mass (see e.g., Eadie & Harris 2016;Deason, Belokurov, & Sanders 2019;Deason et al. 2021;Rodriguez Wimberly et al. 2022;Prudil et al. 2022).
A large number of RRL catalogues have been compiled over the years in existing large sky surveys, which cover a wide range of photometric depths (hence distances) and different regions of the sky.These systematic searches include the Quasar Equatorial Survey Team (QUEST) and the La-Silla QUEST surveys (Vivas et al. 2004;Zinn et al. 2014), the Northern Sky Variability Survey (NSVS; Kinemuchi et al. 2006), the Sloan Digital Sky Survey (SDSS) Stripe 82 (Sesar et al. 2010), the Catalina surveys (Abbas et al. 2014;Drake et al. 2014Drake et al. , 2017;;Torrealba et al. 2015), the Panoramic Survey Telescope And Rapid Response System survey (Pan-STARRS-1, or PS-1; Chambers et al. 2016;Hernitschek et al. 2016;Sesar et al. 2017), the Optical Gravitational Lensing Experiment (OGLE) survey (Soszyński et al. 2016(Soszyński et al. , 2019)), the second and third data releases of the Gaia mission (Holl et al. 2018, and Clementini et al. 2019and Clementini et al. 2023 using the Specific Objects Study pipeline, SOS), the Dark Energy Survey (DES; DES Collaboration 2016; Stringer et al. 2021), and the Zwicky Transient Facility survey (ZTF; Masci et al. 2019;Chen et al. 2020;Huang & Koposov 2022).Only a small subset of these surveys have allowed astronomers to reliably detect RRLs beyond 100 kpc, mostly due to instrumental limitations, while thousands of RRLs are predicted to be found in the halo between 100 and 300 kpc (Sanderson et al. 2017).In this vein, Medina et al. (2018) identified 16 RRL candidates with Galactocentric distances > 100 kpc over a ∼120 deg 2 area using data from the High cadence Transient Survey (HiTS; Förster et al. 2016).The current census of distant RRLs is still likely incomplete, and the outer limits of the Galaxy have yet to be comprehensively mapped.This serves as motivation for larger surveys focusing on the detection of RRLs at large distances.
In this study we introduce the Halo Outskirts With Variable Stars (HOWVAST) survey, with which we aim to extend the reach of known well-characterized outer halo RRL surveys out to ∼270 kpc.In Section 2, we describe the survey strategy, the observations carried out for this study, and the methodology followed for data processing.In Section 3, we provide a detailed description of our RRL selection and classification pipelines, as well as the methods used for the determination of their periods, amplitudes, and distances.Additionally, we contrast our detected RRLs with those from the literature and use these comparisons as an indicator of our detection completeness.In Section 4, we focus our attention on the most distant RRLs in our sample (those with H > 100 kpc).Finally, in Section 5 we study the spatial distribution of our RRLs following a Markov Chain Monte Carlo methodology, and discuss the similarities and differences between our results and studies of other regions of the halo.We conclude this manuscript by summarizing our results and outlining the implications of our findings in Section 6.

DECam data
The data used in this work were obtained as part of three independent campaigns carried out with the Dark Energy Camera (DE-Cam;Flaugher et al. 2015), which is mounted on the 4 m telescope at Cerro Tololo Inter-American Observatory (CTIO).The first campaign corresponds to the HiTS survey, which was originally designed to characterize the early stages of supernovae explosions in real time (Förster et al. 2016).Specifically, we use the data from HiTS that were observed between 2015 February 17 and22 (program 2015A-0608, PI: Förster).This region of the HiTS survey consists of 50 Galactic halo fields (∼150 deg 2 ) and includes 14 fields that were observed in previous HiTS campaigns (see Figure 1).The HiTS 2015 fields were observed up to five times per night, and are located between 137 and 160 deg in right ascension, and −7 and 2.6 deg in declination (Förster et al. 2016).The data were taken mainly in the −band, with 87 s exposures and a cadence of 1.6 hr.Observations in  the −band were performed as well, with individual exposures ranging from 81 to 102 s, which allowed the inclusion of − colours for our analysis.This configuration summed up a total of 20 to 29 epochs in , and from one to ten in , per field.For a more detailed description of HiTS' design, its observing strategy, and a comprehensive review of its characteristics we refer the reader to the work by Förster et al. (2016).It is worth mentioning that these data were not included in the work of Medina et al. (2018), where a previous HiTS campaign (from year 2014) was analyzed.
The second and third observing campaigns took place on halfnights between 2017 August 26 and 31 (programs 2017B-0907, PI: Muñoz, and 2017B-0253, PI: Carlin) and full nights between 2018 April 20 and 23 (2018A-0215, PI: Carlin, and 2018A-0907, PI: Muñoz), in the context of the HOWVAST survey (Medina et al. 2021).For HOWVAST we selected DECam fields to cover a con-siderable range of Galactic latitudes of the MW halo (from −43.0 to 28.7 deg).The footprint of HOWVAST was chosen to avoid a significant overlap with deep large-sky surveys, such as the DES (DES Collaboration 2005, 2016).In the first HOWVAST observing run, we observed 16 DECam fields during four consecutive halfnights (in first quarter moon phase), with 180 s exposures in the −band and a cadence of approximately one hour.The second HOW-VAST campaign consists of 24 fields, separated into two groups of 12 fields at different Galactic latitudes, observed in the −band during four consecutive nights (in first quarter moon phase).As in the first run, the integration times for this campaign were of 180 s, but with a cadence of ∼40 min.This results in a combined area of ∼120 deg 2 mapped in the halo surveyed by HOWVAST, and time series containing from 15 to 30 epochs per star.In addition, we obtained from two to four 240 s-exposure epochs per field in the −band each year in order to facilitate, thanks to the colour information, the identification of RRLs in our analysis.The coordinates of the HOWVAST fields in the equatorial system (J2000), including the number of observations per field, are provided in Table 1, making the distinction between the high-and low-Galactic latitude fields observed during our second campaign.
Combining the footprints of HiTS and HOWVAST, in this work we analysed an area of ∼270 deg 2 to search for RRLs.To increase the sample size of halo RRLs for our posterior analysis, we complement our detections with the RRL catalogue of Medina et al. (2018) (which covers ∼80 deg 2 of additional area).Therefore, in this work we analysed a total of ∼350 deg 2 .The sky coverage of the entire HiTS campaign and HOWVAST is shown in Figure 1 in equatorial and Galactic coordinates.

Pre-processing
The data reduction for this work was performed using the DECam community pipeline (Valdes et al. 2014).
The data from the HiTS 2015 campaign were pre-processed as part of the work by Medina (2017).A catalogue with the sources in these fields was created following the methodology from Medina et al. (2018).We first defined an x, y pixel coordinate system based on the output information generated by the SExtractor photometry software (Bertin & Arnouts 1996).To do this, we selected a reference frame for both the and the −band, for which the observing conditions were closer to optimal (i.e., from a photometric night and with the best seeing).Subsequently, we used the scaling constants found by the HiTS pipeline (Förster et al. 2016) to perform the alignment of the individual observations with respect to the reference.Then, we cross-matched the catalogues aligned in the common x, y coordinate system and rejected sources with fewer than five detections in the −band for the rest of the analysis.In order to keep sources with preliminary indications of variability for further processing, we disregarded sources for which the uncertainties in the mean flux exceeded by more than two times their flux standard deviations.Finally, we applied the x, y pixel to equatorial coordinate transformation using the astrometric solutions derived by Förster et al. (2016), which is computed from the crossmatch with known stars in the USNO catalog (Monet et al. 2003).
To pre-process the HOWVAST observations we adopted an alternative approach, based on the data processing pipeline in development for the Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration et al. 2009;Bosch et al. 2019).This pipeline was used to detect sources in the images, measure aperture fluxes, and perform a source point spread function (PSF) fitting.Because the sources of interest of this work are stars, we use PSF fluxes and magnitudes throughout the HOWVAST data treatment.For the subsequent variable star analysis, we only examined stellar sources with more than 15 data points in their time series, and with flux stellar deviations at least 2.5 times larger than their mean flux uncertainties.

Photometric calibration
In order to account for atmospheric effects affecting the epochs in our time series, we determined a photometric zero-point relative to the reference frame chosen for HiTS and HOWVAST separately.For this, we first computed instrumental magnitudes following where mag inst represents the instrumental magnitude either in the or the filter, Flux is the source flux, exp corresponds to the exposure time, and are the filter-dependent DECam photometric zero-point and first-order extinction coefficient per CCD1 (respectively), and is the airmass at the time of the observations.To calibrate the photometry of the HiTS 2015 fields, we first anchored our instrumental magnitudes to the reference frames in and .For the photometric calibration of the HOWVAST 2017 and 2018 data, we selected reference frames in which the PSF of the sources were minimum, similar to what was done for the HiTS 2015 data preprocessing.We compared our instrumental magnitudes with those in the reference frame, so that mag ref = mag inst + Δ rel , where Δ rel is the zero-point relative to the reference epoch, and mag ref is the object magnitude calibrated with respect to said epoch.
We then calibrated the photometry of the references using the archival data stored in the National Optical-Infrared Astronomy Research Laboratory (NOIRLab) Source Catalog (NSC; Nidever et al. 2021), as all the surveys considered in our work overlap the catalogs published in the second data release of such database.The photometric calibration of the NSC is based on the PS-1 survey, on the Skymapper and the Asteroid Terrestrial-impact Last Alert System (ATLAS) all-sky stellar reference catalogues (Wolf et al. 2018 andTonry et al. 2018, respectively), and on model magnitudes from linear combinations of photometry from catalogs such as the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) and the American Association of Variable Star Observers (AAVSO) Photometric All-Sky Survey (APASS; Henden et al. 2015).For the calibration, we limited the NSC data to star-like sources, with starClass flags larger than 0.85 (a starClass value of 0 is assigned for extended sources in the NSC, and a value of 1 is used for point-like sources).The sample was selected to include the best-matching NSC star within two arcseconds from each of the sources in our catalogue.We used only NSC stars with and magnitudes between 16.5 and 20.5, and magnitude errors smaller than 0.05, in addition to a two sigma clipping process performed over the median magnitude difference to remove outliers.From this comparison, we obtained an additional zero-point and colour term on a chip-by-chip basis, for each DECam field.Therefore, the calibrated magnitudes are given by mag where mag represents the calibrated magnitudes, NSC and NSC are the zero-point and the colour coefficient resulting from the magnitude comparison with the NSC data, respectively, and − is the colour of a given star.The average root mean square of these calibrations (across all fields and CCDs) is approximately 0.02 mag for the HiTS and the 2018 fields (in both bands), and of 0.03 mag for the 2017 fields.We note that these values are sufficiently small for the purpose of our work, but that calibrating our magnitudes onto Gaia synthetic photometry (see Gaia Collaboration et al. 2023) could contribute to an overall decrease if needed.Finally, the mean magnitudes were corrected for reddening using the dust maps of Schlegel, Finkbeiner, & Davis (1998), adopting the extinction coefficients (i.e., = 3.237 ( − ) and = 2.176 ( − )) from Schlafly & Finkbeiner (2011).Magnitude errors were computed by propagation of uncertainties.

Selection of the RR Lyrae
Since the data used in this work were obtained from two different surveys, we note that two slightly different methodologies were adopted to process the data for the search and characterization of RRLs.
The first sample corresponds to the RRLs identified in the HiTS 2015 survey by Medina et al. (2017).When looking for RRL candidates in this survey, objects with a magnitude variation smaller than 0.2 magnitudes were filtered out from the original source catalogues.Although applying this cut undoubtedly hinders the detection of low-amplitude RRLs, we consider this necessary in order to avoid contamination from non-variables in the faint end of the survey, where photometric uncertainties can easily reach 0.1 mag.Additionally, only sources redder than −0.2 in − , and bluer than 0.6 were considered for further processing.These values were chosen to remove potential spurious RRLs falling outside of the instability strip and to reduce the number of candidates to be inspected (as done by, e.g., Vivas et al. 2019).
The period of the sources in the HiTS 2015 catalogues were determined by Medina et al. (2017) based on the generalized version of the Lomb-Scargle period detection routine (GLS; Zechmeister & Kürster 2009), which incorporates a constant in the typical Lomb-Scargle sinusoid fitting procedure.By doing this, the results are overall less susceptible to aliasing and provide a more accurate frequency selection in the power spectrum.To compute the statistics and period selection, the GLS tool within the astroML Python module (VanderPlas et al. 2012) was used.Only sources with periods longer than 4.8 hours (0.2 d) and shorter than 21.6 hours (0.9 d), as well as those for which the GLS statistical level detections were smaller than 0.08, were considered in order to reduce the number of RRL candidates.We note in passing that these cuts should encompass most of the periods observed in RRLs (Catelan & Smith 2015) while avoiding aliasing around the period of one day and re-moving contamination from short period BL Her stars.For comparison, none of the RRL in the Gaia DR3 catalogue (Clementini et al. 2023) have periods shorter than 0.2 d, and only 0.13 per cent of them display periods longer than 0.9 d (0.05 per cent longer than 0.95 d).Finally, the two most significant periods were chosen and further inspected when more than one period were detected and met the aforementioned requirements.The last step for the selection of RRLs in the HiTS 2015 fields was to visually inspect the light curves resulting from the previous cuts and to look for objects with light-curve shapes, periods, and amplitudes typical of RRLs.The search resulted in a total of 95 RRLs in the fields that do not overlap with those inspected by Medina et al. (2018).
The selection criteria for RRLs in HOWVAST data are similar to the one followed for HiTS.We filtered out sources with minimum to maximum magnitude variability smaller than 0.2 magnitudes, and those with − colours clearly differing from the expectations for RRLs.In this case, given that a subset of the HOWVAST fields lie in regions with non-negligible interstellar extinction (at lower Galactic latitudes), we adopted a more generous cut and considered stars with − between −0.25 and 1.0 for further analysis.
To determine periods for the remaining sources, we used the Python package P4J2 , which was specifically designed for period detection on irregularly sampled and heteroscedastic time series, using the Cauchy-Schwarz Quadratic Mutual Information (QMI) as the criterion to be maximized by this routine (Huĳse et al. 2018).We first inspected the two periods with highest likelihoods, as long as they were longer than 0.2 d and shorter than 0.90 d and considered to have high statistical significance (at a 0.01 level).It should be noted that adopting these filters makes the most significant selection filter for HOWVAST's RRLs, namely the period detection criterion, comparable with the one used for HiTS' candidates.In fact, applying the QMI method to the RRL sample of HiTS yields periods that are indistinguishable from the GLS-based ones, with a median difference of the order of 10 −5 d.This similarity is not an unexpected result, as Huĳse et al. (2018) argue that the difference in the ability to recover the period of an RRL between both methods is minimized for time series containing over 20 epochs.We refer the reader to Huĳse et al. (2018) for a detailed description of the similarities and differences of both methods.
Finally, we visually inspected the phased light curves and selected only RRL-like objects, based on their light curve shapes, periods, and amplitudes.For stars exhibiting more than two high probability signals in the power spectrum, we examined the four most likely periods before choosing the star's main period.
The final list of RRLs from the HOWVAST data consists of 397 stars.Thus, by considering both RRLs from HiTS and HOWVAST, in this work we report the detection of a total of 492 RRLs.Their main properties are provided in Table A1.
In the remainder of this work, and for the sake of improving the number statistics of our analysis, we complement our sample with the RRLs from Medina et al. (2018).Hereafter, we refer to this updated catalogue as the full or combined sample.With this, we increase our catalogue size from 492 to a total of 663 RRLs.

Distance determination
We determined the absolute magnitude of our RRLs in the and bands ( and , respectively) using the period-luminositymetallicity relations for PS1 filters from Sesar et al. (2017) and as- where stands for the periods of the RRLs, and the second term takes into account the uncertainty of the zero point of the relations and their metallicity dependence.Given that these relations are only valid for fundamental-mode periods, for RRc stars we "fundamentalize" their periods prior to using Equation 3 by following: where F is the fundamentalized period (Catelan 2009).Heliocentric distances H are then computed through the distance modulus, and their uncertainties determined from error propagation.We note that computing metallicities from the Fourier decomposition of our light curves yields a median [Fe/H] of −1.4 dex (with individual uncertainties of ∼0.3 dex) when using the period-31 -[Fe/H] relation derived by Mullen et al. (2021).This approach, however, is not directly applicable to all of our RRLs due to the sparsely sampled light curves of a significant fraction of the sample.The effects of our metallicity assumption on the resulting distances in Equation 3 are expected to be small.In fact, an offset of 0.5 and 1.0 dex in [Fe/H] would lead to differences in H smaller than 4 and 8 kpc for remote RRLs (> 100 kpc), respectively.Figure 2 displays the heliocentric distance distribution of our RRLs and of those from the combined sample (i.e., considering the RRLs from Medina et al. 2018).Our sample consists of RRLs with H spanning from 7 to ∼270 kpc.Most of these stars lie within 50 kpc (429 RRLs; 87.2 per cent), whereas 54 of them have H between 50 and 100 kpc (10.9 per cent), and 11 (1.8 per cent) lie beyond 100 kpc.We further describe the most distant subsample in Section 4.An overdensity of 16 RRLs near 85 kpc (< >∼ 20.5 or (< >∼ 20.2) is clear from the figure, and is associated with RRLs in the Sextans dwarf spheroidal galaxy (dSph) that were not detected by Medina et al. (2018).We note that, similar to Medina et al. (2018), we find a discrepancy of ∼ 10 kpc between the mean distance to Sextans from RRab-and RRc-only samples (∼ 84 kpc and ∼ 75 kpc, respectively).This offset (of 11 ± 6 per cent) may be caused by the fact that the period-luminosity-metallicity relations from Sesar et al. (2017) are derived for fundamental-mode RRLs only or by differences in the passbands used in our works.Thus, we corrected the distances to the entire RRc stars in our sample (including those shown in Figure 2) by this factor.Vivas et al. (2019) searched for periodic variables in the Sextans dSph and all of our stars are found in their catalog when crossmatching within a radius of 1.5 arcsec, with the exception of the stars HiTS100752-020827 ( H ∼ 80.5 kpc), HiTS101128-  013921 ( H ∼ 82.8 kpc), HiTS101338-015258 ( H ∼ 80.1 kpc), and HiTS101734+001322 ( H ∼ 86.9 kpc), which are new detections.We highlight that our RRL classification matches with that of Vivas et al. (2019) for all the stars in common, and that the mean (absolute) period difference is of 0.008 d.Finally, we determined the period of three RRLs for which no period was reported by Vivas et al. (2019), namely HiTS101456-022025 (RR178 in their work, with a period of 0.4890 d), HiTS101342-021246 (RR142, 0.3999 d), and HiTS101551-015619 (RR192, 0.3232 d).

Completeness and recovery rate
Understanding how effective our methodology is in detecting RRLs and its limitations is crucial for interpreting the results of our search and for our RRL spatial distribution analysis.To quantify the effects of photometric completeness on our RRL detection efficiency, we take advantage of the fact that our data have been included The next step was to generate 5000 artificial RRL light curves mimicking the cadence and the expected photometric uncertainties (as a function of magnitude) of our survey.We then applied our photometric completeness (per epoch) estimates to the time series of a given light curve, based on the likelihood that an RRL would be observed at a given epoch (which is a function of magnitude as well).This typically results in a smaller number of epochs per light curve than in the ideal case (i.e., if the photometric completeness was 100 per cent), especially for fainter magnitudes.We run our selection pipeline and label as recovered the stars that pass our filters (see Section 3.1) and have a computed period within 10 per cent of the real (i.e., simulated) one.We find that, in this ideal scenario, we expect up to 90-95 per cent of the real RRLs to be recovered within 100 kpc.These numbers decrease significantly for sources fainter than > 21.5 (∼ 150 kpc), where the recovery rate drops to < 80 per cent.However, since there are other effects in play for real observations of individual sources (e.g., blending of sources, proximity to bright sources, and their positions in the CCDs), the number of recovered RRLs should be slightly smaller overall than predicted by our idealized simulation.

Comparison with previous surveys
To assess the fraction of RRLs from other surveys that might be missing in our catalog we crossmatch our sample with the RRLs in the Catalina Real-time Transient Survey (CRTS; Drake et al. 2013aDrake et al. , 2014Drake et al. , 2017;;Torrealba et al. 2015), the PS-1 (Sesar et al. 2017) survey, and with the Gaia catalogue generated with the SOS pipeline.For the latter, we use the recently published catalogue based on Gaia DR3 (Clementini et al. 2023), which almost doubles the size of its predecessor, Gaia DR2 (Clementini et al. 2019).We perform our crossmatches using a search radii of 3 arcsec to account for the different pixel scales of these surveys.For these comparisons, we only considered the RRLs from the literature that fall within the footprint of our survey, excluding the DECam's CCD gap regions.We emphasize that, by design, HOWVAST does not significantly overlap the area covered by large surveys such as the DES (Stringer et al. 2021) and the ZTF (Chen et al. 2020;Huang & Koposov 2022) survey.The results of our comparisons, including the few overlapping RRLs from the DES and the ZTF are summarized in Table 2.
We restrict our comparison to a meaningful magnitude (distance) range to avoid using stars in the CRTS, PS-1, or Gaia that could have saturated epochs in our survey ( 20 kpc).When limiting the comparison to RRLs between 20 and 40 kpc, observed at least 70 pixels from the edge of the CCDs (and with detections in the and bands, for Gaia), the number of recovered RRLs is close 70-75 per cent in each case with the exception of the ZTF, where we recover ∼55 per cent of the RRLs.If we expand the distance range to the region spanning from 40 to 80 kpc, we are able to recover between 41 and 67 per cent of the RRLs (in the case of Gaia and the CRTS, respectively).A tentative explanation for this difference, albeit the low number statistics, is the possible contamination in the Gaia SOS catalogue at these distances, as it has been shown that artefacts and spurious detections might be present in crowded areas (e.g., close to the Galactic plane; Holl et al. 2018;Clementini et al. 2019;Rimoldini et al. 2019), in combination with our survey's photometric completeness limitations and RRL selection strategy.Beyond 80 kpc, no RRLs from the CRTS that fulfill our selection cuts are found, while we recover ∼50 per cent of the RRLs from PS-1 and Gaia.
We find that 99 (98) per cent of the stars classified as RRab (RRc) in our catalogue have the same classification in the Gaia catalogue, with a median absolute period difference of 0.0026 (0.0009) d.Similarly, 99 (100) per cent of the ab-type (c-type) RRLs identified in our work have the same classification in the CRTS, in which case the median absolute period difference is of 0.0020 (0.0008) d.Thus, we consider our classification methodology to be robust.
We detect 103, 89, and 87 RRLs that are not listed in the concatenation of the aforementioned catalogues, when crossmatching using a search radius of 1, 3, and 5 arcsec, respectively.These stars are located from ∼7 to 265 kpc in heliocentric distance, and the majority of them (58 per cent) lie in the low Galactic latitude fields of the second HOWVAST campaign.Interestingly, 80 per cent of the 69 RRLs with H ≤ 80 kpc are classified as c-type, which might be a consequence of contamination by blended sources (for RRLs near the Galactic plane) and/or the misclassification of variable objects (e.g., eclipsing binaries).This does not occur for the stars further than 80 kpc, where 60 per cent of the new RRLs are classified as ab-type.We highlight that the majority of the most distant stars in our catalog are new discoveries.We describe this subsample in more detail in Section 4.

Classification and Bailey diagram
In order to classify the RRLs in our catalogue, we adjusted the light curve templates from the SDSS Stripe 82 (Sesar et al. 2010) to our phased light curves.This was performed using the templates in the and bands available in the Python package (VanderPlas & Ivezić 2015).The final classification of RRLs into aband c-subtypes was based on the inspection of the best-fitting templates, their amplitudes, and periods.This resulted in 331 RRab stars (67 per cent), 157 RRc stars (32 per cent), and 4 RRLs (1 per cent) that do not fall in either category, with indications of pulsations in the fundamental-mode and first-overtone.We classified the latter as RRd stars.The relative fractions are similar if we consider the full sample, where 69 per cent and 30 per cent of the stars are classified as RRab and RRc, respectively.The distribution of these stars in the period-amplitude space (Bailey diagram) is shown in Figure 4 colour-coded by type.
The position of RRLs in the Bailey diagram can be used to confirm their classification and as a tool to assess the Oosterhoff type (Oosterhoff 1939) of stellar systems.The Oosterhoff groups are seen as a dichotomy in the mean period and amplitude, and the ratio of RRab and RRc stars in globular clusters.Cluster RRLs can be split into Oosterhoff-I type (of shorter RRab star periods ∼0.55 d and lower RRc star fractions; OoI), Oosterhoff-II (RRab star periods ∼0.65 d and higher RRc star fractions; OoII), and those that lie in between the two regimes (Oosterhoff-intermediate; Oo-int).This dichotomy, however, is not present in most MW dwarf galaxies and field stars (the latter being predominantly OoI).
Figure 4 shows the OoI, OoII, and Oo-int fiducial lines in the Bailey diagram as defined by Fabrizio et al. (2019) for −band RRL amplitudes.To account for the differences between RRL amplitudes in the and −bands from our work and those in the −band, we scale our amplitudes by a factor of 0.90 and 1.21, respectively (Sesar 2012).From the orthogonal proximity of the stars to either of the Oosterhoff groups' fiducial lines, we conclude that most of the RRab stars in our sample could be considered OoI or Oo-int, which is the expected trend for field RRLs.We discuss in more detail the distribution of the most distant RRLs in the Bailey diagram in Section 4.

RR LYRAE STARS BEYOND 100 KPC
We identify 9 RRLs with mean and magnitudes fainter than 20.7 and 20.3, respectively, which corresponds to heliocentric distances larger than 100 kpc.When considering the full sample, the number of RRLs beyond this limit increases to 23.The spatial distribution of these stars is depicted with larger symbols in Figure 3. None of the stars identified in this work are listed in the catalogues used for comparison in Section 3.4.Among these RRLs, three stars are located beyond 200 kpc.The number of epochs in the light curves of these RRLs spans from 22 to 27 in the −band, and from 15 to 28 in the −band.The folded light curves of the newly detected RRLs are shown in Figure A1, and their main properties are summarized in Table 3.
Most stars in this subsample are classified as ab-type (83 per cent), which is not surprising given their relative abundance and the fact that they are easier to identify than RRc or RRd.The location of these distant RRab stars in the Bailey diagram (Figure 4), does not show a strong association with the locus of the fiducial line of the OoII group (nor the OoI group), albeit their tendency for periods of pulsation longer than 0.60 d.In fact, the average period of these RRab stars is 0.66 d, similar to the mean period of distant RRLs found by Medina et al. (2018) and the collection of RRLs in MW ultrafaint dwarf galaxies studied by Vivas et al. (2016) (0.67 d; see also Martínez-Vázquez et al. 2021).Figure 4 shows that the distribution of distant RRLs does not follow the overall trend of RRab stars within 100 kpc.In particular, they are not preferably located near the locus of the OoI group.
High Amplitude Short Period (HASP) variables, that is, those with periods shorter than 0.48 d, and −band amplitudes larger than 0.75, have been interpreted as coming from progenitors or regions in the Galaxy with populations more metal-rich than [Fe/H] = −1.5 (Fiorentino et al. 2015).Therefore, RRLs lying in this region of the Bailey diagram can provide insights in the building of the halo and its progenitors.In fact, most MW dSphs lack HASP variables, while these stars are not rare in the halo and among globular clusters more metal rich than −1.5 dex and massive dwarf galaxies (Fiorentino et al. 2017).In the combined sample, we find 26 RRab stars populating the HASP region, which corresponds to only 6 per cent of our full RRab star sample.The heliocentric distance of these stars (see Section 3.2) ranges between 9 to 45 kpc, thus none of the RRLs in the most distant subgroup lie in the HASP region.The relatively low fraction of HASP RRLs in our sample might be an indication of the contribution of dwarf galaxies to the dual origin of the outer halo (and its dependence on Galactocentric distance), as further discussed in Section 5. We note, however, that additional evidence is required to support this assertion.
To examine the overall consistency of the number of distant RRLs found in our work with the results of previous studies of similar photometric depth, we can perform a direct comparison assuming high completeness out to similar distances.Stringer et al. (2021) detected 6,971 ab-type RRL candidates in the ∼5,000 deg 2 of the DES' footprint, among which 4,569 do not belong to the known substructures and galaxies considered by the authors.Of this subsample, 18 per cent are located beyond 100 kpc from the Sun, which implies a rough density of six distant halo RRLs every 40 deg 2 (or 0.16 per deg 2 ), without accounting for their estimated completeness (expected to be > 70 per cent at ∼ 150 kpc).This number is a factor of three larger than the number of distant RRLs detected in our study (two RRab stars beyond 100 kpc every 40 deg 2 , or 0.056 per deg 2 ).Nonetheless, our density is more consistent with the findings of Stringer et al. (2021) if we only consider their candidates with more than 25 observations in total (considering , , , , and ) and with an RRab score > 0.90 as assigned by their classifier.In this case, the DES RRab star density decreases to 0.083 distant RRLs per square degree (or three RRLs every 40 deg 2 ), in rough agreement with our results.

Potential groups of distant RR Lyrae stars
In Figure 3, we highlight the spatial distribution of the RRL candidates detected beyond 100 kpc.From the figure we identify two groups of stars with similar on-sky positions and heliocentric distances.Associations and groups at large distances (especially at H > 100 kpc) are unlikely to happen by chance from halo stars.Using mock stellar haloes and focusing on RRLs beyond 100 kpc, Sanderson et al. (2017) showed that the median of the minimum angular distance to the nearest star for bound (unbound) RRLs beyond 100 kpc is ∼0.01 deg (3.0 deg), with the majority of minimum separations between 0-0.03 deg (∼0.3-10.0deg).Moreover, Sanderson et al. (2017) found that the closest pairs of RRLs tend to originate from the same building block (regardless of their bound/unbound status).
We analyzed these groups looking for indications of their potential association with known substructures.The main properties of the RRLs in these groups are presented in Table 4.The first group consists of the stars HV153403-321831 and HV152905-315335, two RRab stars located at 135 ± 5 and 144 ± 5 kpc, respectively (with right ascensions of ∼233 deg).These stars are separated by 1.1 deg, which corresponds to ∼3 kpc at H ∼ 140 kpc.The second group, at right ascensions of ∼237.5 deg, comprises the stars HV154834-320810 and HV155407-361645, with distances of 110 ± 4 and 104 ± 4, respectively.Both of these stars are classified as RRab stars.This group shows an angular extension of ∼4.3 deg (or 8 kpc at H ∼ 107 kpc).We note that the stars in both groups are too separated to be considered part of an intact (or not heavily disrupted) bound satellite, but they lie well within the predicted range of separations for unbound (but associated) debris from Sanderson et al. 2017.
As these RRLs might be associated with the ongoing tidal disruption of MW satellites, we inspected the Python library galstreams (Mateu 2023), which collects celestial, distance, proper motion, and radial velocity information for ∼97 distinct stellar streams.Nevertheless, we find no streams within the galstreams database with distances similar to those of our groups (most of the streams close to the position of our groups have H < 40 kpc).Comparing the positions and distances of these groups to the model of the Sagittarius stream by Dierickx & Loeb (2017) shows that only the group with higher right ascension is in proximity to the stream, but at larger distances (most Sagittarius stream stars at this right ascension are located at H ∼ 50 kpc ).Moreover, all of the stars in these groups have latitude-like coordinates in the Sagittarius stream system ( Sgr ; Majewski et al. 2003) larger than 16.6 deg, making their association with the stream unlikely.Therefore, we find no clear indications of associations between our groups of clumped RRLs and known satellites or streams.Nevertheless, we suggest that the association of the stars in these groups is likely.

Radial density model
In this section we address in more detail the radial density distribution of our full sample.
Many observational studies and simulations have suggested that the properties of the radial distribution of stars in the halo is connected with their origin (e.g., Zinn 1993;Vivas & Zinn 2006;Pillepich et al. 2014).From the slope(s) of the radial distribution, for instance, one can infer the existence of an inner halo thought to contain both accreted and formed in-situ stars, and an outer halo, expected to have been formed largely from the accretion of satellites (e.g., Watkins et al. 2009;Bullock & Johnston 2005;Zolotov et al. 2009;Naidu et al. 2020) In order to characterize the spatial distribution of our RRLs in the halo, we follow the methodology of Medina et al. (2018).We adopt two models to represent the data -one assuming a spherical halo, and one for an ellipsoidal halo.To account for the oblateness of the latter, we assume = / = 0.7 (Sesar, Jurić, & Ivezić 2011), where and are the axes in the disc plane and along the vertical direction, respectively.Thus, prior to binning our sample in distance, we transform our computed distances from heliocentric to spheroidal and ellipsoidal Galactocentric distances ( sph and el ) using: where ⊙ stands for the distance from the Sun to the Galactic centre (assumed to be 8 kpc; Gravity Collaboration et al. 2021), and and are the Galactic latitude and longitude, respectively.The adoption of two fixed flattening parameters is motivated by our limitations in constraining the shape of the halo from our survey's small (and non-contiguous) footprint, which would lead to over-interpretations of our results.
We adopt a power-law model to describe the radial density ( GC ) of our halo RRLs, where GC represents either sph or el .Thus, ( GC ) = ⊙ ( GC / ⊙ ) , where ⊙ is the local RRL number density and is the slope of the profile for a simple power law.Given that vast observational evidence suggest the existence of a break in the halo radial density profile between 20 and 35 kpc (e.g., Saha 1985;Watkins et al. 2009;Deason, Belokurov, & Evans 2011;Sesar, Jurić, & Ivezić 2011;Belokurov et al. 2018a;Stringer et al. 2021), we describe the explored regions with simple and broken power-law models.For the latter, a break at break is used, so that, in logarithmic form: where = log ( ⊙ ), and the subindices denote each side of the density profile (i.e., 1 and 2 correspond to the inner and outer density, respectively).
To explore the parameter space and their distribution, we employ emcee (Foreman-Mackey et al. 2013), a Python implementation of the invariant Markov chain Monte Carlo method (MCMC).For this, we leave 2 , 1 , 2 , and break as free parameters, and adopt the priors used by Medina et al. (2018).We find that running emcee with 200 walkers and a chain of 500 steps is sufficient to reach convergence.The selected values are provided in Table 5 and correspond to the median of the marginalized posterior parameter distributions, and their errors represent their 16th and 84th percentiles.The posterior probability distribution for the obtained parameters of the brokenpower-law model of the oblate halo are depicted in Figure A2.We examine our findings and their implications in the following section.

Number density profiles
We report the results of the aforementioned methodology for our catalogue excluding only the subsample of stars belonging to the Sextans dSph.Based on the estimations of our detection efficiency and completeness, we constrain our analysis to RRLs with distances smaller than 145 kpc, the limit at which the idealized RRL recovery rate is ∼80 per cent and beyond which this rate decreases drastically.Additionally, we correct our densities by the RRL detection completeness values obtained in Section 3.3.While our RRL recovery rates from comparisons to other surveys are somewhat lower than the completeness expected based on our idealized estimates, we use the completeness values based on our dataset.We choose this because our data are deeper than any of the literature surveys, and thus expected to be more complete.Furthermore, we wish to avoid the complicating effects of unknown rates of spurious RRL detections and incompleteness in existing surveys.In order to inspect the (in-)homogeneity of the RRLs radial distribution, we follow two approaches: one computing the density profile of our entire sample of RRLs, and one measuring the profile in the four distinct regions covered by our survey.By doing this we could, in principle, directly inspect anisotropies in the halo distribution of stars at relatively small scales, albeit with the additional challenge of low-number statistics.Lastly, for each region, we also included an analysis of the distribution of ab-type RRLs only, as they are typically less affected by contamination and are more uniquely identifiable than c-type and d-type RRLs.
Our results for a spherical halo model and no RRL class distinction are shown in Figure 5, and a summary of the overall results is provided in Table 5.A first look into Figure 5 and Table 5 reveals that, regardless of the region considered, a break is visible between 15 and 25 kpc.This feature is, however, less clear in the fields closer to the Galactic plane (2017 and 2018B, see Figure 1), where more RRLs are detected overall at distances < 20 kpc.This is also the case if we restrict our sample to ab-type RRLs only.Because ab-type RRLs are less prone to contamination and the vast majority of the stars in our sample have Galactic coordinates | | > 2 kpc, the observed difference is likely due to the overall higher density of stars towards the Galactic plane (as a reference, Mateu & Vivas 2018 estimated the thick disc scale height to be ℎ ∼ 0.65 kpc using RRLs).Based on the resulting reduced 2 ( 2 ) values of the models, we conclude that broken-power-law models yield better fits than simple-power-laws.The break in the profile is also clearly observed in the plot of Figure 5, depicting the radial distribution of RRLs when using all the regions combined.In the latter, we see that the underdensities and overdensities at given distance bins visible in the individual region profiles get averaged out, as a consequence of better number statistics.In that case, the best-fit case corresponds to the spherical halo and we find break at 18.1 +2.1 −1.1 kpc, where the profile displays an inner slope of 1 = −1.88+0.23  −0.16 and a steeper outer halo slope of 2 = −4.47+0.11 −0.18 .We find that the break radii from the different regions are con- sistent within their uncertainties (Table 5), whereas 2 displays a larger dispersion overall.For a spherical halo, these values vary between 17.4 +2.7 −1.4 and 21.5 +2.0 −2.1 kpc for the break , and −4.96 +0.51 −0.52 and −4.49+0.47  −0.47 for the outer halo slopes.The largest differences are seen when contrasting the inner slopes 1 , which tend to be steeper for the fields near the Galactic plane.
In Table 6, we display the slopes of the radial distribution of different tracers reported in the literature, and in Figure 6 we con-trast our best-fit parameters with broken-power-law fits from these works.These works generally model the radial distribution adopting (or deriving different halo flattenings, and cover large areas and different regions of the halo.In this regard, we note that our survey maps a smaller area than previous studies, but we are able to better trace the halo beyond 40 kpc.We find that our measured break radii are consistently smaller than those found by Stringer et al. (2021), Das, Williams, & Binney (2016), and Chen et al. ( 2023) (who find break closer to ∼30 kpc), regardless of the adopted halo flatten-Table 5. Simple and broken-power-law parameters from the sampled posterior probability distributions described in Section 5. We report the results from using the RRLs in all of our regions combined, and from the individual areas of the survey, as well as the adopted shape of the halo (spheroidal and oblate, denoted as sph and el , respectively).For most of the studied regions the best-fitting model corresponds to the broken-power-law profiles.stars).Regarding the inner slope values, we observe that our results are higher than those from the literature overall.We attribute this to our inability to map the inner halo and contamination in our sample within 20 kpc (in particular, the misclassification of RRc near the Galactic plane).The measured outer slopes from each of our regions are broadly consistent with previous studies based on RRLs, and our value from using all the fields lies within 2 of the 2 of most of the works shown in Figure 6.Our values of 2 are also consistent with the SPL density profiles of recent works that focus on tracers at large distances but are unsuccessful at accurately mapping the inner halo (and thus, are unable to find a break).This is the case of the studies of RRLs and BHB stars by Medina et al. (2018) and Yu et al. (2024), respectively (which are not shown in Figure 6 for consistency), and most remarkably those of Lopez-Corredoira et al. (2024)'s with K and M giants.We note, however, that our density outer slopes are systematically steeper than those from previous works mapping the distant halo ( GC > 80 kpc) with BPLs using different halo tracers (e.g., Thomas et al. 2018;Fukushima et al. 2019), and the recent RRL-based SPL outer halo exploration by Feng et al. (2024).This might be caused by contamination in their tracer samples (e.g., by blue stragglers in BHB star catalogues), an incomplete census of RRLs with large distances in our catalogue, or intrinsically differ-ent radial distribution for different stellar samples.We highlight the similarity of our outer halo results regardless of whether we use RRab-only samples or all RRL types, which is likely a consequence of the high number fraction of RRab in our sample.

Density profiles in context
The radial density profile of RRLs that we analyse in this work extends to larger radii than any previous study of similar or larger surveyed area.Additionally, the extrapolation of our measured slope beyond the 145 kpc radius used for MCMC modeling provides a good match to the outermost data points in Figure 5, suggesting that the slope we measure is appropriate even in the outermost halo (beyond ∼200 kpc).While the ∼350 deg 2 area covered by our study is much smaller than that used in many other measurements of halo density profiles, our results are consistent with those of most previous studies (see Figure 6).We also note that our sample of RRLs with relatively well-sampled lightcurves is unlikely to suffer significant contamination, and thus constitutes a robust and secure sample of outer halo stars.Figure 7 illustrates the total number of stars beyond a given radius for various power-law outer halo density profiles.For our best-fit oblate halo model, 2 = −4.47, a total of ∼3000 RRLs are expected over the entire sky beyond 100 kpc, or roughly 0.073 per square degree This is comparable to the model-based expectations from Sanderson et al. 2017, which predict ∼2000-6000 RRLs from unbound merger remnants beyond 100 kpc in the MW (see their Figure 5).For our survey of 350 deg 2 , we would thus expect to find 25 RRLs at GC > 100 kpc, while we find a total of 23.Interestingly, the right panel of Figure 7 suggests that the mere identification of two secure RRLs beyond 230 kpc in our 350 deg 2 study rules out outer halo slopes with −5 or −3.5.Over the entire sky, an outer

DISCUSSION AND SUMMARY
We have described the search for RRLs in different directions of the remote MW halo using DECam data from the HiTS and the ongoing HOWVAST surveys.We construct light curves from time series containing from 16 to 38, and from 15 to 32 epochs in the and −bands, respectively.Considering all of the studied fields in HiTS and HOWVAST, we detect a total of 492 RRLs (397 in the HOWVAST fields and 95 in the HiTS fields) in a combined area of ∼270 deg 2 , including at least 90 RRLs not listed in existing surveys.When combining our catalog with that of our previous study (Medina et al. 2018), the number of RRLs increases to 663 over a total area of ∼ 350 deg 2 .The heliocentric distances of our RRLs range between 7 and 270 kpc.We identified 9 new RRLs beyond 100 kpc, which we add to the still small list of well-characterized tracers of the old component of the MW at large distances.
The distribution of distant RRab stars in the Bailey diagram is not preferentially located in a unique Oo-group, and they are rather uniformly distributed in the period-amplitude space.In this subsample, we identified two potential RRL associations containing stars with similar distances (both at around 110 kpc) located within a few degrees from each other, and with a mean period of 0.71 d, consistent with them being linked to the OoII group.This might be an indication of them being accreted material from ultra-faint dwarf galaxies.Moreover, previous studies have shown that neighboring stars at these distances are unlikely to occur by chance (Sesar et al. 2014;Baker & Willman 2015;Sanderson et al. 2017), which makes them potential tracers of known or undiscovered substructures.We found that the position of these groups is inconsistent with those of previously known streams (e.g., the Sagittarius stream), and cannot directly associate these stars with the accretion of UFDs with our data only.We conclude that the stars in these groups are likely associated, and advocate for additional data and follow-up studies to confirm their association and to determine their parent populations.
We characterized the (radial) spatial distribution of our RRLs with power-law profiles, by adopting a spherical halo model ( = 1) and an ellipsoidal model (oblate, with = 0.7).Furthermore, we analyze the density profiles from the RRLs in our entire sample, and from different directions in the halo.For this, we followed an MCMC approach and consider RRLs located at distances <145 kpc from the Galactic centre for our modeling.We found that the profiles are better described by broken-power-laws, as it has been shown by other works in the literature.For the preferred model (spherical halo), our best-fitting results suggest a break in the RRLs distribution at 18.1 +2.1 −1.1 kpc, with an inner slope of −2.05 +0.13 −0.15 , and a steeper outer slope of −4.47 +0.11  −0.18 .Stellar haloes are important testbeds sensitive to various aspects of galaxy formation models, and comparing observations (e.g., the properties of their density profiles) with simulations is an important requirement to draw meaningful conclusions.In recent years, several authors have measured the stellar distribution of MW-like galaxies using sophisticated cosmological simulation suites (e.g., the Illus-trisTNG project; Pillepich et al. 2018).In these simulations, outer halo slopes are typically found to be in the range −5.5 < < −3.5, where recently formed haloes or those with a large fraction of their total stellar mass originating from mergers have shallower slopes, and steeper slopes correspond to quiescent recent accretion histories (Pillepich et al. 2014).Our measured density profiles are remarkably consistent with the results reported by Merritt et al. (2020), who predicted a median outer slope (beyond 20 kpc) of −4.5 for MW-like galaxies using the TNG100 simulation of the IllustrisTNG project (e.g., Pillepich et al. 2018;Nelson et al. 2019).
Better number statistics (with large scale halo surveys of high completeness) resulting from the advent of large and deep photometry campaigns will be crucial to reconstruct a more complete local and global picture of the outer halo's history, structure, inclination, and shape.These studies will permit us, e.g., to assess the amount of variation in stellar density at large radii (where asymmetries are expected to be more evident; e.g.Pandey 2022) and to better characterize the disequilibrium state caused by the dynamical response of the Galactic halo to the infall of massive satellites (observed over thousands of square degrees; see e.g.Conroy et al. 2021;Han et al. 2022;Rozier et al. 2022).In this regard, the upcoming ten-year Rubin Observatory Legacy Survey of Space and   ) with a blue solid line.These models were obtained from the SDSS Stripe 82 RRLs light curve templates (Sesar et al. 2010).The main properties of these RRLs are summarized in Table 3.

Figure A2.
Corner plot of the posterior probability distributions for the broken-power-law profiles described in Section 5.The adopted value of each parameter is the median of the corresponding marginalized distribution, and their uncertainties represent the 16th and 84th percentiles.These parameters are computed for the RRLs in all the fields studied in this work.

Figure 1 .
Figure 1.Spatial distribution of the surveys used for this work, shown in both equatorial and Galactic coordinates.The DECam fields corresponding to the HOWVAST survey are plotted in red and green, while the fields observed by the HiTS survey are shown in blue.

Figure 2 .
Figure 2. Distribution of the heliocentric distances H of the RRLs analysed in our work.The red-filled region represents the distance of the stars detected in this work (as described in Section 3.2), whereas black bars depict the distribution of the full RRL sample, including those from Medina et al. (2018).

Figure 3 .
Figure 3. Spatial distribution of the RRLs analysed in this work, colour-coded by their heliocentric distance H in kpc.The stars from the HiTS campaign analyzed by Medina et al. (2018) are plotted with small black symbols, and the RRLs from this work with distances larger than 100 kpc (described in Section 4) are plotted with large colour-coded circles.An approximation of the footprint of each DECam field observed is shown in gray in the background as a reference.Enlargements of the two regions containing RRLs beyond 100 kpc with similar distances (and potentially associated with each other) are provided as insets.A third enlarged region depicts the overdensity associated with the Sextans dSph.

Figure 4 .
Figure 4. Bailey diagram of the RRLs analysed in this work.RRab stars are displayed with red pentagons, whereas light-blue and brown symbols represent c-type and d-type RRLs, respectively.The amplitudes represent the minimum-to-maximum variation of the fitted light curves in the band, obtained by scaling the and magnitude amplitudes by a factor of 0.90 and 1.21, respectively (Sesar 2012).Stars with estimated heliocentric distances larger than 100 kpc are plotted with black symbols.The dashed regions depict the fiducial lines for OoI, Oo-int, and OoII defined by Fabrizio et al. (2019).Black solid lines delimit the region containing HASP RRab variables (as defined by Fiorentino et al. 2015), which is shown as a grey shaded area.
in an early version of Data Release 3 of the DECam Local Volume Exploration (DELVE) survey (for an overview of DELVE, seeDrlica-Wagner et al. 2021).DELVE includes new data observed specifically for the survey, but also aggregates those new data with all existing DECam data in the public archive.Thus, our HOWVAST RRL data are included in DELVE DR3 processing, wherein all available data are coadded to create deep, stacked images.We estimate our survey's photometric completeness by comparing our single exposure catalogues with DELVE's coadded catalogues.For the latter, we consider point sources only ("high confidence stars," following the definition of EXTENDED_CLASS from Section 4.7 ofAbbott et al. 2021), with DELVE photometric uncertainties < 0.5 mag and number of epochs greater than zero.Then, we crossmatch these coadded catalogs with ours and compute the fraction of recovered stars as a function of magnitude.Our survey recovers over 95 per cent of the DELVE sources brighter than ∼ 22 mag consistently across epochs and the recovery fraction drops to ∼ 50 per cent at ∼ 23 mag.The photometric completeness curves obtained with this strategy are very similar to those of the HiTS survey presented byMartínez-Palomera et al. (2018) (Figure2in their paper), which is not surprising given the similarity in observing conditions, strategies, and designs between both surveys.

Figure 5 .
Figure 5. Binned RRL number density profiles of the regions studied in this work, for a spherical halo model.The region corresponding to the HiTS fields is shown with blue curves, whereas that of HOWVAST 2017 is shown in orange.Green curves depict the high and low Galactic latitude areas surveyed by HOWVAST 2018 (left and right, respectively).The profile resulting from considering the entire studied region is shown in black.The best solution determined via MCMC is shown as a solid line in each panel, and the shaded regions depict the 3 confidence levels.The uncertainty shown for each density bin represents Poisson noise.

Figure 6 .
Figure 6.Comparison between our best-fitting radial density profile parameters for a BPL spherical halo model, and those from other BPL halo studies in the literature.Parameters obtained from the spatial distribution of RRLs only are shown with circles, whereas triangles represent those from studies based on different halo tracers (e.g., BHB stars or K-giants).The markers are colour-coded to illustrate the distance limits of each work, and the error bars depict the uncertainties in the parameters estimation, when available.We denote our values as TW ab for the parameters from our RRab-only MCMC model, and TW all for the model based on our full RRL sample.A horizontal dashed line is used to highlight the results of the latter.The references for the literature works are defined as follows: K08 (Keller et al. 2008), W09 (Watkins et al. 2009), S11 (Sesar, Jurić, & Ivezić 2011), D11 (Deason, Belokurov, & Evans 2011), A12 (Akhter et al. 2012), F14 (Faccioli et al. 2014), Z14 (Zinn et al. 2014), P15 (Pila-Díez et al. 2015), X15 (Xue et al. 2015), D16 (Das, Williams, & Binney 2016), T18 (Thomas et al. 2018), S21 (Stringer et al. 2021), and C23 (Chen et al. 2023).We note that the break values found by K08 (45 kpc), A12 (45 kpc), and T18 (41.4 kpc) are not shown in the upper panel, for illustrative purposes.

Figure 7 .
Figure 7. Expected number of RRLs for density profiles with different slopes (for ellipsoidal halo models).Left: Predictions over the entire sky, showing radial density slopes of = −3, −4, −5, −6, and our measured value of 2 = −4.47.A total of ∼3000 RRLs are predicted at a distance larger than 100 kpc based on our fit, and only ∼350 at distances > 230 kpc.Right: Predicted number of stars in a survey similar to ours, with area of 350 deg 2 .Open circles show the actual number of outer halo RRLs analyzed in our study.Finding two RRLs beyond 230 kpc rules out the most extreme slopes illustrated here, as a shallow slope (e.g., −3.5) would predict hundreds of stars at these distances and steep slopes (e.g., −5.5) would make it unlikely to find any RRLs at GC > 230 kpc in a 350 deg 2 survey.In both panels, grey lines are used to represent 1000 random drawings from the MCMC chains described in Section 5.1.

Figure A1 .
Figure A1.Folded light curves in the and band of our sample of RRLs with H > 100 kpc, and out to ∼270 kpc.For each star, we overplot the best-fitting model from the Python module (VanderPlas & Ivezić 2015) with a blue solid line.These models were obtained from the SDSS Stripe 82 RRLs light curve templates(Sesar et al. 2010).The main properties of these RRLs are summarized in Table3.

Table 1 .
Identification numbers and equatorial coordinates of the DECam fields observed by HOWVAST in 2017 and 2018.We quote the central J2000 coordinates of our fields and the number of observations per field in the and bands (N and N , respectively).Each field covers three square degrees, approximately.The fields of the second campaign are labeled according to their positions with respect to the Galactic plane (2018A and 2018B represent the high-and low-Galactic latitude fields).

Table 2 .
Total number ( ) and recovered RRLs ( ) from previous surveys as a function of heliocentric distance H .This table shows that, while our survey is partially limited to sources with mean magnitudes fainter than 17.5 in and , we are able to recover over 70 per cent of the RRLs beyond 20 kpc.

Table 3 .
Medina et al. (2018)he sample of distant RRLs discussed in Section 4, including equatorial coordinates, mean magnitude and number of epochs per filter, and heliocentric and Galactocentric distances.We flag the stars detected in this work with zeros and those fromMedina et al. (2018)with ones.
* The period and amplitude of pulsation are computed from the photometric band with more observations.

Table 4 .
Main properties of the potential groups of RRLs beyond 100 kpc discussed in Section 4.1.

Table A1 .
Same as Table3but for the entire sample of RRLs analysed in this work.The full table is provided as supplementary material and will be made available at the Strasbourg astronomical Data Center (CDS).The period and amplitude of pulsation are computed from the photometric band with more observations. *