Small-scale stellar haloes: detecting low surface brightness features in the outskirts of Milky Way dwarf satellites

Dwarf galaxies are valuable laboratories for dynamical studies related to dark matter and galaxy evolution, yet it is currently unknown just how physically extended their stellar components are. Satellites orbiting the Galaxy's potential may undergo tidal stripping by the host, or alternatively, may themselves have accreted smaller systems whose debris populates the dwarf's own stellar halo. Evidence of these past interactions, if present, is best searched for in the outskirts of the satellite. However, foreground contamination dominates the signal at these large radial distances, making observation of stars in these regions difficult. In this work, we introduce an updated algorithm for application to Gaia data that identifies candidate member stars of dwarf galaxies, based on spatial, color-magnitude and proper motion information, and which allows for an outer component to the stellar distribution. Our method shows excellent consistency with spectroscopically confirmed members from the literature despite having no requirement for radial velocity information. We apply the algorithm to all $\sim$60 Milky Way dwarf galaxy satellites, and we find 9 dwarfs (Bo\"otes 1, Bo\"otes 3, Draco 2, Grus 2, Segue 1, Sculptor, Tucana 2, Tucana 3, and Ursa Minor) that exhibit evidence for a secondary, low-density outer profile. We identify many member stars which are located beyond 5 half-light radii (and in some cases, beyond 10). We argue these distant stars are likely tracers of dwarf stellar haloes or tidal streams, though ongoing spectroscopic follow-up will be required to determine the origin of these extended stellar populations.


INTRODUCTION
The early stages of galaxy formation involve the accumulation of baryonic and non-baryonic matter at conjunctions of the cosmic web, where the first proto-galaxies are suggested to have formed (Mo et al. 2010).Baryons are accreted into dark matter haloes, and these halos grow hierarchically (White & Rees 1978;Frenk et al. 1988).This general framework applies to galaxies of all sizes, including typical galaxies like the Milky Way (MW;  * ≈ 5 × 10 10  ⊙ ) and Andromeda (M31), down to the least massive systems (ultra-faint dwarfs, UFDs;  * ≲ 10 2 − 5  ⊙ ; Simon 2019).
The debris of past accretions can be traced in the stellar haloes of galaxies, where gravitational accelerations are low and dynamical timescales are long.Substantial observational evidence of accreted substructures is present in our own MW, as our stellar halo is populated with a multitude of surviving dwarf galaxies (∼60 to date; McConnachie 2012) and disrupted satellites in the form of stellar streams (see Grillmair &Carlin 2016 andHelmi 2020 for recent reviews, and Mateu 2023 for the most up-to-date compilation of known stellar streams in the MW halo).
While these hierarchical processes are also expected to apply to ★ E-mail: jaclynjensen@uvic.casmaller scale systems, little observational evidence for stellar haloes around dwarf galaxies currently exists (Deason et al. 2022), and their observational study is complicated by their intrinsically low masses.They are particularly feeble systems: their shallow potential wells are sensitive to a multitude of internal (star formation, stellar feedback) and external processes (mergers, ram pressure stripping, tidal stripping) that can strongly influence their individual morphologies (e.g., see discussion in Higgs et al. 2021).Over a dwarf's lifetime, these environmental and internal factors make it difficult to disentangle the origin of the stars which are present in the dwarf's outskirts.
In addition, the intrinsic faintness of dwarfs means observations of the extended substructure in the dwarf's outskirts are challenging, especially when viewed against the dominant MW foreground populations.
With the advent of the Gaia satellite (Gaia Collaboration et al. 2016) and its wealth of kinematic data for ∼1.5 billion stars in the most recent data release (DR3; Gaia Collaboration et al. 2023), major progress has been made to study the individual resolved stars of MW dwarf galaxy satellites.A significant challenge when studying MW satellites is the ability to effectively and appreciably remove MW foreground contamination.However, Gaia proper motions provide one way of overcoming this challenge, as stars moving with the bulk proper motion of the dwarf galaxy can be separated from MW stars moving in the foreground.Recent works using Bayesian techniques (applied to Gaia photometry and astrometry) have found high success in obtaining dwarf member samples, thereby enabling the measurement of systemic proper motions for most of these dwarf galaxies (e.g., Pace & Li 2019, McConnachie & Venn 2020b,a, Pace et al. 2022, Qi et al. 2022, and Battaglia et al. 2022).
In some of these Bayesian framework examples, a primary assumption for the dwarf model is the spatial distribution of stars.In most cases, the stellar density of the dwarf is assumed to be an exponentially decaying function (e.g., McConnachie & Venn 2020b;Battaglia et al. 2022;Pace et al. 2022).However, as the stellar distribution in the outskirts of dwarf galaxies remains unknown, this presumption may unduly restrict the detectability of stars in the dwarf's stellar halo, should one exist.For example, satellite systems (as a consequence of their accretion into the MW's halo) are under the influence of the MW potential and susceptible to tidal forces, which can act to strip stars from the dwarf satellite.Additionally, recent simulations conducted by Goater et al. (2024) and Deason et al. (2022) showed that dwarf-dwarf mergers may exhibit a secondary, extended radial profile beyond the primary component of the dwarf galaxy.The extended secondary component in this case (depending on the mass ratio of the merger itself) may contain a large fraction of stars which originated from a past accretion event and may even exhibit lengthy stellar distributions similar to tidal tails.
Recently, an exciting observation was made by Chiti et al. (2021) who uncovered evidence of an extended stellar halo in the MW satellite, Tucana 2 (Tuc2).Spectroscopically confirmed members of this UFD were identified out to significantly large radial distances (∼9 circularized half-light radii, equating to a physical distance of 1 kpc).Given (i) the outskirt stars' lower metallicity compared to stars in the main body of the dwarf, (ii) a lack of a velocity gradient across the body of Tuc2, and (iii) the extended orientation of stars which are perpendicular to Tuc2's orbit, the authors argue that these radially distant stars are likely evidence of a past accretion event (i.e., dwarf-dwarf merger).
Motivated by this recent discovery, we initiated a search using Gaia to hunt for evidence of extended stellar structures around all of the known MW dwarf galaxies (and candidate systems).In this work, we expand the Bayesian methods of McConnachie & Venn (2020b, hereafter MV20a) to include extended multi-component substructure in the membership criteria, thus allowing us to locate members at much larger radii than in previous studies.This paper is organized as follows.Section 2 describes the relevant properties of the MW dwarfs we examine and our quality cuts in the Gaia dataset.In Section 3, we discuss the primary differences in our analysis method between this work and the implementation of MV20a.In Section 4, we apply the algorithm and estimate the contamination rates.Section 5 provides a discussion of the systems for which we find evidence of extended populations (i.e., "haloes") of stars.We conclude with a summary in Section 6.

DATA
We target all Milky Way (MW) dwarf galaxy satellites and likely candidates, as maintained in the tables originally published in Mc-Connachie (2012) 1 .However, we intentionally exclude the relatively massive Magellanic clouds, as well as the Sagittarius (Sgr) dwarf galaxy whose outskirts are so extended that they form the basis of an extensive literature by themselves (for the most recent Gaia-based work, see Gaia Collaboration et al. 2018b andKallivayalil et al. 2013 for the Magellanic systems, and Ibata et al. 2020 andRamos et al. 2022 for Sgr).Our dwarf inventory is therefore all other systems out to Leo T (heliocentric distance of ∼450 kpc), including a handful of recently discovered dwarfs.Newly discovered systems we include are: Boötes 5 (Smith et al. 2022, Cerny et al. 2022), DELVE J0155-6815 (Cerny et al. 2021a), Eridanus 4 (Cerny et al. 2021b), Leo Minor and Virgo 2 (Cerny et al. 2022), and Pegasus 4 (Cerny et al. 2023).
In Figure 1, we show the positions of all the candidate and confirmed dwarfs studied in this work in Galactic coordinates, and we list their names and positions in Table 1.
For each dwarf, we use the measured structural properties, distances, radial velocities, and metallicities as provided in the updated tables of McConnachie (2012), with some exceptions.For the newly added systems, we use parameters as reported in the relevant detection papers.We also updated parameters for Boötes 3 and Grus 1, both of which have recent structural measurements/updates (Moskowitz &Walker 2020 andCantu et al. 2021, respectively).The parameters used for the newly added and updated systems are reported in Table A1 in the Appendix.
While all dwarfs have measured distances and half-light radii, not all have measured ellipticities, position angles, radial velocities (RVs) and/or metallicities.We use all the aforementioned information for these systems, when available, to create the model likelihoods (discussed in Section 3) which represent the dwarf galaxy as a whole.For all dwarf parameters used in this work (structural parameters, positions, distance moduli, RVs, and metallicities) with the exception of updates/new dwarfs listed in Table A1, we refer the reader to Table 1 of MV20a.
Our analysis is based on the data available in Gaia early Data Release 3 (eDR3; Gaia Collaboration et al. 2021).Specifically, we use stellar sources with all five astrometric solutions (positions, parallaxes, and proper motions) and photometry (G, B  , and R  magnitudes).For most dwarfs, we select all such sources within 2 • of the center of the dwarf, which corresponds to on-sky distances well beyond 10 half-light radii for any dwarf.However, for larger systems (Antlia 2, Boötes 3, Carina, Crater 2, Fornax, Sculptor, Sextans 1, and Ursa Minor), we use a larger selection radius of 4 − 9 • to better capture any distant members.
To select stars with reliable photometry and astrometry in Gaia eDR3, we require that they pass the same criteria discussed in McConnachie & Venn (2020a, hereafter MV20b).Specifically, the sources must have high-quality astrometry (ruwe < 1.3), as given by Lindegren et al. (2018) and Lindegren et al. (2021a).We de-redden the photometry using the dust maps in Schlegel et al. (1998) and relevant extinction coefficients (see Equation 1 of Gaia Collaboration et al. 2018a).We then correct the parallaxes of the stellar sources by the global zero-point of −0.017 mas (Lindegren et al. 2021b).As part of the membership selection algorithm, we filter the parallaxes to only select stars in the field that are at broadly consistent distances to the dwarf (i.e., stars whose 3- parallax uncertainty places them within the 3- range of the dwarf, inferred from the distance modulus).
Following MV20a and MV20b, we do not require stellar radial velocities (RVs) to identify members.However, we later use these data to validate the method and estimate contamination fractions (Section 4.2) when RVs are available.The radial velocities used in this work are from: (i) the APOGEE Survey (Majewski et al. 2017) using the 17th data release of the Sloan Digital Sky Survey (SDSS; Blanton et al. 2017;Abdurro'uf et al. 2022).The APOGEE instruments (Wilson et al. 2019) are located on the 2.5-m SDSS and 2.5-m DuPont telescopes (Bowen & Vaughan 1973;Gunn et al. 2006) providing a dual-hemisphere coverage of Milky Way dwarf galaxies, whose targeting is described in Zasowski et al. (2013Zasowski et al. ( , 2017)), Beaton et al. (2021), andSantana et al. (2021).In this work we use the radial velocities from the APOGEE survey which are measured as described in Nidever et al. (2015) and Holtzman et al. (in prep).These exist for a limited subset of the dwarfs.The number of APOGEE stars within the Gaia field of each dwarf are reported in the sixth column of Table 1; (ii) a compilation of surveys targeting individual dwarf galaxies.The number of stars per Gaia field from these surveys are reported in the seventh column of Table 1.
In practice, for stars with more than one RV measurement, we take the weighted mean and weighted standard deviation of the multiple RV measurements.

METHODS
Our methodology is an extension of MV20a (inspired by Pace & Li 2019; see other similar implementations in Pace et al. 2022, Qi et al. 2022, and Battaglia et al. 2022), to which we direct the reader for specifics in the approach of the analysis.In summary, we assume every star in a given field is a member of one of two populations: (a) the satellite, or (b) the MW foreground.A star's total likelihood is therefore given as: where L  and L  are the likelihoods of the MW and dwarf models respectively, and   represents the fraction of stars in the dataset that are members of the dwarf.In our Bayesian framework, the two likelihoods in Equation 1 are based on three main components: (i) the stellar color and apparent magnitude (L    ); (ii) the radial distance of the star from the dwarf's center (L  ); (iii) the stellar proper motions (L   ), such that: The three likelihood functions L  , L    , and L  are constructed for both the satellite and MW foreground populations.The final probability of a star being a member of a dwarf is therefore: For each aforementioned likelihood function, we map each as a 2-D "look-up map".The method to create the look-up maps for the MW foreground models are the same as in MV20a.However, our methodology for the dwarf model varies from the original implementation.In the following subsections, we highlight the differences between the MV20a paper and our approach.All dwarf parameters that we use in the satellite model are adopted from Table 1 in MV20a, with the exceptions of the new and updated systems described in Section 2 and listed in Table A1.

Inclusion of the Horizontal Branch to L 𝐶 𝑀 𝐷
The implementation for the CMD satellite model is a likelihood function consisting of two separate likelihood maps: one for the red giant branch (RGB) and one for the horizontal branch (HB).The creation of the RGB likelihood map is performed in the same manner as in MV20a.We additionally include a model for the HB, which was not available in MV20a or MV20b.Our HB implementation is simplified compared to the RGB, since we had found that the synthetic populations that fit the RGB well do not describe the HB population to a sufficient accuracy.In particular, the synthetic HB does not generally extend as blue as necessary to match the data, nor does it reproduce the observed HB shape2 .We highlight this in the central panel of Figure 2, which shows the positions on the CMD for stars in the field of Sculptor and the PADOVA isochrone (Girardi et al. 2002; converted to Gaia photometric bands using Weiler 2018 corrections) that best matches the RGB.
For the HB likelihood function, we choose to represent the HB as a constant value in G-magnitude, ranging between the blue end set by (B  − R  ) = −0.25, and a red end extending to the RGB isochrone.This constant magnitude corresponds to the mean magnitude of an old (12 Gyr) and very metal-poor ([Fe/H] = −2.19)isochrone that is shifted by the dwarf's distance modulus.
In a similar manner to the implementation of the RGB, the likelihood function of the HB is essentially a Gaussian whose sigma is two factors added in quadrature.The first is an inherently assumed width of the stellar population; for the HB, we assume an approximate width of 0.1 mag in Gaia G.The second is the mean errors in (B  − R  ), to account for the errors in magnitude.
After the RGB and HB models are created, they are merged to create one likelihood map, such that each pixel in L    is the maximum likelihood from either HB or RGB model.Then, the map in synthetic stellar populations; nonetheless, these observed effects need to be accounted for.
is normalize such that the integral of the entire likelihood function (L    = L  + L   ) is equal to 1.
The left panel of Figure 2 shows an example of the constructed CMD likelihood for our example dwarf, Sculptor.On the right, we show the Sculptor stars whose probabilities obtained in the algorithm correspond to likely members (  ≥ 20%) in green.The HB model, though certainly simplistic, effectively identifies probable HB members of the dwarf.

Spatial likelihood Modifications to Identify Distant Members
For the dwarf spatial likelihood model, the original approach in MV20a is to create a look-up map, by co-adding many realizations of a 2-D exponential elliptical profile to approximate the dwarf's shape.Each realization is created by randomly sampling the structural parameters of the dwarf (position angle, ellipticity, and halflight radius) using Gaussian uncertainty distributions, or assuming a circular profile shape if not all structural information is available.In contrast, our approach assumes that there are two spatial components to every dwarf.The first (L , ) represents the main body, and corresponds to an exponential stellar density profile accounting for the structural parameters of the dwarf from the literature similar to the methods of MV20a.Here, we expand this expression to include a secondary profile (L , ) which represents a putative outer, extended component.We choose to model this as an additional exponential profile with two unknown parameters.The first is the normalization () and the second is its scale radius (  ).The overall spatial likelihood is therefore given by: where  is the distance of the star from the dwarf's center (accounting for position angle and ellipticity),   is the exponential scale radius of the dwarf (derived from the literature value of the half-light radius,  ℎ , such that  ℎ ≃ 1.68  ),  is the normalization of the outer component to the inner component (such that 0 ≤  < 1), and   is the exponential scale radius of the putative outer component (such that   >   ).If the preferred value for a dwarf is  = 0, this implies that no outer component is necessary in order to explain the data.In addition to requiring   >   , we set a weak prior on   such that it must be smaller than the radius of the field.
Our implementation of the spatial likelihood allows us to search over every dwarf to determine which dwarfs prefer an outer profile to exist in order to explain the data (i.e.,  > 0).However, we do not want to make strong assumptions about the shape of such an extra component, should it exist.Tidal tails can be very long and extended in specific directions (i.e., along the dwarf's orbit); alternatively, there may be an "excess" of stars at all position angles around the dwarf due to an extended halo.To approach the problem agnostically, we adopt two simple assumptions and run the algorithm using both in order to check for consistency.The first is that the projected shape (i.e., ellipticity) of the outer profile is the same as the main body of the dwarf; the second is that the outer component is spherical.For the first,  , =

√︃
2   + (   1− ) 2 where (   ,    ) are the tangent plane positions (, ) of a star, whose axes have been rotated by the position angle.For the circular outer profile,  ,  = √︁  2 +  2 .In contrast to MV20a, the spatial likelihood function for the dwarf now contains unknown parameters.This means we cannot easily con-struct look-up maps for the spatial likelihood as was done previously, but instead we need to calculate the spatial likelihood function direct from Equation 4for each realization.While this allows the easy addition of a second spatial component, it does not allow the inclusion of the uncertainties in the measured parameters for the inner profile.We comment further on the influence of the spatial uncertainties in Section 5.10.
In Figure 3, we show an example of the spatial likelihood implementation for the Sculptor dwarf galaxy using an elliptical outer profile.The left panel shows a 2-D tangent plane map color-coded by the total spatial likelihood (L  ) and the right represents the likelihood in 1-D as a function of elliptical half-light radius.The relevant scale radii of the two components are shown as dotted lines, where the scale radius of the inner profile (r  ) is shown in red and the best-fit scale radius of the outer profile (r  ) is shown in blue for comparison.In the 1-D panel on the right, we show the relative densities for the inner (red) and outer (blue) components, where the outer profile is scaled by the algorithm's estimate for .The associated scale radii are shown as vertical dotted lines in this panel.For dwarfs where the algorithm identifies an outer profile (i.e.,  is non-zero), we can also estimate where the outer profile begins to dominate over the inner profile.We refer to this boundary as the "transition radius" (   ), and for the Sculptor dwarf we highlight this boundary as a black dashed line in Figure 3. Stars beyond this boundary will primarily belong to the outer component, and may have a particularly interesting physical origin.
To confirm internal consistency when using our 2-component L  model to the previous work in MV20a and MV20b, we additionally compare our results to what we refer to as the "1-component" model.In practice, this means that the L  model is fixed per dwarf, as was done in the preceding works (i.e., errors in structural parameters are accounted for, and only one exponential is used to describe the dwarf's density).The difference here is our addition of the HB from Section 3.1 to the CMD model.We find it necessary to compare both the 1-component and 2-component models for internal consistency, and discuss this further in the following subsection.2).The inner half-light exponential radius (  , red dotted), the outer scale radius (  , blue dotted), and the transition radius designating the approximate radii where the fraction of starlight is dominated by the secondary outer profile (   , black dashed) are shown for comparison in both panels.The right plot shows the 1-D spatial likelihood for the inner (red line) and outer (blue line) spatial likelihood profiles.

Proper Motion Estimates and Added Prior
In the original implementation of MV20a and MV20b, the spatial likelihood model per dwarf was fixed, and so only three unknowns were to be solved for.These parameters being: the dwarf's systemic PMs in both directions, and the fraction of stars in the satellite.With the addition of an outer component whose unknowns describe the structure of the outer profile, our 2-component model now has five total variables to be estimated.We maintain the same requirement as the previous implementations to solve for the systemic PM of the dwarf, as it provides an additional method of verification in our methods.
For our 1-component runs, we first confirm that the addition of the HB in L    does not change the results of the proper motion estimate.In the Appendix Figure A1, we compare the relative measurements and 1- errors of our estimates in the 1-component model (black) to the reported PMs in MV20b (grey dashed) and Battaglia et al. (2022, blue square).For the more recently discovered systems mentioned in Section 2, we also plot the reported PM estimates from each dwarf's respective detection papers (pink triangle).In nearly all systems with previous measurements, we find consistency to past works within 1-to 2-.
In both versions of the 2-component runs (elliptical and circular outer profiles), we found that 51 out of the 64 satellites also have consistent PM estimates to MV20b.The 13 systems that fail the 2-component model do not produce consistent PM estimates due to an inability to converge in the spatial parameters  and   .We found that these particular systems have two properties in common: firstly, these dwarfs are the smallest UFDs such that the 1-component model only identifies ≲5 member stars within 6r ℎ .And secondly, these small systems have substantially higher MW contamination fractions within the same area, such that the fraction of members within 6r ℎ is ≲20%.This means that the MW contamination is significantly higher than the signal of the satellite within those radii (a MW fraction of ≳80%).Figure A2 in the Appendix shows this trend, with the systems that perform poorly highlighted in red text.
Given this issue, we decided to adopt a loose PM prior to better constrain L   in the case of the 2-component spatial model, such that the proper motion estimate must be within 5- of the PM of the dwarf derived in the single component case (accounting also for systematic errors; Lindegren et al. 2021a).Even with the addition of the PM prior, we find that these same 13 systems are unable to converge on any preferred values for the secondary component.These systems are: Boo4, Cet3, DES, Hor2, Ind1, LeoMi, Peg3, Pic2, Pisc2, Ret3, Tuc5, Vir1, and Vir2.It is unsurprising that the algorithm fails in these cases.We are essentially attempting to fit 5 parameters to less than 10 stars, and the algorithm is unable to resolve the dwarf's signal against the presence of the MW foreground.For the remainder of this work, we remove these systems from any further analysis.

A Note on Reported Parameters
Using the three likelihood functions (L    , L  , and L  ) for both the MW foreground and satellite models, we search over parameter space for the combination that maximize the total likelihood for a given dwarf.The estimated parameters are the following 5 variables:   * ,   , ,   , and   .These parameters are determined using the python package emcee (Foreman-Mackey et al. 2013).We apply the algorithm to all Gaia fields obtained for the dwarfs listed in Table 1.
Typically, the parameters for the outer profile have asymmetric posterior distributions.Therefore, we opt to report the values corresponding to the mode of the posterior distribution.The mode is determined by creating a histogram of the posterior using very fine, equally spaced bins.
The uncertainties on this estimate are then reported as the high and low bounds of the highest density interval (HDI).The HDI is defined as the shortest interval in a posterior distribution that contains a specified credible mass.Consequently, the mode is always contained within the credible region of the PDF.We elect to use a credible mass of 68%, which is essentially the Bayesian equivalent to  ± 1, or the 68% confidence interval in frequentist statistics.

RESULTS AND VALIDATION
The best-fit estimates for the spatial parameters ( and   ) allow us to identify a handful of dwarfs whose outer profile solutions are nonzero.That is, in the ∼60 dwarfs we explore in this work, we find a total of 9 which are suggested to host extended stellar substructure.In both of the 2-component runs, we find that the inventory of dwarfs where B is non-zero is the same between both assumed outer profile shapes.These dwarfs are: Boötes 1, Boötes 3, Draco 2, Grus 2, Sculptor, Tucana 2, Tucana 3, and Ursa Minor.The preferred parameter values for their outer stellar component in both the elliptical and circular runs are given in Table 2.
A natural consequence of the existence of an outer spatial profile is the transition at which the secondary component dominates the fraction of starlight.Largely, the location of this transition depends on the outer profile parameters,  and   .The location of this transition can be estimated empirically using the spatial likelihood function, and can then be used as a modest proxy to highlight stars that may be considered to be outskirt members.We reference this feature as the "transition boundary" (   ).In Table 2, we report the angular and physical semi-major axis of this transition, as well as the ellipticity of the transition boundary ellipse.In all systems, the position angle of this boundary is the same as that assumed for the dwarf itself (as in Table 1 of MV20a).
Prior to discussing the dwarfs which appear to host outer profiles in detail, we first explore the differences in the probability of membership between the two outer component shape assumptions in the following subsection.We then conduct an analysis to validate our methods using stars with RVs, when available.

Elliptical and Circular Outer Profiles
Our main rationale of this work is to determine, with high fidelity, which dwarf galaxy systems host outskirt stellar substructures.However, we also do not know a priori what the properties of the putative outer component may be.With respect to the projected shape of this component, our strategy has been to adopt two different shapes (one the same shape as the inner regions, i.e. elliptical, while the other is circular).
As reported in Table 2, we observe that the set of dwarfs the algorithm identifies as having two components is the same for both outer profile runs.The conclusion that these dwarfs have more stars in the outskirts than expected from a simple exponential, is apparently robust to our outer profile shape assumption.
In Figure 4, we show the differences in probabilities (  −   ) between these two datasets as a function of  ℎ .For clarity, stars with   ≥ 20% in at least one of these datasets are plotted.Those with   above this threshold in both datasets are shown as filled circles.For stars whose probability is ≥20% in one set but not the other, we plot them as black triangles (elliptical set only) and black + icons (circular set only).
Figure 4 shows that, for stars interior to the transition boundaries, the choice of an elliptical or circular outer profiles does not significantly change the probabilities.To emphasize this, each panel also shows a grey highlighted region, corresponding a change in probability of ±5% between the two runs.We find that >90% of the overall data presented in this figure are found within these bounds.
As seen at larger radii, the greatest probability differences occur at farther radial distances.We find that these differences are a result of the outer profile shape: for example, when examining the majority of the data which appear in only one set versus the other (all black icons), the positive y-axis consists only of the "elliptical only" set, while the negative y-axis are those of the "circular only" set.The differences seen in Figure 4, whereby some stars are more likely to be members in one realisation than the other, is a result of the stars' angle from the major axis; this concept is further corroborated when examining the tangent plane positions of these same stars in Figure 5.As shown in the figure, stars in these mutually exclusive sets are located along the dwarf's major axis ("elliptical only") and minor axis ("circular only").As a final confirmation, we further examined the residuals of the tangent plane spatial likelihoods (L , − L ,  ) and found the major and minor axes are enhanced in their respective model.
At this time, it is not obvious which model best fits these outer profile dwarfs as the structure will no doubt be linked to the stars' physical origin.However, our approach in this study is not to define the shape of dwarf galaxy outskirts explicitly, but instead is to identify whichever stars plausibly belong to these systems (of course, particularly those at large radii).Given our overall motivations, we remain agnostic about which shape should be utilized, and instead assume that any star with a reasonably high probability should be considered a member.For this reason, we decide to take the maximum probability (  ) per star, between all three runs (the single and the two 2-component versions − i.e, elliptical and circular outer profiles) to determine the possibility of membership.As evidenced by the data presented here, member stars may be overlooked when assuming one model over the other.In using the maximum probability of the 3 runs, we ensure all plausible stars are accounted for. .Over 90% of the data shown here have minimal changes in probabilities (±5%) between elliptical and circular models (highlighted by the grey region in each panel).The vertical magenta dotted and brown dashed lines are the transition radii solutions for the circular and elliptical profiles, respectively.Gru2 and Tuc3 are not included in the figure, since their structural parameters are assumed to be circular and there is no difference in probabilities between either 2-component run.
We anticipate spectroscopic follow-up will not only be necessary to confirm that a star is a member of its nearby satellite, but also to confirm their origin (accreted or tidal), and in turn derive the dwarf's stellar halo shape.

Examination of Purity Fractions and Contamination
In the Bayesian framework presented here, a typical method to classify data as a member of one population or the other is to select a decision boundary for the probabilities.Generally, a reasonable boundary is 50%; in our case, >50% implies a dwarf member, and <50% is of the MW population.In this work, we instead choose to adopt a boundary that is inferred from the most likely members, by confirming their membership via radial velocity information (when available), to better interpret a reasonable limit.This approach also allows us to implement a consistency check, by comparing the samples of true MW foreground and true spectroscopically confirmed members to their associated probability (i.e.,   ).
We do this by first calculating what we refer to as the "purity fraction".For each star in a given field with spectroscopic data, we determine whether its measured RV is consistent with being a member of the dwarf.A star is considered a true member if its individual RV is within 3- of the satellite's RV, accounting for the dwarf's velocity dispersion, the star's individual RV error, and the dwarf's systemic RV error in the calculation.RV, errors, and dispersion of the dwarf are taken directly from Table 1 of MV20a or from the updated systems listed in Table A1.The purity fraction is defined as the fraction of stars per probability interval that have RVs consistent with membership of the dwarf, specifically in each probability bin from 0 to 1 (0% to 100%).We note that ∼3 out of every 100 stars may have an RV more than 3- from the satellite, but still be a genuine member, although we will end up classifying it as a non-member based on its RV.
We highlight an example of this validation test as applied to the Sculptor dwarf galaxy in Figure 6.The RVs are taken as the weighted mean of APOGEE DR17 and individual survey measurements, if either exist for a given dwarf (listed in Table 1).For Sculptor, we plot all stars in the field from Gaia eDR3 with no known spectroscopic data (at the time of this publication) in grey in the tangent plane (left) and CMD (center) panels.Colored symbols show stars with RV data separated by probability of membership to Sculptor; stars with   > 10% (dwarf members) are shown as red points, while those with   < 10% (MW foreground) are marked as blue x icons.The right panel shows individual star velocities as a function of elliptical  4. We additionally plot in grey the 5 and 7 ℎ ellipses for each system.For direct comparison to the incidence of the Orphan/Chenab stream to Gru2, and the 300 km s −1 stream to Seg1, we include the approximate stream position (black line; from the python package galstreams; Mateu 2023) and highlight the approximate width ( = 1 • for Orphan/Chenab, and 0.4 • for the 300 km s −1 stream) in grey.
half-light radii (in units of  ℎ ), where the dwarf's systemic RV is highlighted as the black horizontal line.
The overwhelming fraction of stars flagged as members by this probability division (red points in Figure 6) exhibit RVs that are consistent with the dwarf satellite.This remains true, and is particularly apparent, for stars found out to large distances where the MW foreground population dominates the total number of stars in any radial annulus.We find this result to be particularly encouraging, as our kinematic model for the dwarf is independent of RV considerations.
We now estimate the level of contamination in our Sculptor sample by determining the purity fractions.Per probability bin (in steps of 0.1, or 10%), we determine the number of RV stars which are consistent to dwarf membership, divided by the total number (N) of RV stars in that bin.Errors for the purity fraction are reported as to account for Poisson statistics.Sculptor's purity fraction (and relative number of stars to obtain this statistic) are reported in Table 3.
Largely, the fractional purity for Sculptor across all   bins >10% is high (though, small number statistics in multiple bins suggest rather large errors), and a significant number of RV members are successfully identified with the highest probabilities (1731 member stars out of 1781, in the >90% bin).In the >90% bin, we note that there are 50 or so stars which are not members based on their RVs, yet their probabilities in the algorithm are high.We note these ∼50 stars are all located in the innermost regions of Sculptor (i.e., <2 r ℎ ; shown as the red points in the RV panel of Figure 6).Though they indeed have inconsistent RVs to Sculptor, we note that these particular stars are all a) towards the center of the dwarf, b) have consistent proper motions to Sculptor, and c) overlap the same region of the CMD as other Sculptor stars.Given that our model would therefore indicate that these stars are indeed members based on these three criteria, it is anticipated that they would be ascribed high probabilities (>95%, in fact).As evidenced by the purity fraction estimates in Table 3, these data represent the 3% of stars which contaminate the largest probability bin.We expect this type of contamination, and at close elliptical distances we do expect to have a larger number of stars with inconsistent RVs (though the fraction of contamination remains low).
We repeat this process for every dwarf which has a measured systemic RV.After obtaining the relative numbers of true members and the total per probability bin, we then choose to express the purity as a function of probability by marginalizing the data over all satellites.This "purity curve" (as it is hereafter titled) is essential to determining an appropriate decision boundary for dwarf membership based on the probability assigned, with the idea in mind to retain the most complete member sample out to large radii.
Figure 7 shows, for all stars with a radial velocity, the fraction of stars that are RV confirmed members as a function of   .The first panel shows this for all dwarfs (with measured systemic RVs), the middle panel for dwarfs with   ≲ −8 mag, and the right panel for dwarfs with   ≳ −8 mag.
In the idealized case, we might expect that the probability of membership is synonymous with the purity fraction, such that the correlation between purity and probability is one-to-one.However, all panels of Figure 7 show that the   estimated from the algorithm is in fact conservative.For example, we observe that ∼60% of stars with   = 20% are members (whereas, we would expect only 20% to be members).We find that the traditional limit of   = 50% actually results in membership fractions of >85%.
We also stress the unequal distribution of the number of stars in each   bin, resulting in vastly different errors between the intermediate bins versus the first and last ones (corresponding to the <10% and >90%, respectively; also seen in Table 3).This observation continues across all 3 dwarf galaxy groupings.It is first obvious that most stars in the field are not members of dwarfs, and so the lowest bin therefore contains a large number of stars (thereby resulting in small errors).Secondly, most spectroscopic member stars are correctly flagged in the algorithm as members with high   (such that the purity fraction is 0.97 for   > 0.9).For comparison, the raw number of RV stars whose probabilities are intermediate (i.e., 0.1 ≤   < 0.9) is ∼620 stars, while the first and last bins contain ∼3360 and ∼6670, respectively.Given the corresponding purity fractions in the first and last bins (0.09 and 0.97, respectively), we conclude that we appropriately identify the majority of MW and dwarf member stars in a given field, with very little contamination at high   .

DISCUSSION
In this section, we conduct a detailed review of each dwarf for which an outer profile has been detected via our methodology.For each dwarf, we search the literature for three main findings: (1) what is the most distant stellar member observed to date, and how far does each work trace each system; (2) do any previous works observe a break or transition in the stellar density profile, and if so, how do their results directly compare to ours; and (3) what evidence (if any) is there in the literature to argue for/against the dwarf's disruption via tidal influence imparted by the MW?
Throughout this discussion, we note whether previous works find evidence for a bump or break in the stellar density profiles of our 9 systems, and compare their findings to our estimated transition boundaries from Table 2.For easy comparison, we construct surface number density profiles (constructed by taking candidate members stars whose   ≥ 10%) presented in Figure 8.To highlight the benefit of the 2-component model and the detectability of radially distant stars, we also include an exponential function in each panel (red) for reference.
Additionally, we report the number "outskirt" candidate members; these are defined as stars whose radial distances place them farther than the transition boundary in Table 2, and are brighter than Gaia G ≤ 19.5 mag.The spread of members radially and in magnitude are shown in Figure 9, where each point is color-coded by   .For each dwarf, we also highlight the transition boundary (   ) as arrows to show the approximate "outskirt" region.
The following section contains common themes and questions for each dwarf.Complementing figures for this section are Figure 5, Figure 8, and Figure 9.

Boötes 1
Boo1 is one of the most luminous (  ≈ −6.Of these spectroscopic data in the literature, the most distant, confirmed members of Boo1 have been traced out to ∼4 ℎ (Longeard et al. 2022;Waller et al. 2023).Filion & Wyse (2021) also obtained photometric candidate members by probing the nearby blue horizontal branch (BHB) population, and found these BHBs occupy an extended outer envelope of the dwarf.For comparison, we detect 61 stars on the HB, and find the spatial distribution of these stars largely resembles that of Boo1 as presented in Filion & Wyse (2021) (see their Figure 6).The farthest BHB in their study was located at a radial distance of ∼10 ℎ , though this star does not have consistent proper motions to the dwarf and does not appear as a member in this work.We do confirm that all 29 spectroscopic members from Longeard et al. (2022) are correctly identified as members, and with high probabilities of membership (  > 50%).For comparison to previous works, the most distant member identified via our algorithm is at an elliptical radius of ∼9 ℎ from the center of Boo1, and does show consistent PMs to other spectroscopically confirmed members of Longeard et al. (2022).
As reported in Table 2, we find two possible solutions for the transition boundary in the elliptical (  , ) and circular (  ,  ) 2-component models.These values correspond to   , = 0.34 • (∼1.8 ℎ ) and   ,  = 0.57 • (∼3 ℎ ), respectively.Multiple works have previously reported transitions in Boo1's stellar components.Koposov et al. (2011) probed stars in the inner half-light radius and detected two separate kinematically "hot" and "cold" samples.The authors suggest the nature of these two samples are clearly segregated by different velocity dispersions and metallicities in each population.In a more extended spectroscopic study, Longeard et al. (2022) claim there is evidence of a break in the velocity dispersion and metallicity profiles.This break is located at ∼2 ℎ , which corresponds to our solution for   , .As evidenced by the stellar density plot in Figure 8, we also observe that our stellar density profiles break from the exponential at approximately this same position.
The tidal influence of Boo1 has been suspected in multiple previous works, such as Roderick et al. (2016) who observed a distinct overdensity near Boo1's estimated tidal radius (reported as 0.54 • and approximated via a King density profile).We note here that this particular report also coincides with our solution for   ,  .However, recent works in Battaglia et al. (2022) and Pace et al. (2022) have concluded that the orbit of Boo1 has as pericenter of ∼35 kpc, and N-body simulations from Read et al. (2006) have found that tidal shocking and stripping may be negligible for systems whose pericenters are >35 kpc.It is therefore questionable if tides are afflicting the morphology of Boo1.
At present, no work has confirmed if Boo1 is actively disrupting.Considering that (1) the spectroscopic data in Longeard et al. (2022) shows a moderate velocity gradient, (2) studies argue for a relatively small pericenter, and (3) the dwarf's elongation aligns with the direction of orbit (Pace et al. 2022), a natural conclusion is that tides could have affected the distribution of the outermost stars in Boo1.
However, there is some evidence favoring the alternative model (an extended stellar halo induced by a merger), or even a combination of the two scenarios.Firstly, Longeard et al. (2022) report a metallicity gradient for Boo1, which has only ever been observed in one other UFD (Tucana 2, see below).Pace et al. (2022) also argue that MW dwarfs with appreciably low central densities within the half-light radius, compared to the average density of the MW at the dwarf's pericenter, can be used as a proxy to identify likely disrupting dwarfs.The authors report that Boo1's density ratio is larger than that of known disrupting systems (Ant2, Boo3, Tuc3, and Sagittarius), placing it more in accordance with other typical MW satellites.And finally, abundances from Waller et al. (2023) showed that the outer stars (selected via the algorithm and shown to have consistent RVs to Boo1) have the same chemical properties as inner ones.This result may indicate that stars from the inner regions of Boo1 have been stripped to the outskirts, but additionally it cannot rule It is unclear at present if the outskirt members of Boo1 are due to the effect of tides, accreted members, or both.To conclude that Boo1's halo is more consistent with one scenario over the other would require ruling out the possibility of an extended halo from a contributing accretion, via evidence of a separate lower mass galaxy (e.g. a lower [/Fe] knee).Therefore, more spectroscopic observations with detailed abundances will be necessary to conclude the dynamics in this intriguing UFD.For future campaigns of outer profile stars, we identify 34 outskirt candidates whose magnitudes are G ≤ 19.5 mag and radial distance is greater than    .Only 9 of these stars have reported RVs in the literature.

Boötes 3
Boo3 was first identified in SDSS by Grillmair (2009); yet since its discovery, it has remained relatively unstudied.Compared to other UFDs, it is very typical in absolute magnitude (  ≈ −5.75 mag; Correnti et al. 2009) and luminosity (  ≈ 1.74 × 10 4 L ⊙ ).However, it is a particularly diffuse and large ( ℎ = 33.03′ = 0.56 • ; Moskowitz & Walker 2020) system.At a distance of ∼46 kpc (Carlin & Sand 2018), its physical half-light radius (∼450 pc) is relatively extensive for a UFD.Boo3's faint luminosity and diffuse nature in part explains the lack of substantial literature on this dwarf.
Since its discovery, Boo3 has been argued to be the progenitor (or at a minimum, associated with) the Styx stellar stream (Grillmair 2009;Carlin et al. 2009;Carlin & Sand 2018).A number of results support this conclusion: firstly, Styx overlaps Boo3 in the matched filter map of Grillmair (2009) at a distance corresponding to ∼45 kpc (i.e., the same distance of Boo3; Grillmair 2009).And secondly, Boo3 appears misshapen, double-lobed, and very diffuse (Grillmair 2009).These morphological characteristics suggest that the system is not in dynamic equilibrium and is being actively disrupted.
Currently, there are no publications that report any positions or kinematics of resolved stars in Styx.Only the distance to the stream and its length on the sky have been estimated, while orbit properties of Boo3 have been established and used for comparison to the stream's path.Grillmair (2009) 2022) also strongly supports a disrupted dwarf.Both studies find Boo3's orbit to be highly eccentric ( ≈ 0.85 and 0.95, respectively) and to have a small pericenter distance (7 − 9 kpc).It is important to note that only two MW dwarfs are known to have pericentric orbits of <10 kpc.Aside from Boo3, Tuc3's pericenter is also extremely close to Galactic Center, and is the only UFD known to be tidally disrupting (Li et al. 2018b) At such a moderate distance, it is surprising that no resolved Styx members have yet been identified and kinematically associated to Boo3, though the evidence appears to favor that Boo3 is likely a disrupting dwarf.The most distant spectroscopic member of Boo3 to date is located at 2.56 ℎ (Carlin et al. 2009;Carlin & Sand 2018), which we find falls within our estimated transition boundaries of   , = 3.02 ℎ and   ,  = 3.11 ℎ .We successfully identify all spectroscopic members reported in Carlin & Sand (2018), save for 2 stars whose updated proper motions in Gaia eDR3 suggest they are not consistent with membership to Boo3.Though the globular clusters M3 and NGC 5466 (and its stream; see Jensen et al. 2021) are close to Boo3 in projection, we find no contamination from these systems as their proper motions and CMDs are unequivocally distinct from Boo3 stars.
At present, the outskirts of Boo3 remain largely unexplored.In this work, we report 65 candidate stars (with no known spectroscopic follow-up) whose distances are greater than    and have magnitudes ≤19.5 mag.

Draco 2
Dra2 is a very small and nearby collection of stars whose estimated distance and brightness is ∼21 kpc and   ≈ 180 L ⊙ (  ≈ −0.8 mag; Longeard et al. 2018) respectively.With a half-light radius of a mere 3 ′ (0.05 • , or ∼20 pc) and a sufficiently low (marginally resolved) velocity dispersion (≲6 km s −1 at the 95% confidence interval), it has been a subject of controversy as to whether it is a true, dynamically cold UFD, or alternatively, a star cluster.
Though Dra2 is very nearby, spectroscopic follow-up of this intriguing system has been fairly limited due to a lack of bright (g < 19 mag) stars.First discovered in the Pan-STARRS 1 (PS1) 3 survey (Chambers et al. 2016) by Laevens et al. (2015), only two works have published further results on the system.Martin et al. (2016a) and Longeard et al. (2018) obtained a total of 14 spectroscopic members reported for Dra2.Of these stars bright enough to be included in Gaia eDR3, we report that we successfully identify all 12 spectroscopic members successfully via our 2-component algorithm.
The true nature of Dra2 is still largely speculated.Longeard et al. (2018) suggest that Dra2 seemingly conforms with other UFD trends in its absolute magnitude, metallicity, and larger size than most globular clusters, but ultimately it lies on the boundary between either object.Given the small number statistics, Longeard et al. (2018) are unable to deduce a metallicity dispersion for Dra2 which would yield further evidence towards a dwarf scenario.In addition, the marginally resolved velocity dispersion is very cold for a UFD, but is still comparable to other known dwarf systems like Boo1 (2.4 km s −1 ; Koposov et al. 2011) and Seg1 (3.7 km s −1 ; Simon et al. 2011).
If Dra2 is indeed a cluster, the approximate Jacobi tidal radius is only 10 pc (according to Longeard et al. 2018 and assuming the mass is solely stellar), suggesting that the stars we find at many half-light radii would be members that have migrated to significant distances as a result of tides.However, as shown in Figure 5, our sample of outer members are largely perpendicular to Dra2's elongation (also perpendicular to the direction of Dra2's orbit).This particular situation is contrary to the expected distribution of tidal tails, and therefore not likely to be initiated from a tidal stream.
However, multiple studies argue that Dra2 could be in the midst of tidal disruption.Namely, Pace et al. (2022) and Longeard et al. (2018) conclude that it is currently at pericenter (∼21 kpc), or close enough to experience tidal disturbance from the MW.Without a resolved velocity dispersion, Pace et al. (2022) note that the central density ratio argument (as suggested for Boo1) cannot be utilized for Dra2.However, they note that if the velocity dispersion is sufficiently low (∼1 km s −1 ), Dra2's morphology is likely caused by tidal disruption.
Currently, the most distant spectroscopic member of Dra2 is located at 2.6 ℎ (∼8 ′ ; Longeard et al. 2018).In our work presented here, we identify members up to 9.15 ℎ , and contrary to past works, we identify an additional, bright (G-band = 17.5 mag) RGB star at this large elliptical distance.This particular star was not identified in our 1-component runs, suggesting that the spatial likelihood without the free parameters of an outer profile was too restrictive to obtain a complete sample.Follow-up of this individual star may lead to a better understanding of Dra2, as it is located perpendicular to the system's orbit (as in the stellar features of Tuc2, see below).Further, if this most distant star is bound to the system, it could be argued that Dra2 (under the assumption of a dwarf galaxy nature) could have a sizeable DM halo.
Regardless, Dra2 lacks a substantial sample of bright enough RGBs for spectroscopic follow-up (save the one most radially distant example), as many would-be targets in the outskirts lie at G-band > 19 mag.In this work, we only find 4 stars brighter than G ≤ 19.5 mag, and only 2 of these are beyond the estimated    = 0.21 • (4.2 ℎ , or 80 pc).
As we show in Figure 5, the OC stream clearly overlaps Gru2 (grey highlight), but additionally the calculated transition boundary reaches the stream's center.Given that our method relies on proper motions, but not parallaxes/distances or RVs, it is very likely that the OC stream affects our determination of an outer profile for Gru2.It can also be noted that Gru2 is the only dwarf with no stars beyond our calculated    .Though Gru2's outer profile detection may be compromised by an unrelated structure in the field, we note that the algorithm identifies only member stars out to 5 ℎ (all of which are also members in our 1-component model).In comparison, previous work in Simon et al. (2020) estimated a tidal radius for Gru2 to be roughly 300 pc, or ∼3 ℎ .According to the central density to MW pericenter density ratio calculated in Pace et al. (2022), Gru2 is potentially tidally disrupting.The spectroscopic members traced in Simon et al. 2020 extend out to 2 ℎ , and so it is possible that some radially distant stars are evidence of this disruption.However, given the OC stream's width and overlap with these stars, it is certainly possible that the outer stars we identify are indeed members of the OC stream and not Gru2.
In dynamical models constructed in Koposov et al. (2022), the authors claim that, while Gru2 itself is not the progenitor of the OC stream, its orbit coincides with recent interactions with the OC progenitor.These close passages suggest that Gru2 may have been accreted with the OC progenitor, which is more similar in mass to a classical dwarf system.Determining the influence of tides on Grus2 is complicated by the interplay with not only the OC stream, but additionally the LMC.Orbit projections in Battaglia et al. (2022) suggest that, while Gru2 was not accreted with the LMC, it has been interacting/potentially captured by the Magellanic group within the past 200 Myrs.This is further corroborated by results in Santos-Santos et al. (2021) who reported that Gru2 is a potential associate of the LMC.Given the relatively close pericenters (24 − 31 kpc) determined by Pace et al. (2022) and Battaglia et al. (2022), and the dwarf's low central density, tides may have been important in shaping this system.Suffice to say, the complicated interplay of Gru2, the LMC, and the OC stream make for intriguing future work.

Sculptor
Of the many MW satellites, Scl was the first dwarf spheroidal (dSph) discovered (Shapley 1938), thereby earning its title as a classical dwarf.It is the fourth brightest dSph at   ≈ −10.82 mag (  ≈ 1.8 × 10 6 L ⊙ ; Muñoz et al. 2018), preceded only by Fornax, Leo1, and Sagittarius.It is a typical-sized system with a physical halflight radius of ∼310 pc ( ℎ = 12.33 ′ = 0.21 • ) and a large velocity dispersion (  ∼ 9.2 km s −1 ; Walker et al. 2009a) indicative of the presence of dark matter.
Given the brightness and angular size of this system, Scl's resolved stellar populations have long been studied, yielding a wealth of information and detail into this particular system.Multiple studies have explored the chemo-dynamics of Scl and concur that it consists largely of old (>10 Gyrs  et al. 2022).Much like Boo1 and UMi (see below), Scl's innermost regions (<2 ℎ ) are composed of a metal-rich and kinematically cold population that transitions to a metal-poor, kinematically hot population (e.g., Tolstoy et al. 2004;Battaglia et al. 2008).The transition from the inner to the outer population occurs at ∼0.2 • (∼0.95 ℎ ), and is potentially a result of outside-in star formation.
Regarding kinematic details, little evidence has been found to suggest that Scl is in the process of losing stellar mass via tides.Irwin & Hatzidimitriou 1995) and did not resolve an RV gradient (see also Martínez-García et al. 2022).At a heliocentric distance of 86 kpc (Muñoz et al. 2018), MW tidal influence is not purportedly large enough to strongly distort Scl's morphology.Recent simulations from Iorio et al. (2019) show that the stellar component of Scl is not obviously affected by tides, but up to 60% of the DM halo may have been stripped by present-day.
In comparison, recent work by Sestito et al. (2023a) utilizes our candidate members for Scl and obtained spectroscopic follow-up for the 2 most radially distant (out to 10 ℎ ) and relatively bright stars in our sample.The authors confirm that these distant stars are bona fide members of Scl based on their metallicities and RVs.As this finding expands the radial stellar distribution of Scl out to ∼3 kpc, the authors also explored the stellar density profile and its derivative − Γ(  ) − using the full sample of candidates in this present work and compared them to a tidal disruption simulation from Peñarrubia et al. (2008).With our candidate members, Sestito et al. (2023a) report a sharp deviation from an exponential profile in the stellar density in Scl (also shown in this work in Figure 8), located at about 25 ′ .According to the tidal model of Peñarrubia et al. (2008), tides can be responsible for an outer excess over the initial exponential profile, also called at a "kink" radius.This is followed by a region where Γ(  ) approaches a power-law with a slope of −4, before The "kink" in Γ(  ) is located at ∼25 ′ (∼2 ℎ ), corresponding with a profile excess also seen in the density profiles of Irwin & Hatzidimitriou (1995) and Figure 8.In our work, we find the transition boundary solutions are   , = 0.84 • (∼4 ℎ ) and   ,  = 1.06 • (∼5.1 ℎ ).These are much larger than the "kink" radius from Sestito et al. (2023a), yet are less than the nominal tidal radius from Irwin & Hatzidimitriou (1995, and less than the "break" radius from Sestito et al. 2023a which was found to be ∼80 ′ = 6.5 ℎ ).Seemingly, our proposed values do not correspond to these estimates.
Given the number of candidate members identified in this work, we also are able to produce an isophote contour map for Scl to highlight any morphological distortions in the outskirts (Figure 10).The contours are made with 2 ′ x2 ′ pixels, smoothed with a Gaussian kernel ( = 2 ′ ).Contour levels are set to [0.0025, 0.005, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.8] of the maximum counts.The central surface brightness corresponds to 29.8 mag/arcsec 2 , while the outermost contour corresponds to 23.3 mag/arcsec 2 .In agreement with past works, we find that Scl's morphology does not appear to be strongly distorted, but is rather nicely elliptical (particularly for the innermost contours).

Segue 1
At a heliocentric distance of 23 kpc (Belokurov et al. 2007a), Seg1 is one of the nearest dwarf galaxies.Its size is relatively compact ( ℎ = 3.95 ′ = 0.07 • , or ∼30 pc; Muñoz et al. 2018) and similar in luminosity to other MW UFDs (  ≈ −1.5 mag;   ≈ 340 L ⊙ ; Martin et al. 2008).In the discovery paper of Belokurov et al. (2007a), this system was proposed as a globular cluster containing very few stars.However, spectroscopic follow-up conducted by Geha et al. (2009) and Martinez et al. (2011) suggested that Seg1's internal kinematics are consistent with a dwarf galaxy.The former determined that Seg1 has an exceptionally large dark matter content (M/L ≈ 1340 − 2440 M ⊙ /L ⊙ ) while the latter claimed this system may be one of the highest DM density systems in the Local Group.
As the focus of many individual spectroscopic studies, Seg1 has consistently proven to be a particularly unique system.In addition to its high M/L ratio, Seg1 is one of the most metal-poor dwarfs in the MW ([Fe/H] = −2.7;Frebel et al. 2014;Norris et al. 2010), comparable to other UFDs like Tuc2, Hor1, and Boo2 (Simon 2019).Intriguingly, Frebel et al. (2014) observed a large metallicity dispersion in this system ( [/ ] = 0.95), yet all stars were found to have high -enhancement regardless of metallicity (see also Vargas et al. 2013).Given the high [/Fe] ratios across all measured stars in this study, the authors deduced that Seg1 lacks enrichment from Type Ia supernovae (SNe) which pollute the dwarf's ISM with iron over time, and that its star formation must have only lasted for a few hundred Myrs.Enriched only by the core-collapse of massive stars, Seg1's exceptionally distinctive chemical trends suggest it is a candidate "first galaxy" (Webster et al. 2016), or one that experienced a short burst of star formation, and has remained since unchanged.
Seg1 is located in a particularly advantageous position for spectroscopic study.Its close proximity, in addition to its high Galactic latitude (+50 • ) means that a) stars down to ∼1 mag below the main sequence turn-off are accessible for spectroscopic follow-up (Frebel et al. 2014), and b) contaminating MW foreground stars can be effectively removed.As such, previous work by Simon et al. (2011) obtained a fairly complete sample of radial velocity measurements for 98.2% of candidate member stars within 10 ′ (or ∼2.5 ℎ ) of Seg1.In an independent study, Norris et al. (2010) targeted stars out to ∼9 ℎ in the g-band range of [17, 21.4] mag.Additionally, Seg1's field was partially covered by APOGEE DR17.The relatively large number of stars with radial velocities in and around the dwarf (348 stars in total; see Table 1) alongside the high completeness achieved in previous studies for the inner regions of Seg1 provides us with a more complete dataset to test the robustness of our algorithm detection.
Using the samples from Norris et al. (2010) and Simon et al. (2011), we compare our probabilities estimated via the algorithm to their velocity-defined members.For the Simon et al. (2011) sample, of which 30 stars are bright enough to be in Gaia eDR3 and the most distant member is located at 2.9 ℎ , we obtain consistent membership probabilities for all stars.For the Norris et al. (2010) sample, we find similar consistency such that all but 2 of their spectroscopic members are also found as members in the algorithm.These 2 remaining cases were found in this work to have inconsistent PMs to the rest of Seg1, and thus have been flagged as non-members.
Thus far, the most radially distant spectroscopic member of this system is 5.25 ℎ .For comparison, our most distant member (with no present spectroscopic follow-up) is located at 6.8 ℎ .We determine the radius of the transition boundary is   , = 0.11 • (∼1.57 ℎ ) and   ,  = 0.14 • (∼2 ℎ ) in the elliptical and circular cases, respectively.For this well-studied system, we only identify 2 candidates beyond    with G ≤ 19.5.These are located at   = 5.5 and 6.3 ℎ , and G = 17.7 and 18.0 mag respectively.Though Seg1 has been extensively studied, it is at present unclear if it is currently undergoing tidal disruption.A photometric study conducted by Niederste-Ostholt et al. (2009) claimed to have found a 1 • extension of stars they deemed "tidal debris".However, a more detailed analysis in Geha et al. (2009) reports no kinematic tracers nor any other clear evidence of tidal disruption from their sample.Geha et al. (2009) also determine a Jacobi radii of 220 pc (or 7.4 ℎ ) where stars in Seg1 are likely to be unbound.We confirm in our work that no member stars (in Gaia eDR3) are beyond this radius.
Searching for tidal signatures in Seg1 is further complicated by the multitude of stellar streams intersecting in the field of the dwarf.These are: Slidr, Sagittarius, Gaia-10, Orphan/Chenab, and the 300 km s −1 stream (Mateu 2023; see Figure 5).In this less than fortunate circumstance, we do find sources of contamination from these substructures.6 stars with measured RVs are flagged as Seg1 members, but their origin is clearly the 300 km s −1 stream (at 18 kpc in heliocentric distances, and approximately the same proper motions; Frebel et al. 2013;Fu et al. 2018) while 8 more originate from a lower RV structure (likely Sagittarius, see Geha et al. 2009).However, we reiterate that the algorithm performs well at minimizing the number of false negatives, though there is contamination of false positives from other substructures (especially when they have similar proper motions to the dwarf) than would be ideal.Given the detailed spectroscopic follow-up thus far of the system, it is indeed impressive that only a few interlopers from these stellar streams managed to contaminate our sample.

Tucana 2
Upon first glance, Tuc2's characteristics are fairly typial compared to other UFDs.Its intrinsic brightness (  ≈ −3.9 mag;   ≈ 3150 L ⊙ ) and physical size ( ℎ = 9.83 ′ = 0.16 • , or ∼160 pc; Bechtol et al. 2015) places it squarely in the locus of other MW satellites in the size-luminosity relation (see Figure 2 of Simon 2019 and Figure 3 of Bechtol et al. 2015).Velocity dispersion measurements also support Tuc2's classification as a dwarf galaxy ( ∼ 3.8 km s −1 ; Chiti et al. 2023), and its mean metallicity falls well below the proposed "metallicity floor" for most MW globular clusters (Harris 1996(Harris , 2010 version; version;Wan et al. 2020 and references therein).A somewhat unique characteristic is its low metallicity, suggesting Tuc2 is one of the most metal-poor systems in the MW ([Fe/H] = −2.71dex; Ji et al. 2016b;Chiti et al. 2018), on average.Despite these parallels to other UFDs, Tuc2 has been shown to be a particularly interesting dwarf.In 2021, Chiti et al. searched the outskirts of Tuc2 for spectroscopic members and reported two compelling results.Firstly, they detected the most radially distant Tuc2 member, located at 1.1 kpc (10.71 ℎ in elliptical radii).This observation was (at the time) the most distant star in a UFD, preceded only by spectroscopic members of other dwarfs up to 4 ℎ (Frebel et al. 2016;Norris et al. 2010).Secondly, the authors identified 7 RGBs in Tuc2's outskirts, located at distances >2 ℎ .The authors suggested these stars were members of an extended stellar component which originated from the dwarf.As to the physical motivation for these stars to be displaced to such large radii, Chiti et al. (2021Chiti et al. ( , 2023) ) argue that their location conflicts with tidal stripping.Instead, Chiti et al. (2021Chiti et al. ( , 2023) ) suggest that this structure is a result of an early galactic merger, whose remnant accreted stars remain in the outskirts of Tuc2 (corroborated by the outside-in formation scenario in dwarfdwarf mergers; Benítez-Llambay et al. 2016).This initial discovery is significant to our understanding of dynamics of stars in low-mass systems, particularly as we do not yet have a census of dwarfs whose outskirts show evidence for this type of scenario.Chiti et al. (2021Chiti et al. ( , 2023) ) propose three main possibilities for the extended substructure in Tuc2.One hypothesis is that Tuc2's structure could be the result of bursty feedback or energetic SNe.By comparing the innermost stellar chemistries to outskirt members, Chiti et al. (2023) found that the metallicities and chemical abundances are similar to other UFDs.While UFDs in their formation may have experienced early bursty feedback (e.g., Wheeler et al. 2019), the conclusion that Tuc2 seems to have had similar enrichment compared to other UFD systems suggests that the energies of past SNe would make the extension in Tuc2 a much more common feature than currently observed.In agreement with this conclusion, the authors find that the chemical abundance trends favor a low mass progenitor, and therefore does not suggest particularly energetic SNe that could act to displace these stars.
The second (and most common conclusion) is a tidal scenario.However, the authors' observations find no conclusive evidence for tides.When modelling its tidal disruption, Chiti et al. (2021) determined that Tuc2's tidal debris would be observed following the direction of orbit, as is typical for many other dwarf and globular cluster streams.Even when including the time-varying potential of the LMC in their orbits, they find that Tuc2's observed stellar extension is indeed located perpendicular to the proposed direction of tidal debris.The authors further confirmed in later high-resolution spectroscopic studies that there is no clear RV gradient in Tuc2 that would be indicative of disruption (Chiti et al. 2023).
Also contrary to tides is Tuc2's current distance (58 kpc; Bechtol et al. 2015), but most importantly, its orbit places its pericenter between ∼35 kpc (Battaglia et al. 2022) and ∼45 kpc (Pace et al. 2022), depending on the model for MW mass (both including the LMC).In either case, this pericenter is not particularly small, and indeed Pace et al. (2022) do not find a particularly large density ratio for this system that would suggest it is actively tidally disrupting.
The lack of evidence for both the previous scenarios (i.e., tidal origins or expulsion due to energetic SNe) led Chiti et al. (2021Chiti et al. ( , 2023) ) to this final hypothesis: that Tuc2's morphology and steep metallicity gradient could be the remnants of an ancient merger.Simulations of a UFD (similar to Tuc2) in Tarumi et al. (2021) showed that the early (first ∼100s Myrs) merging of two galaxies can produce extended features in the direction of the merger, which could explain the perpendicular stellar extension in Tuc2.These simulations also showed that accreted stars from the merger can be deposited in the outskirts of the central dwarf, as the central dwarf is dynamically heated which forms an extended stellar distribution.Similarly, simulations by Benítez-Llambay et al. (2016) find that mergers in dwarfs can induce metallicity (and age) gradients in the surviving dwarf galaxy.Arguably, this could be why Chiti et al. (2023)  The nullification of the previous two hypothesis presents a strong case for the dwarf-dwarf merger remnant scenario, though direct observational evidence is clearly difficult to ascertain.Nonetheless, in our methods presented here we successfully identify all spectroscopic members in the literature as members of Tuc2.We find only one member star with a particularly low probability (7%); however, this star is evidently an AGB, and considering we do not actively model AGB in our methods, the CMD likelihood is very low for this star.
Given the expansive work conducted by previous studies, very few stars in our sample are candidates for follow-up.The most radially distant star we identify is the same as Chiti et al. (2021) at a distance of   = 10.71ℎ .Considering that the transition boundaries are   , = 0.70 • (∼4.38 ℎ ) and   ,  = 0.54 • (∼3.38 ℎ ), there are only 4 candidates in the outer profile regime which have not yet been observed.We note that 3 of these stars are HBs.

Tucana 3
Tuc3 is a fairly small ( ℎ = 6 ′ = 0.1 • , or ∼45 pc) and faint (  ≈ 800 L ⊙ ;   ≈ −2.4 mag; Drlica-Wagner et al. 2015) collection of stars located near the LMC.Discovered in DES alongside Gru2, the first observation of the dwarf was reported simultaneous to the ∼4 • -long tidal tails emanating from the dwarf's center.Besides the aforementioned Dra2 and Seg1 dwarf galaxies, Tuc3 is also one of the closest dwarfs (∼25 kpc).Its proximity has enabled multiple high-resolution spectroscopic studies such that its stellar stream and chemical trends have been studied in detail (e.g., Simon et al. 2017;Hansen et al. 2017;Li et al. 2018b;Marshall et al. 2019).
Though the presence of tidal tails is a strong indication of tidal disruption, the velocity dispersion of Tuc3 is smaller than typically associated with dwarf galaxies (an upper limit reported as   ∼ 0.1 km s −1 in Simon et al. 2017).However, Simon et al. (2017) argued that the low velocity dispersion could have resulted from the removal of stars which then acts to decrease the velocity dispersion in the system.Indeed, Marshall et al. (2019) observed a metallicity dispersion in Tuc3 − a result in favor of a dwarf galaxy origin since the potential wells in globular clusters are not typically large enough to retain the gas ejected from SNe.At present, its dwarf nature is not ultimately confirmed but certainly likely.
From these spectroscopic studies, it has also been found that Tuc3 is one of only two r-process enhanced UFDs in the MW (Hansen et al. 2017).Marshall et al. (2019) showed that this enhancement must have occurred before Tuc3's most recent pericentric passage, as these same enrichment features are observed in the core and most distant stars in the tails.
We identify most spectroscopic members; however, we find 7 members with substantially low probabilities that are not included in our member lists.Given their proper motions and position on the CMD, they indeed appear to be consistent with membership to Tuc3.We note that these 7 stars are the most radially distant of the sample, resulting in low spatial probabilities.This result is a limit of using a centrally concentrated outer component, rather than having a prior that would be more consistent to a longer stream structure.Nonetheless, we find that the most distant spectroscopic member within our probability limit is at 14.25 ℎ , though the actual most distant spectroscopic member to date is located at 16.63 ℎ .
The reported transition boundary for Tuc3 is    = 0.34 • = 3.4 ℎ .Though we know Tuc3 to be tidally disrupting, it is unclear if this boundary is physically motivated.Close inspection of the stellar density profile in Figure 8 of Drlica-Wagner et al. (2015) shows a break at this radius, which at minimum appears to be a transition from the Tuc3 core to the tidal tails.We note however that a break can also be seen in our stellar density profiles in Figure 8.Given our transition boundary, we identify 18 stars (5 of which are HBs) in the outskirts with reasonable magnitudes that have not been already spectroscopically observed.Pace et al. (2022) and Battaglia et al. (2022) updated kinematic information for Tuc3 and both reported very small pericenter distances (∼1 − 4 kpc), a very radial orbit, and high eccentricities (≈ 0.85 − 0.95).Pace et al. (2022) additionally note that Tuc3 is a system whose central density ratio is indicative of likely tidal influence.Interestingly, Tuc3 has also been proposed to be intertwined with the LMC system, as Erkal et al. (2018) reported that Tuc3 passed within tens of kpcs to the LMC some ∼75 Myrs ago.

Ursa Minor
The final system where we find evidence for an outer profile is the classical system, UMi.UMi resides at a relatively large heliocentric distance (∼76 kpc; Bellazzini et al. 2002), yet it was one of the first discovered MW satellites (Wilson 1955), owing largely to its luminosity,   ≈ −9.03 mag (  ≈ 3.6 × 10 5 L ⊙ ; Muñoz et al. 2018).It is a sufficiently large system ( ℎ = 17.32 ′ = 0.28 • or 380 pc; Muñoz et al. 2018) with a large velocity dispersion (8.6 km s −1 ; Spencer et al. 2018) similar to other classical dwarfs.
It has been speculated in many studies that tides may affect UMi; this was concluded from (i) morphological asymmetries, (ii) a relatively large ellipticity ( = 0.55; Muñoz et al. 2018), and (iii) stellar clumps that are offset from the dwarf's centroid, and appear aligned with UMi's orbital direction (e.g., Olszewski & Aaronson 1985;Irwin & Hatzidimitriou 1995).To better constrain the stellar density of UMi out to its nominal tidal radius, Palma et al. (2003) conducted a large-area (∼3 • radius) photometric survey of UMi, and confirmed two main results: the first being that UMi's isocontours appear Sshaped rather than elliptical (see also Irwin & Hatzidimitriou 1995), and the second was that UMi has two main peaks in stellar densities which are offset from the dwarf's center.Though the photometric data do not confirm if the most distant stars are bound or unbound to UMi, the authors propose their presence and UMi's morphology are likely as a result of tides.
Multiple properties of UMi that would be further indicators of tidal influence remain speculated.For example, UMi's most recent orbit estimates seemingly do not suggest a substantially close pericenter such that tides can strongly affect the dwarf (see simulations in Read et al. 2006).Pericenters estimated from Pace et al. (2022) and Battaglia et al. (2022) range from 35 − 55 kpc, depending on the potential used.Secondly, UMi's Jacobi tidal radius 3 estimates range 3 We note that the tidal radius assumes dynamical equilibrium.In the case of UMi, it has been suggested that the dwarf recently passed apocenter (see Sestito et al. 2023b and Figure 6 in Martínez-García et al. 2022) and may be reasonably relaxed such that the tidal radius could be useful.from ∼76 ′ (1.3 • , or ∼500 pc; Palma et al. 2003;Piatek et al. 2005) up to 4.4 • (or 6 kpc; Pace et al. 2020).We further note that Pace et al. (2022) do not find a substantially low central density ratio such that UMi is a tidal candidate, and though a velocity gradient in UMi has been detected in multiple studies (e.g., Hargreaves et al. 1994), the trend is identified along the minor axis, suggesting a modest internal rotation and not necessarily tides.
Spectroscopic study of UMi's outskirts is already in progress.Sestito et al. (2023b) obtained spectroscopic follow-up for 5 of the most distant and bright stellar candidates identified in this work, ranging from distances of 5.1 to 11.7 ℎ , confirming that all targets are indeed UMi members.According to previous studies, the most distant member was located at 6.8 ℎ ; with this work, we extend the outskirts of UMi out from ∼2.6 kpc out to ∼4.4 kpc.We find virtually no contamination in this sample, as we find only 4 contaminants in the RV sample (i.e., they show inconsistent RVs but are flagged as members in the algorithm).We note the contamination rate for UMi is only 4 stars out of the total 436 that have RV measurements.
In this work, we determined the transition boundaries in UMi are equal to   , = 1 • (or 3.46 ℎ ) and   ,  = 1.29 • (or 4.47 ℎ ).We note the circular solution coincides not only with a previous estimate for the tidal radius (from Piatek et al. 2005 andPalma et al. 2003 at 1.3 • ) but also to a break in the outskirts of the stellar density profile in Palma et al. (2003), located at ∼80 ′ (see their Figure 14, and our stellar density profile in Figure 8).In comparison, Sestito et al. (2023b) report the "kink" radius for UMi is ∼30 ′ (or 1.7 ℎ ) while the "break" is ∼225 ′ (or 13 ℎ ).
In the same manner as Scl, we also created an isophote contour map for UMi using our candidate members from this work.Figure 11 shows contour levels from [0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.8] of the maximum counts, where the central surface brightness corresponds to 30.8 mag/arcsec 2 and the outermost contour corresponds to 25.8 mag/arcsec 2 .Particularly in the central contours, there is some evidence of distortion (i.e., less clearly elliptical) moreso than in Figure 10 for Scl.
Aside from the 5 candidates in Sestito et al. (2023b), we find that there are 29 candidate outskirt(>    ) members with no current spectroscopic follow-up.

Additional Considerations
The algorithm has proven to be robust in identifying (i) candidate members of MW dwarf galaxies and (ii) systems that exhibit evidence of outer substructure, in agreement with much of the literature.In this section, we emphasize three main considerations: • The 9 systems identified in this work may not be the only systems in our sample with two components, but they are the only systems where the second component is identifiable using this technique with Gaia data.
• There are many systems with stars at relatively large radii (>5 ℎ ) that do not require two components to explain the overall stellar distribution.
• We have determined that the detection of an outer profile in our 9 systems is robust against uncertainties in the structural parameters (ellipticity, position angle, and half-light radius).We have confirmed that all systems use the most recent updates in structural parameters.
To address the first of these points, we note that our detection of an outer profile in each dwarf is limited by the density and physical extent of the structure.For example, it will be more difficult to detect a secondary component around small systems that contain only a handful of stars.As such, we cannot say there exist only 9 systems with extended stellar features.We anticipate the total number of dwarfs with outer profiles will likely increase with upcoming Gaia data releases and future large-sky surveys like LSST.
The second consideration relates to the fact that stars can exist at large radius from the host galaxy without requiring the galaxy to have a secondary component.For example, recent follow-up by Roederer et al. (2023) identified 5 member stars in Sextans at distances between 3.5 − 10 ℎ .In comparison, our 1-component model of Sextans does similarly identify distant candidate members, but only out to 6.8 ℎ .This is also the case in Fornax for which we find candidate members out to ∼7 ℎ (2.1 • ), or the same radial distance as stars reported in Yang et al. (2022).In these cases, a secondary component is not necessary to find member out to beyond 5 ℎ , and the existence of members at these radii is consistent with these systems having a single structural component.
And finally, we consider the limitations of the outer component as it relates to the assumed structural parameters.As discussed in Section 3.2, a consequence of the 2-component spatial model is that we cannot easily create a likelihood map that accounts for uncertainties in the dwarf's spatial parameters.A particular case arose with Gru1, as recent structural parameters for this system have been significantly changed.Specifically, Gru1's new  ℎ is in fact 2.5 times larger, and the position angle is now rotated into the next quadrant (i.e., it now is oriented West of North instead of East of North) resulting in a vastly different orientation and size in the likelihood model.In using the previous values for Gru1, we found that the algorithm determined that Gru1 had an outer profile, albeit the transition boundary was only slightly larger than the new measurement for Gru1's half-light radius.Additionally, the scale of the outer profile () was found to be much larger than in any of the other detections.
Determining the stability of the outer profile is an important outcome of this work, and so we thoroughly examined our list of outer profile systems to confirm that our detections are not hindered by large structural uncertainties.In doing so, we were able to verify that: (i) The uncertainties in position angle, ellipticity, and half-light radius are all sufficiently smaller than the changed percent error in Gru1 (or for position angle, a small enough change such that the orientation is not in the opposing quadrant); (ii) If updated structural information exists, we confirmed that using these new values do not produce a different membership list for these dwarfs, even if the individual probabilities may be different; and (iii) The updated values do not invalidate the detection of an outer profile.
Of the dwarfs where we detected outer profiles, the only systems with recently updated spatial parameters were Boo1 and Gru1.We found no change in Boo1's results when using the updated parameters in Longeard et al. (2022).

CONCLUSIONS & SUMMARY
In this work, we present an updated algorithm for the detection of dwarf galaxy stellar members using Gaia eDR3 photometry and astrometry.Our application differs from previous works in that we allow for the detection of an outer component in the stellar distribution (if there is one), enabling us to develop a census of dwarfs which likely host extended structure.We applied the algorithm to all known MW dwarf galaxy satellites (and dwarf candidates).Of these, we detect outer profiles in a total of 9 systems; these are: Boo1, Boo3, Dra2, Gru2, Seg1, Scl, Tuc2, Tuc3, and UMi.
We have shown that the presence of an extended profile in our model does successfully acquire additional members in the outskirts of these dwarfs.Many of these stars would have likely been missed without this change, as evidenced by the excess of stars seen in the stellar density profiles of Figure 8.The inclusion of these distant members indicates a divergence from the fiducial 1-component exponential function and at distances >    .The only exception is Gru2, which we discuss is limited by the presence of the OC stream that overlaps the dwarf in projection (and inconveniently has the same proper motions as Gru2).
For systems with an extended profile, we note that the addition of a secondary component does not change the vast majority of stellar members in the interior (<    ).Primarily, the change is evident in the outskirts where we observe that the 1-component case restricts the detectability of distant members.Using our membership lists for Scl and UMi, Sestito et al. (2023a,b) showed that the 1-and 2-component data produce exceptionally similar inner stellar density profiles, save for the evident truncation of stars in more distant radial bins.As the 1-component model has been shown to limit these detections, we report that our adjustments merely allow for the inclusion of additional stars at large radii, allowing us to probe further into the outskirts of these satellites.
Of the systems that we identify as having extended structure, a handful in the literature have already been reported as possibly disrupted systems (Boo1, Boo3, Tuc3, Scl, UMi; see references for each in Section 5).Encouragingly, we also find evidence of an extended feature in Tuc2.We therefore conclude that the algorithm performs well at locating systems with extended haloes and tidal disturbances.We again note however that we remain agnostic regarding the origin of these outer features, and that it will be up to future follow-up campaigns to make this conclusion.
An additional concept we report here is the boundary at which the spatial likelihood transitions from the inner exponential to the outer, which may prove to be an interesting feature.This feature corresponds to known breaks/bumps in the stellar density profiles of other works as explored in Section 5.As Figure 8 shows, the transition boundary nicely corresponds with the divergence of the stellar density from a singular exponential.In some clear cases (like Boo3 and Tuc3), we note the substantial excess of stars beyond this boundary appear as a stellar density bump.We argue it may be a useful limit to examine the central dwarf members compared to the stellar halo outskirt stars, for exploration in future works.
We also showed in this work that the algorithm performs exceptionally well at removing MW foreground contamination from our samples.By comparing our membership probabilities to the RVconfirmed dwarf members for which we have data, we confirmed that we obtain reasonable purity down to   ≥ 10% and determined that a purity of >85% can be obtained for stars with   ≥ 50%.However, in fields where the MW foreground contamination dominates significantly and the system in question is particularly small, we find that the algorithm is unable to effectively converge on spatial parameters for the outer profile (see Section 3.4).
This algorithm proves to be very efficient at detecting true members of dwarf galaxies, even those at large radial distances.Already, our source list of individual resolved members has proven useful for detailed study of these interesting dwarfs (Boo1, UMa1, and ComaBer in Waller et al. 2023;Boo5 in Smith et al. 2022;Ret2 in Hayes et al. 2023;Scl in Sestito et al. 2023a;and UMi in Sestito et al. 2023b).With the availability of deeper photometric campaigns (e.g., LSST) paired with the impressive precision and increasing accuracy of Gaia kinematics in each new data release, we anticipate the detectability of dwarf stellar haloes will be exceptionally promising.We look forward to disentangling the origins of these extended features, with the advent of the Gemini High-resolution Optical SpecTrograph (GHOST; Pazder et al. 2016) in the near future.

Figure 1 .
Figure 1.Galactic coordinates of the MW dwarf galaxy satellites used in this work, overlaid on the Gaia stellar density map (image credit: ESA/Gaia/DPAC).

Figure 2 .
Figure 2.An example of the CMD likelihood (L   ) calculated for the Sculptor dwarf galaxy.The left panel shows the CMD likelihood map, consisting of two models: one for the horizontal branch (L   ) and another for the red giant branch (L  ), as described in the text.The middle panel is the CMD of Sculptor field stars, with the corresponding PADOVA isochrone overlain in red.The RGB of this model approximates the data well; however it is clear that the HB does not extend as blueward, and misses many Sculptor HB stars.Stars in the right panel are those with   ≥ 20%, using our CMD likelihood map in the left panel, which clearly captures the full extent of the HB and exhibits little contamination.

Figure 3 .
Figure 3.An example of the spatial likelihood (L  ) for the Sculptor dwarf galaxy when assuming an elliptical outer profile.The 2-D spatial likelihood in the tangent plane is calculated using the values ( and   ) solved for by the algorithm (see values in Table2).The inner half-light exponential radius (  , red dotted), the outer scale radius (  , blue dotted), and the transition radius designating the approximate radii where the fraction of starlight is dominated by the secondary outer profile (   , black dashed) are shown for comparison in both panels.The right plot shows the 1-D spatial likelihood for the inner (red line) and outer (blue line) spatial likelihood profiles.

Figure 4 .
Figure 4. Differences between calculated probabilities in the elliptical and circular 2-component runs, for dwarfs with detectable outer profiles.In all panels, three subsets of the data are shown: stars which appear only in the elliptical dataset are shown as black triangles (  > 20%), stars which are only in the circular dataset are shown as black + icons (   > 20%), and stars in common between both profile runs are highlighted as colored points (both   and    > 20%).Over 90% of the data shown here have minimal changes in probabilities (±5%) between elliptical and circular models (highlighted by the grey region in each panel).The vertical magenta dotted and brown dashed lines are the transition radii solutions for the circular and elliptical profiles, respectively.Gru2 and Tuc3 are not included in the figure, since their structural parameters are assumed to be circular and there is no difference in probabilities between either 2-component run.

Figure 5 .
Figure 5. Tangent plane plots for all dwarfs with identified outer profiles.The relevant radii representing the transition boundary solutions (  ,  and   , ), and the three subsets of data (elliptical only, circular only, and stars in common) are represented the same as in Figure4.We additionally plot in grey the 5 and 7 ℎ ellipses for each system.For direct comparison to the incidence of the Orphan/Chenab stream to Gru2, and the 300 km s −1 stream to Seg1, we include the approximate stream position (black line; from the python package galstreams; Mateu 2023) and highlight the approximate width ( = 1 • for Orphan/Chenab, and 0.4 • for the 300 km s −1 stream) in grey.

Figure 6 .
Figure 6.An example of the radial velocity purity test, for the Sculptor dwarf.Grey points in the tangent plane (left) and CMD (middle) represent all field stars analyzed in this work.Stars with radial velocity information available are shown as colored points in all three panels; red indicates stars which we consider as dwarf members (  > 10%) are marked as points, while foreground MW contamination (  < 10%) are shown as faint x's.We show these radial velocities as a function of distance in the right panel, where Sculptor's systemic radial velocity is represented as the black horizontal line for comparison.

Figure 7 .
Figure 7.The fractional purity curves provided by the radial velocities of stars observed each dwarf, compiled together for both of the 2-component runs.These curves have been compiled for (a) all dwarfs, (b) classicals only, and (c) ultra-faints only.Green shading highlights the fractional purity regimes where statistically, at least 1 in 2 stars is correctly identified as a member (>50% fractional purity).
3 mag;   ≈ 2.9 × 10 4 L ⊙ ; Muñoz et al. 2018) and relatively nearby (66 kpc; Dall'Ora et al. 2006) UFDs.With a physical half-light radius of ∼200 pc (10 ′ , or 0.17 • ; Longeard et al. 2022), it is comparable in  ℎ to other MW classical systems.Boo1 has been the subject of many recent studies, resulting in a significant number of spectroscopic members (>70 stars; see compilation in Waller et al. 2023).While its size and chemical trends [/Fe]-plane suggest it is a more massive UFD, Frebel et al. (2016) determined Boo1's abundances of neutron-capture elements are much lower than observed in MW halo stars, a trait unique to UFDs (Frebel & Norris 2015).

Figure 8 .
Figure 8. Stellar density profiles for all dwarfs with identified outer profiles.Each plot is made with stars whose   ≥ 10%.Vertical lines represent important radii: grey is each dwarf's half-light radius, and magenta dotted and brown dashed represent the transition radii (  ,  and   , , respectively).For comparison, we additionally plot the 1-component exponential function in red to highlight the excess of dwarf members we detect using the 2-component model.

Figure 9 .
Figure 9. Magnitude versus distance for dwarf member stars with   > 10%.Arrows in each panel represent    to highlight the outskirts, and are color-coded the same as Figure 5.
first attempted to constrain the orbit of Styx (via an orbit-fitting method applied to the on-sky path), and was able to give first estimates for the PMs of Styx stars at the position of Boo3.This requires the assumption of either prograde or retrograde motion for the stream, but recent proper motion estimates of Boo3 in Pace et al. (2022) and Battaglia et al. (2022) would be consistent with Styx if the stream is in a retrograde orbit, providing further evidence that Boo3 may be the progenitor of Styx.The updated orbit of Boo3 provided by Pace et al. (2022) and Battaglia et al. ( . Boo3 allegedly recently passed pericenter (∼140 Myrs ago; Carlin & Sand 2018) and, according to Pace et al. 2022, does have a relatively low density ratio (on the order of Sagittarius and Tuc3), which they argue is in favor of a tidal disruption scenario.
) metal-poor ([Fe/H] = −1.8dex; Tolstoy et al. 2023) stars with well-mixed enrichment (e.g., Hill et al. 2019).Star formation in Scl is proposed to have had one epoch which ceased 8 − 10 Gyrs ago and lasted ≲1.4 Gyrs (Bettinelli et al. 2019); de los Reyes For example, both Pace et al. (2022) and Battaglia et al. (2022) find large pericentric distances for Scl, in the range of 45 − 65 kpc, depending largely on the assumed MW mass.Additionally, Tolstoy et al. (2023) recently consolidated archive spectra for Scl RGBs out to the nominal tidal radius (∼6.2 ℎ via a King density profile;
Figure A2.A comparison of the dwarf signal versus the MW foreground contamination, for systems with less than 20 high probability members (  > 50%) in the 1-component runs.We find that the addition of free-parameters in the spatial likelihood results in a failure to converge in the 2-component runs due to a) the small number of member stars to begin with (≲5), and b) the signal in the field dominates by a factor of >5.The systems where the 2-component algorithm fails are highlighted in red text.These dwarfs are excluded from the rest of the analysis.

Table 1 .
The Milky Way dwarf satellites explored in this work (1; marked with * if no systemic RV is reported for the dwarf) and their equatorial positions in the sky (2, 3).Galactic coordinates are included in (4, 5).We highlight the number of stars in a given dwarf field which also contain RV information from APOGEE DR17 (6), as well as the number of stars with RVs available in focused dwarf studies (7).Relevant references for (7) are listed in the final column (8).

Table 1 -
continued from previous page.

Table 2 .
Dwarfs with secondary outer profiles, as determined in both the elliptical and circular runs.The transition boundary (   , see text) is defined as the approximate boundary where the outer profile dominates the total fraction of starlight.We provide parameters describing this boundary by the semi-major axis and ellipticity of the best-fit ellipse.

Table 3 .
Purity fraction estimates for the Sculptor dwarf, as described in the text.Each probability bin is in increments of approximately 10%.Note that there are no stars in our sample in the 0.4 ≤   < 0.5 bin (i.e., no data with measured RVs are assigned probabilities in this range).