Open Access
Issue
A&A
Volume 690, October 2024
Article Number A166
Number of page(s) 20
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202450351
Published online 08 October 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The Milky Way’s (MW) stellar halo provides valuable information about the Galaxy’s history since it contains most of the metal-poor, [Fe/H] < −1, and oldest stars with an age > 10 Gyr (e.g. Beers & Christlieb 2005; Frebel & Norris 2015; Das et al. 2020). Although it is the most extended baryonic component of the Galaxy, reaching out to ≈200 kpc, its stellar content only contributes around 2% to the total stellar mass of the Galaxy (e.g. Deason et al. 2019). Moreover, the stellar halo is not smooth, but is comprised of substructures from globular clusters (e.g. Myeong et al. 2018; Massari et al. 2019; Forbes 2020; Limberg et al. 2022), to debris from past accretion events (e.g. Helmi et al. 1999; Koppelman et al. 2019; Naidu et al. 2020; Yuan et al. 2020; Limberg et al. 2021; Malhan et al. 2022), and satellite dwarf galaxies (e.g. McConnachie 2012). Throughout the paper, we refer to the Galactic stellar halo as the ‘halo’ unless otherwise explicitly stated when referring to the dark matter halo of the MW.

Determining the structure of the MW halo involves counting stars across the sky. With the current photometric surveys, this task is limited by dust extinction, the survey’s magnitude limit, and its footprint, which are encoded into the survey selection function. Knowing the survey selection function is crucial for a rigorous approach to measuring the radial density profile of the halo in order to interpret observations with theoretical models. For instance, a non-smooth halo would imply a recent merging activity whereas a smooth one would reveal a more quiescent formation history (Johnston et al. 2008; Deason et al. 2013). We expect both, with varying degrees of contributions, in the hierarchical formation scenario (e.g. Searle & Zinn 1978; White & Frenk 1991; Springel et al. 2006).

Despite the difficulties in obtaining a sample of stars that homogeneously covers the Galaxy, numerous studies mapped the halo density profile, albeit with a not-well-established consensus on the shape of the profile. The radial density profile has typically been described by a single power law, double power law, or Einasto profile (e.g. Deason et al. 2011; Xue et al. 2015; Hernitschek et al. 2018; Wu et al. 2022a) with the choice adopted partially being due to the sample and survey limitations. A halo described by a double power law also defines the radius at which the inner and outer slopes change, known as the break radius, rbr.

Deason et al. (2013) showed that a broken radial density profile is a signature of apocenter pileups of accreted dwarf galaxies. They studied the radial density profiles of idealised simulations with accreted satellites from Bullock & Johnston (2005) and argued that the presence of a break radius results from a massive merger event, that is, a merger that contributes the most to the total stellar halo mass. They concluded that the presence of a break radius in the MW stellar halo density profile is, therefore, a signature of a massive merger event (see e.g. Gilmore et al. 2002; Brook et al. 2003; Meza et al. 2005 for other signatures of the same merger). Recently, the accreted dwarf that dominates the MW stellar halo, and likely causes the observed broken profile, was identified as the Gaia-Sausage/Enceladus (GSE; Belokurov et al. 2018; Helmi et al. 2018; Haywood et al. 2018) and is estimated to constitute up to three-fourths of the inner (rgc < 30 kpc) stellar halo (e.g. Wu et al. 2022b, but see also Lane et al. 2023). Interestingly, Han et al. (2022) and Yang et al. (2022), using different methodologies, found evidence for a doubly broken halo density profile, and associated the two break radii with the apocenters of the GSE orbit (Naidu et al. 2021).

Since then, the GSE debris properties have been extensively studied taking advantage of the Gaia data combined with large spectroscopic surveys. Its stellar mass, peak [Fe/H], and time of merger range from 0.15 − 1 × 109 M (e.g. Feuillet et al. 2020; Naidu et al. 2020; Limberg et al. 2022; Lane et al. 2023), −1.45 < [Fe/H] < −1.17 (e.g. Amarante et al. 2020; Bird et al. 2021; Liu et al. 2022), and 8−11 Gyr (e.g. Gallart et al. 2019; Montalbán et al. 2021; Xiang & Rix 2022), respectively. The differences arise from distinct selection criteria and the samples used (see e.g. Carrillo et al. 2024 for a discussion). Despite the GSE merger event being estimated to have happened 10 Gyr ago, part of its debris is argued as being not completely phase-mixed. For instance, the Hercules-Aquila Cloud South and North (HAC-S, HAC-N, Belokurov et al. 2007a; Watkins et al. 2009; Simion et al. 2014) and Virgo OverDensity (VOD, Vivas et al. 2001; Newberg et al. 2002; Jurić et al. 2008) – overdensities observed in the stellar halo – are chemodynamically consistent with GSE stars (Simion et al. 2019; Perottoni et al. 2022). These substructures, except for HAC-N, have been argued to be a consequence of the GSE merger as seen in pure N-body models (Naidu et al. 2021). On the other hand, Donlon et al. (2020, 2024) argue that the presence of these substructures indicates a more recent time for the merger, ∼2−3 Gyr ago.

While the GSE debris is likely to dominate the stellar halo within the break radius, at larger distances, a more heterogeneous stellar halo, built from the debris of several other minor accretions, is found. The outer halo also comprises a myriad of known dwarf-satellite galaxies bound or in the process of being engulfed by the MW (e.g. Belokurov et al. 2007b; Koposov et al. 2015; Pace et al. 2022; Cerny et al. 2023). For instance, ongoing accretion events are observed for the Saggitarius dwarf (Sgr, e.g. Ibata et al. 1994; Johnston et al. 1995), and Small and Large Magellanic Clouds (SMC, LMC, e.g. Besla et al. 2007). Due to their relatively massive nature, pure N-body simulations of these two systems put into question the ‘Galaxy in equilibrium’ approximation. The stellar content and the dark matter halo of the MW are strongly perturbed due to the interaction with Sgr (e.g. Purcell et al. 2011; Gómez et al. 2013; Laporte et al. 2018a, 2019; Carr et al. 2022; Borbolato et al. 2024), and the LMC (e.g. Weinberg & Blitz 2006; Laporte et al. 2018b; Petersen & Peñarrubia 2020; Garavito-Camargo et al. 2019; Vasiliev 2023).

The LMC should induce a gravitational wake in the MW dark matter distribution and cause sloshing in the stellar halo in a classical Chandrasekhar (1943) sense from stars, and dark matter, being deflected behind the LMC’s orbit. Thus, the LMC’s orbit wake should create an overdensity in the stellar halo, known as the ‘transient wake’ (Weinberg 1998; Laporte et al. 2018b; Garavito-Camargo et al. 2019; Petersen & Peñarrubia 2020; Erkal et al. 2021; Garavito-Camargo et al. 2021). In addition to that, the movement of the MW-LMC system barycenter (Gómez et al. 2015) can also create a signature in the dark matter halo and consequently be imprinted on the stellar component, namely the ‘collective response’ (Weinberg 1989, 1998; Garavito-Camargo et al. 2019; Cunningham et al. 2020; Petersen & Peñarrubia 2021; Rozier et al. 2022). Numerical work also predicts an overdensity of LMC-stripped stars at the track of the LMC orbit (Petersen et al. 2022).

Recently the Pisces overdensity (Sesar et al. 2007; Watkins et al. 2009; Kollmeier et al. 2009; Nie et al. 2015) has been interpreted as a signature of the local wake on the stellar halo (Belokurov et al. 2019). This has been corroborated further with K giants by Conroy et al. (2021) which also suggested a detection of the collective response in star counts. To further complicate matters, retrograde debris has also been associated with GSE debris in the region close to Pisces (Chandra et al. 2023a).

Given these various contributions to the stellar halo, namely different accretion events and dynamical processes, a single-density profile would fall short of describing its complex structure. Ideally, by probing the whole sky with a homogeneous sample, one can measure any spatial dependence on the density profile. In fact, Hernitschek et al. (2018) has done it using RR-Lyrae from the Pan-STARRS1 survey (Sesar et al. 2017a) and fitted the stellar halo density profile for several bins in Galactic longitude, l and Galactic latitude, b, with 30° and 60° width, respectively. After removing known overdensities, they concluded that the stellar halo is well fit either by a single power law or an Einasto profile in all directions, which at face value would be considered at odds with a broken profile and with an out-of-equilibrium MW stellar halo.

Probing the stellar halo at large distances is contingent on measuring accurate stellar distances. Currently, this is mainly done with “standard candles”. These are stars with relatively well-defined distance-module relations, such as red giant branch (RGB), RR Lyrae, and blue horizontal branch (BHB) stars. In this work, we rely on BHB stars as they are in general associated with the MW stellar halo due to their old age (Dotter et al. 2010) and low metallicity (Santucci et al. 2015).

BHB stars have a relatively well-defined absolute magnitude range (e.g. Sirko et al. 2004) and they can be identified photometrically in the colour range, −0.3 < g − r < 0 (e.g. Yanny et al. 2000; Deason et al. 2011). Although spectroscopy is necessary to confirm BHB stars (e.g. Vickers et al. 2012; Barbosa et al. 2022), an appropriate colour selection can identify BHB stars with high confidence (e.g. Belokurov & Koposov 2016; Li et al. 2019; Starkenburg et al. 2019).

In this work, we use Legacy Survey data release 9 to probe the stellar halo density profile with BHB candidates. However, our measurement does not rely on identifying individual BHB stars. Similarly to Deason et al. (2011), we developed a methodology which takes into account the survey footprint, and for a given spatial selection, we measure the BHB fraction from the colour-colour distribution in a probabilistic fashion. This allows us to robustly estimate the halo density profile with its associated errors.

This paper is organised as follows: in Section 2 we present the Legacy Survey data selection and cleaning. In Section 3 we describe the methodology to calculate the BHB star fraction at a given spatial and magnitude selection, and how we model the radial density profile. We proceed by showing the density profiles for the MW, known overdensities regions, and several patches in the sky with high angular resolution, compared to previous studies, in Section 4. We discuss the implications of our findings for the structure of the MW stellar inner and outer halo in Section 5. Finally, in Section 6 we summarise our work and discuss future prospects.

2. Data

In this section, we briefly describe the Legacy Survey, the colour selection to discriminate BHB star candidates and possible sources of contamination.

2.1. Legacy Survey data selection

We use data from the Legacy Survey data release 9 (LSDR9, Dey et al. 2019). LSDR9 covers a total area of 14 000 deg2, in the region −18 ° < δ < +84° and galactic latitude |b|> 18°. The data contain measurements of the g, r, z photometric bands reaching a magnitude depth of approximately 24, 23.5, and 22.5 respectively. Its photometry coverage makes it an excellent database for photometrically identifying BHB star candidates.

We initially select objects within the extinction corrected colour1 intervals −0.4 < (g − r) < 0.2, −0.7 < r − z < 0.1, and within the apparent magnitude range 12 < g < 242. The first step of the data processing consisted of applying the LSDR9 bitmasks3 12 and 13 which remove pixels touching a large galaxy and globular cluster, respectively. This procedure removes objects within nearby galaxies from the Siena Galaxy Atlas 2020 (Moustakas et al. 2023) and stars close to globular clusters from LSDR9 catalogue. We also applied the g, r, z bands offset correction given by Zhou et al. (2023). They showed the presence of an offset between LSDR9 with Gaia synthetic photometry which is also strongly correlated with sky position. We stress the importance of this correction, as without it, any colour-based study can be severely biased.

2.2. Colour-colour space

Figure 1 shows the (g − r)×(r − z) plane of our initial sample. BHB stars are well known to lie within −0.3 < (g − r) < 0 (e.g. Deason et al. 2011; Belokurov & Koposov 2016) and this is our starting point to define the BHB fraction in a probabilistic fashion. Li et al. (2019) defined the boundaries for BHB stars using (g − r) and (r − z) colours, and here we re-write their relation defining the grz colour:

g r z = 1.07163 ( g r ) 5 1.42272 ( g r ) 4 + 0.69476 ( g r ) 3 0.12911 ( g r ) 2 + 0.66993 ( g r ) 0.11368 ( r z ) , $$ \begin{aligned} grz =&1.07163(g-r)^5 - 1.42272(g-r)^4\nonumber \\&+0.69476(g-r)^3 - 0.12911(g-r)^2\nonumber \\&+0.66993(g-r) - 0.11368 - (r-z), \end{aligned} $$(1)

thumbnail Fig. 1.

(g − r)×(r − z) diagram for the first selected LSDR9 objects. The blue dotted lines are the boundaries for selecting BHB star candidates following Li et al. (2019).

where BHB stars are scattered around grz ≈ 0.

As seen in Figure 1, the region associated with BHB stars, within the dashed blue lines, is within the limits of previous studies. This is reassuring as, although we use the same photometric band, the calibration of LSDR9 could be significantly different to offset the BHB relation. Compared to previous studies, we only slightly modified the range on (g − r) colour. Thus, throughout the paper, we only consider stars satisfying the following colour criteria:

0.14 < g r z < 0.07 0.30 < ( g r ) < 0.05 . $$ \begin{aligned}&-0.14 < grz < 0.07 \nonumber \\&-0.30 < (g-r) < -0.05. \end{aligned} $$(2)

2.3. Possible sources of contamination

The colour-colour selection defined in Section 2.2 can be contaminated by blue stragglers (BS), quasars, and white dwarfs (see e.g. Vickers et al. 2012). In the next section, we explicitly take into account the BS contamination, and here we apply some criteria to minimise the number of quasars and white dwarfs contaminants.

At faint magnitudes and towards bluer colours, there is an increasing number of quasars among all the observed objects. Moreover, quasars have a large scatter over the colour space due to their different spectral energy distributions and redshifts (see e.g. Ross et al. 2012). To remove potential quasar objects, we use the fluxes on the LSDR9 r-band and the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) photometric band W1, as the quasar’s dust emission makes them extremely bright in the r-band. We include the criterion −0.09 < FW1/Fr < 0.45 to select objects of interest, where FW1 and Fr are the observed fluxes in W1 and r, respectively. This selection is motivated by the fact that typical quasar objects have FW1/Fr ≳ 0.45.

To estimate the amount of contamination from white dwarfs, we used the Gaia EDR3 white dwarf catalogue (Gentile Fusillo et al. 2021) and verified that they have negligible contamination in the BHB region of this work. For instance, white dwarfs start to populate the selection box for g > 17 but mostly in the redder part of the grz colour. See Appendix A for details.

We also removed stars that lie within 5 times the effective radius of known dwarf galaxies from McConnachie (2012) catalogue, except for the Large and Small Magellanic Clouds (LMC, SMC) which we removed stars within an angular separation of 12.5° and 6.5° from their centres, respectively. Finally, we also mask the region with | B sgr | < 10.5 ° $ |\tilde{B}_{\mathrm{sgr}}| < 10.5\circ $, where B sgr $ \tilde{B}_{\mathrm{sgr}} $ is the Sgr stream coordinate system as defined in Belokurov et al. (2014).

3. Methodology

In the previous section, we defined the colour-colour selection, Equation (2), of the main sample which contains most of the BHB stars, and removed contamination from other sources, namely quasars and white dwarfs. We now proceed with the methodology to obtain the fraction of BHB stars.

3.1. BHB stars fraction

For any given selection within our sample, the BHB stars’ contribution to the grz distribution can be modelled as a Gaussian component with mean μ = 0 and dispersion, σ, due to the errors and scatter of the relation given by Equation (1). We illustrate this in Figure 2 where we show the grz histogram for stars within two different g intervals. In the first example, in the magnitude range 17 < g < 18, the grz colour distribution has a double peak shape, each associated with the BHB and BS components. Thus, we model the grz probability density distribution, P(grz), as a two-component Gaussian mixture with an additional exponential component representing other contaminant sources missed in the data cleaning of the previous section, that is:

thumbnail Fig. 2.

grz-colour distribution. The left panels show grz colour, see Equation (1), as a function (g − r) at different g ranges. The right panels show the P(grz), in grey, of the stars within the selection box on the left panel. The corresponding fit to the data following Equation (3) as well as the BHB, BS and background component are represented in green, blue, orange and magenta, respectively. These two examples are part of the exercise illustrated in the left column of Figure 3 which assumes flat priors for the free parameters. Throughout this paper, we do use Gaussian priors for the BHB stellar component when calculating the BHB stellar densities, see Section 3.1 for details.

P ( g r z ) = w 1 P ( g r z | BHB ) + w 2 P ( g r z | BS ) + w 3 P ( g r z | OTH ) = w 1 G ( g r z , μ 1 , σ 1 ) + w 2 G ( g r z , μ 2 , σ 2 ) + w 3 exp ( β g r z ) , $$ \begin{aligned} P(grz)&= { w}_{1} P(grz|\mathrm{BHB}) + { w}_2 P(grz|\mathrm{BS}) + { w}_3 P(grz|\mathrm{OTH})\nonumber \\&= { w}_{1} G(grz,\mu _{1}, \sigma _{1}) + { w}_2 G(grz,\mu _{2}, \sigma _2) + { w}_3 \exp (\beta grz), \end{aligned} $$(3)

where G(x, μi, σi) is the probability density function (PDF) of the Normal distribution for the BHB and BS components with mean μi and standard deviation σi. The background exponential is characterised by the parameter β. The contribution fraction of each component is given by wi and their sum equals 1. Therefore, there are a total of eight free parameters.

Selecting BHB star candidates using a probabilistic model based on the star’s photometry (e.g. Deason et al. 2014; Thomas et al. 2018; Starkenburg et al. 2019; Yu et al. 2024) or based on spectroscopic parameters (Sirko et al. 2004; Xue et al. 2011; Santucci et al. 2015; Barbosa et al. 2022) and model their spatial distribution has been extensively explored. However, these approaches fail to select BHB stars at the faint end of the magnitude range, where the photometric or spectroscopic uncertainties are large and the overlap between BS and BHB is significant. The methodology presented in this work lies in estimating the number (or fraction) of BHB stars at any given selection within our sample, instead of selecting individual objects based on their probabilities (see, e.g. Deason et al. 2011; Fukushima et al. 2019).

We estimate the number of BHB stars by fitting the distribution given by Equation (3) and sample the model with the Affine Invariant Markov chain Monte Carlo (MCMC) sampler, implemented in the emcee package (Foreman-Mackey et al. 2013). This allows us to obtain an error estimation from the posterior distribution of the eight free parameters of Equation (3). For each wbhb, i posterior, we calculate the number of BHB stars, Nbhb, i, as:

N bhb , i = w bhb , i N P , i , $$ \begin{aligned} N_{\mathrm{bhb},i} = { w}_{\mathrm{bhb},i}*N_{\mathrm{P},i}, \end{aligned} $$(4)

where NP is a random number drawn from a Poisson distribution with observed objects λ = Ntot, Ntot is the total number of objects within a given magnitude and/or spatial selection, and i is the posterior index. From the posterior distribution, we use the 50th percentile as the value of Nbhb and the difference between the 84th and 16th percentiles as its error4. We now proceed with a few tests to verify any dependence on colours and g-band for the BHB Gaussian component, (μbhb, σbhb), which could affect the modelling in the next section.

3.2. BHB star fraction test

We select stars with −0.30 < (g − r) <  − 0.05 and fit the grz distribution using Equation (3), including the errors of grz, for several intervals in g-band. Here we use flat priors for all the parameters. The right column of Figure 2 shows two examples of the fitting. While for the intervals shown the BHB stellar component is centred around grz ≈ 0, this is not the case for the fainter intervals, g > 19, where the flat priors fail to provide a good fit for the data. In fact, for the interval 15 < g < 19, μbhb, and σbhb are relatively constant and they only deviate for stars at the very bright and faint end. This is seen in Figure 3 which shows the results of the μ and σ for the BHB and BS components as well as their fractions for several g-band intervals. At the faint magnitude range, the errors in the grz colours are larger and thus strongly affect the BHB Gaussian distribution.

thumbnail Fig. 3.

Left column: the variation of the free parameters while fitting Equation (3) as a function of g assuming flat priors. While the horizontal lines correspond to the g interval on which the fit was performed, the vertical lines correspond to the 16th–84th interval of the posterior distribution. The BHB stellar parameters are relatively constant within 15 < g < 19 which motivated the use of Gaussian priors throughout this work, see Section 3.1 for discussion. Right column: now with Gaussian priors on the BHB component and within 14.5 < r < 21.5. The decrease in the BHB fraction as a function of magnitude is expected due to the intrinsic density profile of the MW halo.

Santucci et al. (2015) showed that BHB stars identified spectroscopically have a shift in the (g − r), being redder towards larger distances, due to an age difference of ∼2 Gyr. Motivated by this, we repeated the previous exercise for two colour intervals, −0.3 < (g − r) <  − 0.175 and −0.175 < (g − r) <  − 0.05, and verified that, for P(grz) given by Equation (3), μbhb and σbhb are consistent between the two intervals.

Based on these results, for the analysis described in the next sections, we complement the log-likelihood from Equation (3) by Gaussian priors on μbhb centred on 0.0005 with a width of 0.0005, and on σbhb, centred on 0.014 with a width of 0.001. We assume flat priors for μbs, σbs, β, and the fractions wi. Finally, in the next section, we included the criteria:

14.5 < r < 21.5 , $$ \begin{aligned} 14.5 < r < 21.5, \end{aligned} $$(5)

avoiding saturation problems of LSDR9 photometry and completeness in the bright and faint ends, respectively. In the right column of Figure 3 we show how the free parameters vary with the Gaussian prior and including the criteria in Equation (5). Now, the fraction of the BHB component decreases towards fainter magnitudes, as expected due to the density profile of the MW stellar halo (see similar exercise in e.g. Deason et al. 2011).

3.3. BHB stars density

We now proceed and describe how we calculate the radial density profile for any given volume within the LSDR9 footprint. We use the BHB absolute magnitude relation from Belokurov & Koposov (2016):

M g , bhb = 0.398 0.392 ( g r ) 0 + 2.729 ( g r ) 0 2 + 29.1128 ( g r ) 0 3 + 113.569 ( g r ) 0 4 , $$ \begin{aligned} \begin{aligned} M_{g,\mathrm{bhb}} =&0.398 - 0.392(g-r)_0 + 2.729(g-r)_0^2 \\&+29.1128(g-r)_0^3 + 113.569(g-r)_0^4, \end{aligned} \end{aligned} $$(6)

convert to heliocentric distance using the distance-modulus equation, and calculate the Cartesian Galactocentric coordinates, (x, y, z), assuming the Sun’s position (X, Y, Z) = (−8.2, 0, 0.027) kpc (McMillan 2017; Bennett & Bovy 2019).

The stellar halo of the MW is also known to be flattened (e.g. Zinn et al. 2014; Iorio et al. 2018; Hernitschek et al. 2018; Wu et al. 2022a; Yang et al. 2022), therefore we calculate the density profiles as a function of the flattened radius:

r q = x 2 + y 2 + ( z q ( r ) ) 2 , $$ \begin{aligned} r_q = \sqrt{x^2 + y^2 + \left(\frac{z}{q(r)}\right)^2}, \end{aligned} $$(7)

with q being the flattening parameter. Although it is consistently found that the inner halo is flattened while the outer halo is more spherical, there is no consensus on how flat the MW halo truly is and whether its flattening varies with Galactocentric radius or not. It is out of the scope of the current work to fit for the parameter q, thus we assume q = 0.77 for r < 30 kpc and q = 0.99 otherwise. This is motivated by results found in Hernitschek et al. (2018) and is in good agreement with the range of values found in the literature. We also verified that the results discussed below hold for different values of q that can be found in the literature.

A necessary ingredient for spatial density estimation is not only the estimate of the number of stars but also the volume covered by the survey. We can explicitly calculate the volume encompassed by a sample defined within a Galactocentric distance-module range, that is:

a < g 0 M g , bhb < b , $$ \begin{aligned} a < g_0 - M_{g,\mathrm{bhb}} < b, \end{aligned} $$(8)

where g0 is the extinction corrected apparent magnitude in the g-band and Mg, bhb is given by Equation (6). For the same volume, we can also estimate Nbhb and its error as described in Section 3.1 thus providing us with the estimate of the BHB number density, nbhb, and its error, nbhb, er, within the volume constrained by Equation (8) and the LSDR9 footprint.

3.4. Radial density profiles

The methodology described in the previous sections allows us to construct the radial density profile for the full LSDR9 footprint as well as for any direction or patch in the sky. Throughout this work, we fit the number density profiles with a double power law model:

n ( r q ) = { n r q α in , if r q < r br n r br α in ( r q r br ) α out , otherwise . $$ \begin{aligned} n(r_{\rm q}) = {\left\{ \begin{array}{ll} n_{\odot } r_{\rm q}^{\alpha _{\rm in}},&\mathrm{if}\; r_{\rm q} < r_{\rm br}\\ n_{\odot } r_{\rm br}^{\alpha _{\rm in}}(r_{\rm q}-r_{\rm br})^{\alpha _{\rm out}},&\mathrm{otherwise}. \end{array}\right.} \end{aligned} $$(9)

where n is the local density normalisation, rbr is the break radius, and αin/αout is the inner/outer slope.

In the next section, we construct the density profiles calculating the BHB density and its error estimated as in Section 3.1 at several bins in Galactocentric distance. Ideally, it is preferred to perform a fitting using all the available data points and not binned data, but our methodology relies on calculating the probabilistic fraction of BHB stars, and not identifying them individually. Therefore, we follow the approach of Anders et al. (2021) and estimate the free parameters of Equation (9) minimising the χ2 function using the binned profile and taking into account the error in nbhb. We use flat priors for all the four free parameters: n, rbr, αin, and αout, and emcee to sample from the posterior distributions. We take the median values as the parameter values, and the errors are from the 16th and 84th percentiles.

4. Results

The top left panel in Figure 4 shows the Mollweide projection in Galactic coordinates, (, b) of the stars satisfying the data selection criteria described in Section 2. We apply the methodology described in Section 3 to this data, and explore the Galactocentric radial density profile of the MW, as observed with LSDR9, and of regions with known identified substructures in Section 4.1, proceed to a higher angular resolution radial profile in Section 4.2, and depict the anisotropy of the stellar halo in Section 4.3.

thumbnail Fig. 4.

Top left panel: the Galactic Mollweide projection of LSDR9 stellar objects selected using Equation (2). The highlighted regions represent the location of known substructures: Sagittarius stream (pink), HAC-S (yellow), HAC-N (red), Pisces (green), VOD (blue), and OVO (cyan). Top right panel: the Galactocentric BHB stellar density profile of the whole LSDR9 footprint and without HAC-S, HAC-N, and VOD regions as the solid and dotted markers, respectively. Bottom panels: the Galactocentric radial density profiles for the Galactic regions associated with HAC-N, HAC-S, VOD, and Pisces. The crosses correspond to twenty and ten radial bins equally spaced in logarithm within 5 < rq/kpc < 120 for LSDR9 footprint and overdensities, respectively. We also show the resulting fit for 100 samples drawn from the posterior distribution as the solid lines.

4.1. The radial density profile of the MW and known overdensities

The top right panel in Figure 4 shows the Galactocentric radial density profile for the whole LSDR9 survey sky coverage. The BHB stars average profile has a clear break at r q = 19 . 1 1.8 + 1.7 $ r_{\mathrm{q}} = 19.1^{+1.7}_{-1.8} $ kpc with inner and outer slope ( α in , α out ) = ( 2 . 85 0.14 + 0.17 , 4 . 55 0.10 + 0.11 ) $ (\alpha_{\mathrm{in}}, \alpha_{\mathrm{out}}) = (-2.85^{+0.17}_{-0.14}, -4.55^{+0.11}_{-0.10}) $. The location of the break and the inner and outer slopes are within the range of previously identified broken power law profiles in the MW halo, see for example Table 6 in Medina et al. (2024) for a compilation of literature stellar halo density profiles. We defer the discussion and comparison with previous density profile estimations to Section 5.

The MW stellar halo harbours several substructures and overdensities that could potentially alter its density profile shape and some of them have been linked to previous merger events. We identify the regions associated with Hercules Aquila Cloud North & South (HAC-N & HAC-S), Virgo Overdensity (VO), and Pisces overdensity, whose coordinates are given in Table 1, in the top panel of Figure 4 as the coloured regions. We proceed and characterise the regions associated with the overdensities in the inner halo, rq < 30 kpc, HAC-N, HAC-S and VOD, and outer halo, Pisces. The second and third rows in Figure 4 show their radial density profile. The coloured markers correspond to the selected regions and the black solid line to the full LSDR9 profile, namely the average profile.

Table 1.

Galactic coordinates and heliocentric distance used in this work for Hecules-Aquila Cloud North and South (HAC-N and HAC-S, Belokurov et al. 2007a; Watkins et al. 2009; Sesar et al. 2010; Simion et al. 2014), Virgo OverDensity (VOD, Jurić et al. 2008; Bonaca et al. 2012), and Pisces (Nie et al. 2015; Belokurov et al. 2019).

HAC-N, HAC-S, and VOD have a broadly similar profile with a clear break radius, at the same position as the average profile, rq ≈ 20 kpc, and their slopes are within the errors as those estimated for the average sky. VOD shows an excess above the average density at all radii, while HAC-S and HAC-N only have excess within rq < 20 kpc, that is within the distance range previously reported in the literature. Moreover, HAC-N and HAC-S have a steeper outer slope compared to the average profile and VOD.

As HAC-N, HAC-S, and VOD have been associated with GSE, which likely dominates the stellar halo mass content, and they show a clear signature of a broken profile, one could argue that the global profile is broken due to the presence of these regions. We thus mask these three regions, as defined in Table 1, and recalculate the global density profile. However, after removing HAC-N, HAC-S, and VOD, the shape of the halo is still well represented by a broken power law, dotted markers in Figure 4 top right panel, at odds with Hernitschek et al. (2017) whose methodology favours a single power law after removing these substructures. We stress, however, that we probe the halo from 5 to 120 kpc, while Hernitschek et al. (2017) only for rq > 20 kpc. Nonetheless, the outer slope of our profile is in agreement with the slope these authors found.

Among the regions with substructures, the Pisces region is the only one that shows no evidence of a break in its density profile. This is seen both in Figure 4 but also in measurements of the inner and outer slopes, which are consistent within uncertainties, as well as by the fact that the break radius is not well constrained. Pisces overdensity was identified at large radii, and here we detect excess in its profile at rq > 50 kpc. At its direction, there is also an overdensity at rq ≈ 8 kpc not previously reported. Furthermore, the density profile in its direction is less steep than the average outer stellar halo density profile.

The behaviour of the density profile in four different regions hints at an anisotropic density profile of the MW stellar halo which could be “direction-dependent” rather than being described by a single density profile stratified on simple surfaces like ellipsoids. Motivated by this, we proceed to study the radial density profile for several regions in the sky as seen in the Galactocentric angular coordinate system.

4.2. The radial density profile at distinct Galactocentric angular coordinates

In this subsection, we show the radial density profile in several Galactocentric directions. We convert from Galactocentric Cartesian (x, y, z) to spherical coordinates and define the spherical angles Θgc and Φgc as the Galactocentric longitude and latitudes, respectively:

Θ gc = arctan z x 2 + y 2 , Φ gc = atan 2 ( y , x ) . $$ \begin{aligned} \Theta _{\rm gc} = \arctan {\frac{z}{\sqrt{x^2+y^2}}},&\Phi _{\rm gc} = \mathrm{atan2}(y,x). \end{aligned} $$(10)

We first show the density profile for the four quadrants of the Galactocentric projected sky in Figure 5. We do not remove the substructures located within each quadrant. While the four quadrants have signatures of a break and similar outer slopes, within the errors, they differ in the steepness of the inner slopes, αin. For instance, the profile of the second (Φgc > 0, Θgc > 0) and fourth (Φgc > 0, Θgc > 0) quadrants are not as steep as the first (Φgc < 0, Θgc > 0) and third (Φgc < 0, Θgc < 0) ones.

thumbnail Fig. 5.

Galactocentric radial density profiles for the four quadrants of the Galactocentric projected sky. The crosses correspond to ten radial bins equally spaced in logarithm within 5 < rq/kpc < 120. We also show the resulting fit for 100 samples drawn from the posterior distribution as the brown solid lines. The histogram corresponds to the posterior distribution of the break radius, rbr, and the dashed vertical lines are the corresponding 16th and 84th percentiles. The density profile of the full LSDR9 footprint, HAC-N, HAC-S, VOD, and Pisces are shown as the black, red, yellow, blue, and green lines, respectively.

There are also noticeable differences when comparing the profile of the regions associated with substructures to the Galactocentric quadrants in which they are located. HAC-N shows a distinct inner density profile, within the break radius, compared to the second quadrant, where it is located. This is shown in the left panel of Figure 5, where the brown markers represent the radial density for the quadrant and the red line is the best fit for HAC-N, as in the previous section. This quadrant’s profile closely follows the average profile (solid black line) in the inner parts, even without masking HAC-N when measuring it, in contrast to HAC-N. This reinforces the presence of HAC-N as a genuine overdensity, as this region neither follows the global profile nor that of the quadrant it is located.

The third quadrant, shown in the bottom left of Figure 5, includes two overdensities regions: HAC-S and Pisces, whose profiles are shown as the yellow and green lines respectively. While in the inner parts, the quadrant’s profile shows overall overdensity compared to the average one, its outskirts follows closely the average profile, in contrast to the Pisces region.

The first quadrant, top right panel of Figure 5, harbours the region associated with VOD. Its inner part follows closely the average profile, in contrast to the VOD profile (blue solid line), whereas the outer part also shows some overall underdensity. Finally, the fourth quadrant, which harbours the LMC and SMC but no known overdensity, differs only in the inner parts to the average profile having a less steep inner slope.

The differences and similarities between each quadrant, together with the regions associated with each overdensity, point to a directional dependence on the radial density profile which is generally not taken into account in previous studies. We investigate further the apparent asymmetry and pixelised the projected Galactocentric sky using the Python package healpy (Górski et al. 2005); healpy. We use nside = 2 which divides the sky into 48 patches (healpixels) each with an equal projected area of 859.44 deg2. We note that the LSDR9 footprint does not cover all the patches and a fraction of them are not fully covered. Nevertheless, the density for each region defined by a healpixel takes into account the LSDR9 footprint.

We focus on the regions which contain known overdensities with a coherent coverage in Galactocentric radius and show their radial density profile in Figure 6 (see Appendix B for the profile of all patches). We identify the patch location in the Mollweide projection of the Galactocentric coordinates in the top right corner of each panel and its “nested ID” number5. The purple crosses are the average density profile for the highlighted patch and the solid black line corresponds to the best fit for the LSDR9 footprint. For each region, we fit a double power law profile and show the interquartile of the posterior distribution of the break radius as the dashed vertical lines. We also explicitly highlight the locations associated with known substructures: HAC-N (red), HAC-S (yellow), Pisces (green), VOD/Outer Virgo Overdensity (OVO) (blue), and Sgr (pink). The less transparent the colour, the lower the contribution of the given substructure in a given pixel on the sky6. For Sgr we use the trail and leading arms distance-relation from Hernitschek et al. (2017, their Tables 4 and 5). We define OVO as a 20° region centred in right ascension and declination, RA = 207 ° ,Dec = −7°, respectively (Sesar et al. 2017b).

thumbnail Fig. 6.

Radial density profile for several patches in the Galactocentric sky covering an area of 859.44 deg2. The purple crosses correspond to the radial densities of the highlighted region shown in the Mollweide map of each panel. The purple solid lines correspond to the double power law fit using 100 samples drawn from the parameters’ posterior distribution, and the dashed vertical lines correspond to the 16th and 84th percentiles of the rbr posterior distribution. The black solid line is the best fit performed on the whole LSDR9 footprint. Each panel also shows the contribution of known substructures in a given region: HAC-N (red), VOD (blue), HAC-S (yellow), Pisces (green), and Sgr stream (pink). The darker the highlighted region is, the stronger the contribution of the substructure on the highlighted Galactocentric region.

Figure 6 shows the strength of the methodology developed in Section 3.1. By assigning a probabilistic BHB number for a given direction, we have a more detailed stellar halo profile compared to previous studies and can study any angular difference between the profiles. A similar approach was made by Hernitschek et al. (2018) as they probed the halo in several slices in (l, b). However, they sliced in patches which cover an area in the sky almost twice as large as the ones depicted here, removed HAC and VOD overdensities, and excluded the inner parts of the Galaxy, rq < 20 kpc. We now detail the features in those regions which include the known overdensities.

The top left and middle panels in Figure 6 show the density profile of two patches which include the HAC-N overdensity, as defined by the coordinates given in Table 1. The former contains most of the overdensity, as highlighted by the darker red band compared to the latter, and its profile shows a clear break radius at r < 20 kpc. At larger radii, the profile has a steeper decrease compared to the average LSDR9, but we note that the last two radial bins have a large uncertainty. On the other hand, the radial profile of the patch shown in the middle panel is consistent with a single power law and also shows an overall underdensity at radii greater than 20 kpc.

The top right panel in Figure 6 shows one of the sky patches that includes VOD. In this direction, the Sgr leading and trailing arms are expected to contribute at rq ≈ 35 kpc and rq ≈ 60 kpc, respectively, (as in Hernitschek et al. 2017). Although we have removed stars associated with Sgr (Section 2.3), there is still a slight overdensity at the distances associated with Sgr’s leading arm. Moreover, the overdensity associated with VOD seems to be at a slightly larger distance than what was previously reported, and closer to the break radius, namely at rq≈ 20 kpc.

The bottom row in Figure 6 shows the density profile in three patches of the Galactocentric sky which contain HAC-S and Pisces overdensities. In this coordinate system, they overlap in angular space, a feature which is missed in the Galactic coordinate system (see the top left panel in Figure 4). One can visualise it as being at the centre of the MW, looking towards a direction in the sky and seeing through two overdensities located at different radii. Firstly, the inner parts, located within the break radius, consistently show an overdensity associated with HAC-S. Secondly, the overdensity associated with Pisces seems to be mainly located in the patch centred at Φgc ≈ 67.5 deg and Θgc ≈ −41.8 deg7, and shown in the bottom right panel. This panel shows an overdensity through all radii, due to HAC-S and Pisces, and its outer slope is as steep as the average halo. These facts have implications on the interpretation of the possible origin scenarios of Pisces, as we discuss in Section 5.

4.3. The anisotropy of the halo density profile

The anisotropy of the halo density profile is depicted in Figure 7. The top and middle panels show the inner and outer slopes, respectively, of the double power law model fit for regions shown in Figure B.1. We only show the regions with more than four radial density bins available. The crosses indicate the patches where |αout|−|αin|< 1.5, the orange-coloured pixels in the bottom panel. These patches do not show strong evidence for a broken power law profile, or the least their difference is lower than what is observed in the global profile.

thumbnail Fig. 7.

Mollweide projection of the Galactocentric coordinates (Φgc, Θgc). The projection is divided into a grid of 48 pixels (nside = 2) with equal area. The blue, red, yellow, and orange regions delineate the regions associated with VOD, HAC-N, HAC-S, and Pisces at their typical distances. Top panel: the inner slope, αin, measured at rq = 15 kpc, of the two power law model. Middle panel: the outer slope, αout, measured at rq = 50 kpc, of the two power law model. Bottom panel: the difference between αin and αout. These panels highlight the anisotropy of the MW stellar halo profile, where several patches of the sky have distinct slopes compared to the average halo. We only show patches with more than four density points in the radial binning. The crosses and “x” indicate the regions with |αout|−|αin|< 1.5, and |αout|−|αin|> 1.5 with error in the break radius greater than 25%, respectively.

The inner parts of the halo show significant differences in the profile across the Galactocentric sky, as indicated by αin. For instance, the regions associated with HAC-N, HAC-S, and VOD have a flatter inner profile compared to the average LSDR9 profile, as indicated by the lighter-coloured pixels shown in the top panel, in contrast to most of the patches, except for regions in the fourth quadrant.

The regions marked with “x” indicate where |αout|−|αin|> 1.5, white- to purple-coloured pixels in the bottom panel, but the error in the break radius is larger than 25%. Interestingly, this is observed in the regions where HAC-S and Pisces regions overlap (as shown in Figure 6), whereas the patch which includes only Pisces does not show evidence of a broken profile. The other three patches with significant differences between inner and outer profiles, but with rbr poorly constrained, are either close to the Galactocentric plane, or near Sgr stream or LMC stellar halo contaminants.

Recently, Medina et al. (2024) showed how the density profile varies for four small patches, area ≲80 deg2, in the sky. They also highlighted the dependence of the halo density profile on the regions probed by a survey and the stellar tracer used. The results presented here do confirm significant differences in the stratified stellar halo, for rather larger sky patches, when compared to the global average. This has not been shown previously, as most studies either focused on a particular patch of the sky or did not probe the halo with a similar Galactocentric radial extent as done here.

5. Discussion

In this current work, we found that the average density profile of the MW halo, as traced by the BHB stars, is described by a double power law (Figure 4, top right panel) corroborating with a plethora of similar studies. Medina et al. (2024) argues that the apparent lack of consensus on the exact profile shape, such as either single- or double power-law, can be understood by the tracer being used, the magnitude range available, and the spatial selection used, such as the minimum/maximum Galactocentric radius in the density profile construction. In addition to that, and most importantly, the discrepancy found in the literature is also caused by the heavily anisotropic halo shown here.

Overall, the studies that find evidence for a single power law profile do not cover the inner parts of the Galaxy, namely rgc ≲ 20 kpc. Those which include the inner halo, in general, do find evidence for a double-power law profile with outer slopes with an overall agreement with the slope of the single-power law studies (see e.g. recent discussion in Yu et al. 2024; Medina et al. 2024 and their compiled list of previous works on the MW halo density profile).

In Section 4, we have shown evidence for a broken power law (BPL) density profile for the BHB stars in LSDR9 with a break radius at 20 kpc. The BPL profile is also observed when looking at the directions of HAC-N, HAC-S, and VOD overdensities and in some patches in the Galactocentric sky. Moreover, we do find significant anisotropy in the inner and outer halo as several patches have distinct slopes, see Figures 7 and B.2, where the latter also compares with slopes found in previous works. Moreover, some regions have a profile consistent with a single power law, for example in the direction of the Pisces overdensity.

We proceed and discuss how our results can shed light on the understanding of the structure of the stellar halo. We separate the discussion into two sections; in Section 5.1 we discuss the inner halo which is mainly associated with the GSE debris, rgc < 20 kpc, and in Section 5.2 we focus on the outer halo which hosts imprints of the LMC passage.

5.1. The inner halo and the GSE debris

Deason et al. (2013) made the connection between the presence of a broken profile and the merger history of the MW using Bullock & Johnston (2005) simulations. The break radius corresponds to the apocenter pile-up of the debris from a past massive merger. The inner part of the halo, within the break radius, is dominated by the debris of the main merger and thus has a less steep density profile compared to the outer parts, dominated by minor accretion events. As they observed this in the MW, they concluded the MW had gone through a main accretion event in its early stages. This was yet another strong evidence among others (e.g. Chiba & Beers 2000; Gilmore et al. 2002; Meza et al. 2005), of the GSE merger event later confirmed with Gaia data (Belokurov et al. 2018; Helmi et al. 2018; Haywood et al. 2018).

More recently and thanks to Gaia and large spectroscopic surveys, the halo density profile is also built with halo star candidates selected based on stellar chemistry and kinematics (Han et al. 2022; Yang et al. 2022; Lane et al. 2023). This methodology relies on the fact that the GSE debris is the major contributor to the stellar halo mass (e.g. Naidu et al. 2020). Intriguingly, this approach has also provided evidence for a doubly broken halo, namely a three-power law model, when selecting GSE star candidates chemically (Han et al. 2022) or with a metallicity-kinematic decomposition (Yang et al. 2022). Each break radii is associated with the apocenters of GSE orbit during its merger (Naidu et al. 2021).

In addition to that, the GSE debris has been associated with HAC-N, HAC-S, and VOD substructures. They have been shown to share similar dynamics Simion et al. (2019) and chemodynamical (Perottoni et al. 2022) signatures as the GSE near the Solar neighbourhood. In Section 4.1, the radial density profile of these regions is shown for the first time for a large distance extension, 5 < rq/kpc < 120, which allows showing clear evidence for the presence of a broken radial profile. These regions share a similar rbr ranging from ∼22 − 24 kpc, similar inner and outer slopes, and a clear signature of the overdensity in the inner parts, where these substructures have been identified. The location of their break radius is also ∼1 − 2σ larger than the global average. The break radius of the halo profile has been associated with the apocenter of the GSE dwarf orbit (e.g. Naidu et al. 2022; Han et al. 2022). Thus, the difference between the break radius of the overdensities and the global average is unexpected in the scenario where the GSE dominates the inner halo and these overdensities are part of its debris. We also show in Figure 5 that these overdensities’ profiles are also consistent with the Galactocentric quadrants where they are located. In Section 4.2 we further explored these regions in smaller patches in the Galactocentric sky reinforcing the signature of the overdensities.

The density profile for these specific overdensity regions associated with the GSE merger serves as a constraint to idealised simulations that mimic this merger event. For instance, Naidu et al. (2021) presented an N-body model, without gas and star formation, which reproduces several properties of the observed GSE debris, including HAC-S and VOD overdensities, but does not reproduce the substructure in the HAC-N region. The latter is also an important constraint, and its overdensity is also observed in previous studies (e.g. Huang & Koposov 2022), due to its similarity with HAC-S and VOD. In simulations with gas and star formation, either tailored (Amarante et al. 2022) or in fully cosmological context (Orkney et al. 2023), the debris of GSE-like mergers fully phase mixes in the inner parts after a few Gyr and do not show evidence for clear overdensities like HAC-N, HAC-S, and VOD.

To complicate matters further, Orkney et al. (2023) has also shown evidence that MW-like galaxies in the Auriga simulations (Grand et al. 2017) can have stellar halos with high radial velocity anisotropy, that is like the GSE component in the MW, being built with up to three main dwarf satellites of similar mass. These accreted stars have enough similarity in chemistry and dynamical space making it almost impossible to distinguish them apart solely on [Fe/H], [α/Fe], and with dynamics. Moreover, Donlon & Newberg (2023) discusses the possibility of multiple dwarfs being able to explain the chemodynamics of the inner halo, where the GSE is claimed to be the main contributor. We note, however, that in situ contamination can contaminate halo-selected samples (da Silva & Smiljanic 2023).

It remains to be seen, whether the presence of the HAC-N, HAC-S, and VOD can be indeed a feature of a unique merger event which happened 10 Gyr ago, or if any of these overdensities are indeed signatures of a more recent merger as argued by Donlon et al. (2020, 2024).

5.2. The outer halo and the LMC wake

To have a general picture of the outer halo, we divided the Galactocentric sphere into 192 equal-sized area pixels, using healpy Python package with HEALPIX pixelisation nside = 4, and calculated the halo density at 50 < rq/kpc < 120 using the methodology described in Section 3.1. We then calculate the density contrast and show it in Figure 8. The LMC past trajectory is shown as the purple line (Vasiliev 2023)8, with the thicker part corresponding to the orbit within the distance range shown in the map.

thumbnail Fig. 8.

Density contrast of MW BHB stars in the stellar halo at 50 < rq/kpc < 120 seen in Galactocentric coordinates (ϕgc, θgc). The Mollweide projection is divided into a grid of 192 pixels (nside = 4) with equal area. We delineate the approximate Galactocentric regions associated with OVO, Pisces, and Northern OverDensity as the cyan, green, and magenta lines, respectively. The present-day position of the LMC and SMC are shown as the purple scatter points. LMC’s past trajectory (see Vasiliev 2023 for details) and Sgr stream are represented by the solid purple and dashed black lines, respectively.

Although we removed stars within 12.5° from the LMC centre, there is still an overdensity seen in the pixels surrounding LMC’s current position. These are likely contaminants of LMC star members, as the LMC stellar halo can extend up to 30° from the LMC as seen with BHB stars (Belokurov & Koposov 2016) and RR Lyrae (Petersen et al. 2022).

It is also evident from Figure 8 the presence of several regions with a density contrast close to −1. Due to the steep density profile at these distances, αout = −4.55, the average density is dominated by the regions with overdensities. In fact, ∼55% of the BHB stars, as estimated with our methodology, is within ∼20% of the pixels. If the BHB stars followed a purely Poisson distribution, the same amount of stars would be within ∼40% of the pixels. This reinforces the highly anisotropic stellar halo of the MW being dominated by a few patches with overdensities.

The Pisces overdensity, as highlighted approximately by the green circle in Figure 8, has been recently associated with the transient dark matter wake caused by the LMC passage. Its signature has been associated with an overdensity of RR Lyrae (Belokurov et al. 2019), K-giant stars (Conroy et al. 2021), main-sequence stars (Suzuki et al. 2024), and is now also observed with the probabilistic BHB star mapping developed in this work. This transient response of MW’s (dark) halo to the LMC is akin to the Chandrasekhar dynamical friction wake (Chandrasekhar 1943) and traces the orbital path of the LMC in the halo at large distances (e.g. Garavito-Camargo et al. 2019). The shape of this transient wake signal should also be proportional to the mass and velocity of the LMC (and the DM density of the ambient medium). Moreover, Rozier et al. (2022) showed that the detailed geometry of the local wake is also sensitive to the underlying anisotropy, tangential or radial, of the tracers. Our BHB stars map can thus be used to measure the extent of the transient wake response. To this end, we study the BHB star density variations along the sky patches crossed by the LMC’s past orbit from Vasiliev (2023) for 50 < r/kpc < 120.

The top and bottom panels of Figure 9 show how the density varies along the Magellanic Stream Coordinates9 (Nidever et al. 2008). For this exercise, we coarsened the angular resolution, compared to Figure 8, to reduce the Poisson noise. Despite the relatively large patches, it clearly shows the density peak at the Pisces region, point B1, and how it decreases going further away in LMS and BMS. From Figure 9, we can estimate the upper limit of the Pisces overdensity extension to be 60° along and 30° across the LMC’s orbit, which would correspond to a wake width of ∼32 kpc at ∼70 kpc. Future density profile maps along the orbital path of the LMC, as those shown in Figure 9, to characterise the transient wake signal will help further constrain models of the MW-LMC interaction and provide an important test for the nature of dark matter on Galactic scales.

thumbnail Fig. 9.

BHB star density of the sky patches crossed by the LMC’s past orbit (Vasiliev 2023) as a function of the (LMS, BMS) coordinates (Nidever et al. 2008). Each patch has an equal projected area of 859.44 deg2 (healpy using nside = 2) and the density is calculated within 50 < r/kpc < 120. The 16th–84th percentiles of the average BHB density within this distance is shown by the shaded green line. Top panel: the density profile along LMC’s orbit from point A to the pixel containing LMC. At B1, the pixel dominated by Pisces overdensity shows a clear signature of the overdensity in contrast to regions further away from it at LMS < 230°. Bottom panel: the density profile as a function of BMS from point B0 to B1. At these large distances, and due to the steepness of the profile, the average density is dominated by the regions with the overdensities as illustrated in both panels. The lighter purple points correspond to the density at smaller patches (healpy using nside = 4) across the LMC’s orbit.

We note in passing that the Pisces overdensity region may also overlap with the Magellanic Stellar Stream (MSS, Chandra et al. 2023b) and may also contain GSE debris. As discussed previously, VOD and HAC are associated with the GSE merger, and the former shows an overall overdensity at larger radii (see Figure 4 and Vivas et al. 2001). As HAC-S and Pisces overlap in angular position on the Galactocentric coordinate system (see e.g. Figure 6 bottom right panel), it is intuitive to expect some contribution to the GSE debris on this region, potentially biasing the effect of the wake signal. Kinematics together with chemical abundance information could help disentangle the MSS from the wake of the LMC and the contribution, if any, of the GSE debris to Pisces overdensity.

Simulations of the MW-LMC interaction also predict the barycenter shift of the MW-LMC system (e.g. Gómez et al. 2015; Garavito-Camargo et al. 2019), namely the gravitational collective response or “global wake”. This effect has been observed kinematically (Erkal et al. 2021; Petersen & Peñarrubia 2021; Yaaqib et al. 2024; Chandra et al. 2024) and the density dipole associated with the collective response (Weinberg 1998; Laporte et al. 2018b; Petersen & Peñarrubia 2020) should also, in principle, be observed in star counts (Garavito-Camargo et al. 2019). Recently, Conroy et al. (2021) have reported a possible detection of this signal using K giant stars. Given that BHBs come in greater numbers, we attempt to measure the signal from the dipole in star counts.

We show in Figure 8 the approximate region in the Galactocentric sky of the region defined by Conroy et al. (2021), magenta dashed-line, which we refer to as the Northern OverDensity. The BHB density contrast map shows no apparent significant overdensity in this region, but there is an overdensity just to the north of the Northern OverDensity. Part of this region is located in the region associated with OVO, indicated by the cyan circle, and relatively close to the Sgr Stream, the latter being mostly removed. Nonetheless, this may be the remains of Sgr at very large distances (e.g. Sesar et al. 2017a).

Figure 10 shows the total density contrast measured in the Pisces overdensity region and the Northern OverDensity which are associated with the transient and global wake respectively. The purple contours show the interquartile of the density contrast measured obtained from the posteriors of the density estimation. While Pisces overdensity shows a clear signature where the density contrast is > 0.6, the signal measured in the Northern OverDensity has a large error and is consistent with being negligible. While the dipole signal in stellar density is observed in N-body simulations (e.g. Garavito-Camargo et al. 2019), the current BHB star data falls short of making a statistically meaningful measurement of this signal. This may need to await tracers that come in larger numbers in order to detect this signal, such as Main Sequence Turn Off stars.

thumbnail Fig. 10.

BHB density contrast at 50 < rq/kpc < 120 in the Pisces overdensity and the Northen overdensity regions, which are associated with the local wake and the collective response to the LMC orbit, respectively. The purple contours to the shaded region correspond to the 16th, 50th and 84th confidence intervals calculated from the BHB fraction posteriors (see text for details). This result can not confirm the collective response of the MW stellar halo to the LMC passage. For comparison, we show previous measurements using K giant stars (Conroy et al. 2021) and predictions from MW-LMC interaction models (Garavito-Camargo et al. 2019).

We also show the results predicted from Garavito-Camargo et al. (2019) simulations which have significantly lower signatures in the density map. This can be caused by the fact that their stellar halo started as a smooth component and in the MW the stellar halo is heterogeneous and thus the LMC can enhance some previously existing overdensity(ies). Indeed as shown in this study and others (e.g. Bell et al. 2008), the stellar halo beyond 50 kpc is not smooth but inhabited by substructure which due to the long orbital timescales have not had time to fully phase mix over a Hubble time (e.g. Johnston et al. 2008).

6. Conclusion

We used LSDR9 photometric survey to study the density profile of the MW stellar halo. We developed a novel approach to modelling the g − r and r − z colour distribution allowing us to compute the fraction of BHB stars at any given direction. With this approach, we can robustly estimate the errors of the densities. We summarise our main results when probing the stellar halo within 5 < rq/kpc < 120:

  • The global stellar halo shows a broken power law density profile with inner and outer slopes, and break radius, α in = 2 . 45 0.14 + 0.17 $ \alpha_{\mathrm{in}}=-2.45^{+0.17}_{-0.14} $, α in = 4 . 55 0.1 + 0.11 $ \alpha_{\mathrm{in}}=-4.55^{+0.11}_{-0.1} $, r br = 19 . 15 1.8 + 1.7 $ r_{\mathrm{br}}=19.15^{+1.7}_{-1.8} $ kpc, respectively;

  • HAC-S, HAC-N, and VOD Galactic structures show a significant break in stellar density profiles at rq ≈ 20 kpc, with compatible density inner slopes (Figure 4). The presence of these three overdensities regions remains to be seen in tailored GSE-like merger models;

  • The Pisces overdensity region is the only region where we see a density profile consistent with a single power law, in contrast to the previously mentioned overdensity regions;

  • Even after excluding regions containing major overdensities we observe a density profile with a break which seems to indicate that the GSE debris is present in most directions in the Galaxy or alternatively another structure is present in the Galaxy that also has a broken radial density profile with similar break radius;

  • We find strong evidence for an anisotropic stellar halo, indicated by significant differences in the steepness of the density profile when looking at different Galactocentric patches (Figure 7) and by ∼20% of the sky containing ∼55% of BHB stars at r > 50 kpc (Section 5.2);

  • The Pisces overdensity spans 60° along and 30° across the LMC’s orbital path in the sky (Figure 9);

  • If thought to be associated with the transient wake response of the MW to the LMC, then the Pisces overdensity extension would correspond to a wake width of ∼32 kpc at ∼70 kpc (Section 5.2). Such a measurement may help set constraints on the mass and models of the MW-LMC system and also on the nature of dark matter;

  • While we confirm the overdensity in the Pisces region, associated with the local wake induced by the LMC passage, our results do not show a statistically significant signature of the collective response in density (Figure 10).

The methodology developed in this work is a prelude to what will be possible with forthcoming photometric surveys such as the Chinese Space Station Telescope wide imaging survey (CSST, Zhan 2011; Cao et al. 2018), the Large Synoptic Survey Telescope (LSST, LSST Dark Energy Science Collaboration 2012), and Euclid space-based survey (Laureijs et al. 2011). CSST will cover 17 500 deg2 (≈40% of the sky), observe in the NUV, ugriz photometric bands, reaching a limiting magnitude of r ≈ 26.0 and r ≈ 27.2 for the wide and deep imaging survey (Qu et al. 2023), respectively. The LSST main survey will cover 30 000 deg2 and reach up to r ≈ 27.5 at the end of its mission. The depth to be reached by these surveys, with lower magnitude errors compared to LSDR9, will allow us to probe the stellar halo fraction at a higher spatial angular resolution and larger distances. This will refine the measurement of the gravitational wake width in the MW halo due to the LMC, as shown in Figure 9, particularly when extending this work’s methodology to Main Sequence Turn Off stars. These stars are not only more numerous than BHB stars (e.g. Bell et al. 2010), but they will reduce the Poisson noise and allow for finer measurements of density variations across/along the wake and across overdensities (Sanderson et al. 2019).

Finally, we present in Appendix C a catalogue of BHB star candidates with their respective probabilities calculated from Equation (3).

Data availability

Full Table C.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/690/A166 and at https://zenodo.org/records/12155488


1

LSDR9 uses Schlegel et al. (1998). Galactic extinction maps and provides the coefficients at https://www.legacysurvey.org/dr9/catalogs/#galactic-extinction-coefficients

2

All photometric magnitudes are extinction corrected thus we opted not to write them with the subscript “0”.

4

For a Gaussian distribution this relation corresponds to 1σ.

5

The healpy package orders the healpixels either in a nested or ring scheme.

6

The contribution of the substructure to a pixel is computed as a fraction of the pixel area that intersects with the substructure area given in Table 1.

7

Nested healpix number 33, for nside = 2.

8

We made use of the LMC-MW interaction example of the AGAMA python library (Vasiliev 2019).

9

In this reference system, LMC is located at LMS ≈ 0°.

Acknowledgments

The authors wish to thank the anonymous referee for their comments which helped improve the quality of this work. J.A. and C.L. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 852839). SK acknowledges support from the Science & Technology Facilities Council (STFC) grant ST/Y001001/1. J.A. thanks Guilherme Limberg, Hélio Perottoni, Friedrich Anders, and Matt Orkney for conversations that contributed to improve the quality of this work. This paper made use of the Whole Sky Database (wsdb) created and maintained by Sergey Koposov at the Institute of Astronomy, Cambridge with financial support from the Science & Technology Facilities Council (STFC) and the European Research Council (ERC). Some of the results in this paper have been derived using the healpy and HEALPix packages. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 12120101003, # 11433005). The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

References

  1. Amarante, J. A. S., Smith, M. C., & Boeche, C. 2020, MNRAS, 492, 3816 [NASA ADS] [CrossRef] [Google Scholar]
  2. Amarante, J. A. S., Debattista, V. P., Beraldo e Silva, L., Laporte, C. F. P., & Deg, N. 2022, ApJ, 937, 12 [NASA ADS] [CrossRef] [Google Scholar]
  3. Anders, F., Cantat-Gaudin, T., Quadrino-Lodoso, I., et al. 2021, A&A, 645, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Barbosa, F. O., Santucci, R. M., Rossi, S., et al. 2022, ApJ, 940, 30 [NASA ADS] [CrossRef] [Google Scholar]
  5. Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 680, 295 [Google Scholar]
  7. Bell, E. F., Xue, X. X., Rix, H.-W., Ruhland, C., & Hogg, D. W. 2010, AJ, 140, 1850 [Google Scholar]
  8. Belokurov, V., & Koposov, S. E. 2016, MNRAS, 456, 602 [NASA ADS] [CrossRef] [Google Scholar]
  9. Belokurov, V., Evans, N. W., Bell, E. F., et al. 2007a, ApJ, 657, L89 [CrossRef] [Google Scholar]
  10. Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007b, ApJ, 654, 897 [Google Scholar]
  11. Belokurov, V., Koposov, S. E., Evans, N. W., et al. 2014, MNRAS, 437, 116 [NASA ADS] [CrossRef] [Google Scholar]
  12. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611 [Google Scholar]
  13. Belokurov, V., Deason, A. J., Erkal, D., et al. 2019, MNRAS, 488, L47 [Google Scholar]
  14. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [Google Scholar]
  15. Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668, 949 [Google Scholar]
  16. Bird, S. A., Xue, X.-X., Liu, C., et al. 2021, ApJ, 919, 66 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bonaca, A., Jurić, M., Ivezić, Ž., et al. 2012, AJ, 143, 105 [NASA ADS] [CrossRef] [Google Scholar]
  18. Borbolato, L., Perottoni, H. D., Rossi, S., et al. 2024, ApJ, 960, 52 [NASA ADS] [CrossRef] [Google Scholar]
  19. Brook, C. B., Kawata, D., Gibson, B. K., & Flynn, C. 2003, ApJ, 585, L125 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931 [Google Scholar]
  21. Cao, Y., Gong, Y., Meng, X.-M., et al. 2018, MNRAS, 480, 2178 [Google Scholar]
  22. Carr, C., Johnston, K. V., Laporte, C. F. P., & Ness, M. K. 2022, MNRAS, 516, 5067 [NASA ADS] [CrossRef] [Google Scholar]
  23. Carrillo, A., Deason, A. J., Fattahi, A., Callingham, T. M., & Grand, R. J. J. 2024, MNRAS, 527, 2165 [Google Scholar]
  24. Cerny, W., Martínez-Vázquez, C. E., Drlica-Wagner, A., et al. 2023, ApJ, 953, 1 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chandra, V., Naidu, R. P., Conroy, C., et al. 2023a, ApJ, 951, 26 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chandra, V., Naidu, R. P., Conroy, C., et al. 2023b, ApJ, 956, 110 [CrossRef] [Google Scholar]
  27. Chandra, V., Naidu, R. P., Conroy, C., et al. 2024, ArXiv e-prints [arXiv:2406.01676] [Google Scholar]
  28. Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
  29. Chiba, M., & Beers, T. C. 2000, AJ, 119, 2843 [NASA ADS] [CrossRef] [Google Scholar]
  30. Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Nature, 592, 534 [NASA ADS] [CrossRef] [Google Scholar]
  31. Cunningham, E. C., Garavito-Camargo, N., Deason, A. J., et al. 2020, ApJ, 898, 4 [CrossRef] [Google Scholar]
  32. Das, P., Hawkins, K., & Jofré, P. 2020, MNRAS, 493, 5195 [NASA ADS] [CrossRef] [Google Scholar]
  33. da Silva, A. R., & Smiljanic, R. 2023, A&A, 677, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903 [NASA ADS] [CrossRef] [Google Scholar]
  35. Deason, A. J., Belokurov, V., Evans, N. W., & Johnston, K. V. 2013, ApJ, 763, 113 [Google Scholar]
  36. Deason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30 [NASA ADS] [CrossRef] [Google Scholar]
  37. Deason, A. J., Belokurov, V., & Koposov, S. E. 2018, ApJ, 852, 118 [NASA ADS] [CrossRef] [Google Scholar]
  38. Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426 [NASA ADS] [CrossRef] [Google Scholar]
  39. Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168 [Google Scholar]
  40. Donlon, T., & Newberg, H. J. 2023, ApJ, 944, 169 [NASA ADS] [CrossRef] [Google Scholar]
  41. Donlon, T., Newberg, H. J., Sanderson, R., & Widrow, L. M. 2020, ApJ, 902, 119 [NASA ADS] [CrossRef] [Google Scholar]
  42. Donlon, T., Newberg, H. J., Sanderson, R., et al. 2024, MNRAS, 531, 1422 [NASA ADS] [CrossRef] [Google Scholar]
  43. Dotter, A., Sarajedini, A., Anderson, J., et al. 2010, ApJ, 708, 698 [Google Scholar]
  44. Erkal, D., Deason, A. J., Belokurov, V., et al. 2021, MNRAS, 506, 2677 [NASA ADS] [CrossRef] [Google Scholar]
  45. Feuillet, D. K., Feltzing, S., Sahlholdt, C. L., & Casagrande, L. 2020, MNRAS, 497, 109 [Google Scholar]
  46. Forbes, D. A. 2020, MNRAS, 493, 847 [Google Scholar]
  47. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  48. Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631 [NASA ADS] [CrossRef] [Google Scholar]
  49. Fukushima, T., Chiba, M., Homma, D., et al. 2018, PASJ, 70, 69 [NASA ADS] [CrossRef] [Google Scholar]
  50. Fukushima, T., Chiba, M., Tanaka, M., et al. 2019, PASJ, 71, 72 [NASA ADS] [CrossRef] [Google Scholar]
  51. Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, Nat. Astron., 3, 932 [NASA ADS] [CrossRef] [Google Scholar]
  52. Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51 [NASA ADS] [CrossRef] [Google Scholar]
  53. Garavito-Camargo, N., Patel, E., Besla, G., et al. 2021, ApJ, 923, 140 [NASA ADS] [CrossRef] [Google Scholar]
  54. Gentile Fusillo, N. P., Tremblay, P. E., Cukanovaite, E., et al. 2021, MNRAS, 508, 3877 [NASA ADS] [CrossRef] [Google Scholar]
  55. Gilmore, G., Wyse, R. F. G., & Norris, J. E. 2002, ApJ, 574, L39 [NASA ADS] [CrossRef] [Google Scholar]
  56. Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [Google Scholar]
  57. Gómez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802, 128 [CrossRef] [Google Scholar]
  58. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
  59. Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179 [NASA ADS] [Google Scholar]
  60. Han, J. J., Conroy, C., Johnson, B. D., et al. 2022, AJ, 164, 249 [NASA ADS] [CrossRef] [Google Scholar]
  61. Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113 [Google Scholar]
  62. Helmi, A., White, S. D. M., de Zeeuw, P. T., & Zhao, H. 1999, Nature, 402, 53 [Google Scholar]
  63. Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85 [Google Scholar]
  64. Hernitschek, N., Sesar, B., Rix, H.-W., et al. 2017, ApJ, 850, 96 [NASA ADS] [CrossRef] [Google Scholar]
  65. Hernitschek, N., Cohen, J. G., Rix, H.-W., et al. 2018, ApJ, 859, 31 [NASA ADS] [CrossRef] [Google Scholar]
  66. Huang, K.-W., & Koposov, S. E. 2022, MNRAS, 510, 3575 [NASA ADS] [CrossRef] [Google Scholar]
  67. Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194 [Google Scholar]
  68. Iorio, G., Belokurov, V., Erkal, D., et al. 2018, MNRAS, 474, 2142 [NASA ADS] [CrossRef] [Google Scholar]
  69. Johnston, K. V., Spergel, D. N., & Hernquist, L. 1995, ApJ, 451, 598 [NASA ADS] [CrossRef] [Google Scholar]
  70. Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936 [Google Scholar]
  71. Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [Google Scholar]
  72. Kollmeier, J. A., Gould, A., Shectman, S., et al. 2009, ApJ, 705, L158 [NASA ADS] [CrossRef] [Google Scholar]
  73. Koposov, S. E., Belokurov, V., Torrealba, G., & Evans, N. W. 2015, ApJ, 805, 130 [NASA ADS] [CrossRef] [Google Scholar]
  74. Koppelman, H. H., Helmi, A., Massari, D., Price-Whelan, A. M., & Starkenburg, T. K. 2019, A&A, 631, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Lane, J. M. M., Bovy, J., & Mackereth, J. T. 2023, MNRAS, 526, 1209 [NASA ADS] [CrossRef] [Google Scholar]
  76. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018a, MNRAS, 481, 286 [Google Scholar]
  77. Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V., & Garavito-Camargo, N. 2018b, MNRAS, 473, 1218 [NASA ADS] [CrossRef] [Google Scholar]
  78. Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134 [Google Scholar]
  79. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193] [Google Scholar]
  80. Li, T. S., Koposov, S. E., Zucker, D. B., et al. 2019, MNRAS, 490, 3508 [NASA ADS] [CrossRef] [Google Scholar]
  81. Limberg, G., Rossi, S., Beers, T. C., et al. 2021, ApJ, 907, 10 [Google Scholar]
  82. Limberg, G., Souza, S. O., Pérez-Villegas, A., et al. 2022, ApJ, 935, 109 [NASA ADS] [CrossRef] [Google Scholar]
  83. Liu, G., Huang, Y., Bird, S. A., et al. 2022, MNRAS, 517, 2787 [NASA ADS] [CrossRef] [Google Scholar]
  84. LSST Dark Energy Science Collaboration 2012, ArXiv e-prints [arXiv:1211.0310] [Google Scholar]
  85. Malhan, K., Ibata, R. A., Sharma, S., et al. 2022, ApJ, 926, 107 [NASA ADS] [CrossRef] [Google Scholar]
  86. Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. McConnachie, A. W. 2012, AJ, 144, 4 [Google Scholar]
  88. McMillan, P. J. 2017, MNRAS, 465, 76 [Google Scholar]
  89. Medina, G. E., Muñoz, R. R., Vivas, A. K., et al. 2018, ApJ, 855, 43 [NASA ADS] [CrossRef] [Google Scholar]
  90. Medina, G. E., Muñoz, R. R., Carlin, J. L., et al. 2024, MNRAS, 531, 4762 [NASA ADS] [CrossRef] [Google Scholar]
  91. Meza, A., Navarro, J. F., Abadi, M. G., & Steinmetz, M. 2005, MNRAS, 359, 93 [CrossRef] [Google Scholar]
  92. Montalbán, J., Mackereth, J. T., Miglio, A., et al. 2021, Nat. Astron., 5, 640 [Google Scholar]
  93. Moustakas, J., Lang, D., Dey, A., et al. 2023, ApJS, 269, 3 [NASA ADS] [CrossRef] [Google Scholar]
  94. Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, ApJ, 863, L28 [NASA ADS] [CrossRef] [Google Scholar]
  95. Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48 [Google Scholar]
  96. Naidu, R. P., Conroy, C., Bonaca, A., et al. 2021, ApJ, 923, 92 [NASA ADS] [CrossRef] [Google Scholar]
  97. Naidu, R. P., Ji, A. P., Conroy, C., et al. 2022, ApJ, 926, L36 [CrossRef] [Google Scholar]
  98. Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245 [Google Scholar]
  99. Nidever, D. L., Majewski, S. R., & Burton, W. B. 2008, Astrophys. Space Sci. Proc., 5, 243 [NASA ADS] [CrossRef] [Google Scholar]
  100. Nie, J. D., Smith, M. C., Belokurov, V., et al. 2015, ApJ, 810, 153 [NASA ADS] [CrossRef] [Google Scholar]
  101. Orkney, M. D. A., Laporte, C. F. P., Grand, R. J. J., et al. 2023, MNRAS, 525, 683 [NASA ADS] [CrossRef] [Google Scholar]
  102. Pace, A. B., Erkal, D., & Li, T. S. 2022, ApJ, 940, 136 [NASA ADS] [CrossRef] [Google Scholar]
  103. Perottoni, H. D., Limberg, G., Amarante, J. A. S., et al. 2022, ApJ, 936, L2 [NASA ADS] [CrossRef] [Google Scholar]
  104. Petersen, M. S., & Peñarrubia, J. 2020, MNRAS, 494, L11 [NASA ADS] [CrossRef] [Google Scholar]
  105. Petersen, M. S., & Peñarrubia, J. 2021, Nat. Astron., 5, 251 [NASA ADS] [CrossRef] [Google Scholar]
  106. Petersen, M. S., Peñarrubia, J., & Jones, E. 2022, MNRAS, 514, 1266 [CrossRef] [Google Scholar]
  107. Pieres, A., Girardi, L., Balbinot, E., et al. 2020, MNRAS, 497, 1547 [NASA ADS] [CrossRef] [Google Scholar]
  108. Pila-Díez, B., de Jong, J. T. A., Kuijken, K., van der Burg, R. F. J., & Hoekstra, H. 2015, A&A, 579, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Purcell, C. W., Bullock, J. S., Tollerud, E. J., Rocha, M., & Chakrabarti, S. 2011, Nature, 477, 301 [Google Scholar]
  110. Qu, H., Yuan, Z., Doliva-Dolinsky, A., et al. 2023, MNRAS, 523, 876 [NASA ADS] [CrossRef] [Google Scholar]
  111. Ross, N. P., Myers, A. D., Sheldon, E. S., et al. 2012, ApJS, 199, 3 [NASA ADS] [CrossRef] [Google Scholar]
  112. Rozier, S., Famaey, B., Siebert, A., et al. 2022, ApJ, 933, 113 [NASA ADS] [CrossRef] [Google Scholar]
  113. Sanderson, R., Carlin, J., Cunningham, E., et al. 2019, BAAS, 51, 347 [NASA ADS] [Google Scholar]
  114. Santucci, R. M., Beers, T. C., Placco, V. M., et al. 2015, ApJ, 813, L16 [NASA ADS] [CrossRef] [Google Scholar]
  115. Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [Google Scholar]
  116. Searle, L., & Zinn, R. 1978, ApJ, 225, 357 [Google Scholar]
  117. Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, AJ, 134, 2236 [NASA ADS] [CrossRef] [Google Scholar]
  118. Sesar, B., Ivezić, Ž., Grammer, S. H., et al. 2010, ApJ, 708, 717 [Google Scholar]
  119. Sesar, B., Hernitschek, N., Mitrović, S., et al. 2017a, AJ, 153, 204 [NASA ADS] [CrossRef] [Google Scholar]
  120. Sesar, B., Hernitschek, N., Dierickx, M. I. P., Fardal, M. A., & Rix, H.-W. 2017b, ApJ, 844, L4 [Google Scholar]
  121. Simion, I. T., Belokurov, V., Irwin, M., & Koposov, S. E. 2014, MNRAS, 440, 161 [NASA ADS] [CrossRef] [Google Scholar]
  122. Simion, I. T., Belokurov, V., & Koposov, S. E. 2019, MNRAS, 482, 921 [NASA ADS] [CrossRef] [Google Scholar]
  123. Sirko, E., Goodman, J., Knapp, G. R., et al. 2004, AJ, 127, 899 [NASA ADS] [CrossRef] [Google Scholar]
  124. Springel, V., Frenk, C. S., & White, S. D. M. 2006, Nature, 440, 1137 [NASA ADS] [CrossRef] [Google Scholar]
  125. Starkenburg, E., Youakim, K., Martin, N., et al. 2019, MNRAS, 490, 5757 [Google Scholar]
  126. Stringer, K. M., Drlica-Wagner, A., Macri, L., et al. 2021, ApJ, 911, 109 [Google Scholar]
  127. Suzuki, Y., Chiba, M., Komiyama, Y., et al. 2024, PASJ, 76, 205 [NASA ADS] [CrossRef] [Google Scholar]
  128. Thomas, G. F., McConnachie, A. W., Ibata, R. A., et al. 2018, MNRAS, 481, 5223 [NASA ADS] [CrossRef] [Google Scholar]
  129. Vasiliev, E. 2019, MNRAS, 482, 1525 [Google Scholar]
  130. Vasiliev, E. 2023, Galaxies, 11, 59 [NASA ADS] [CrossRef] [Google Scholar]
  131. Vickers, J. J., Grebel, E. K., & Huxor, A. P. 2012, AJ, 143, 86 [NASA ADS] [CrossRef] [Google Scholar]
  132. Vivas, A. K., Zinn, R., Andrews, P., et al. 2001, ApJ, 554, L33 [NASA ADS] [CrossRef] [Google Scholar]
  133. Watkins, L. L., Evans, N. W., Belokurov, V., et al. 2009, MNRAS, 398, 1757 [NASA ADS] [CrossRef] [Google Scholar]
  134. Weinberg, M. D. 1989, MNRAS, 239, 549 [NASA ADS] [CrossRef] [Google Scholar]
  135. Weinberg, M. D. 1998, MNRAS, 299, 499 [Google Scholar]
  136. Weinberg, M. D., & Blitz, L. 2006, ApJ, 641, L33 [NASA ADS] [CrossRef] [Google Scholar]
  137. White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52 [Google Scholar]
  138. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
  139. Wu, W., Zhao, G., Xue, X.-X., Pei, W., & Yang, C. 2022a, AJ, 164, 41 [CrossRef] [Google Scholar]
  140. Wu, W., Zhao, G., Xue, X.-X., Bird, S. A., & Yang, C. 2022b, ApJ, 924, 23 [NASA ADS] [CrossRef] [Google Scholar]
  141. Xiang, M., & Rix, H.-W. 2022, Nature, 603, 599 [NASA ADS] [CrossRef] [Google Scholar]
  142. Xue, X.-X., Rix, H.-W., Yanny, B., et al. 2011, ApJ, 738, 79 [NASA ADS] [CrossRef] [Google Scholar]
  143. Xue, X.-X., Rix, H.-W., Ma, Z., et al. 2015, ApJ, 809, 144 [Google Scholar]
  144. Yaaqib, R., Petersen, M. S., & Peñarrubia, J. 2024, MNRAS, 531, 3524 [NASA ADS] [CrossRef] [Google Scholar]
  145. Yang, C., Zhu, L., Tahmasebzadeh, B., Xue, X.-X., & Liu, C. 2022, AJ, 164, 241 [NASA ADS] [CrossRef] [Google Scholar]
  146. Yanny, B., Newberg, H. J., Kent, S., et al. 2000, ApJ, 540, 825 [NASA ADS] [CrossRef] [Google Scholar]
  147. Yu, F., Li, T. S., Speagle, J. S., et al. 2024, ArXiv e-prints [arXiv:2402.00104] [Google Scholar]
  148. Yuan, Z., Myeong, G. C., Beers, T. C., et al. 2020, ApJ, 891, 39 [NASA ADS] [CrossRef] [Google Scholar]
  149. Zhan, H. 2011, Sci. Sin. Phys. Mech. Astron., 41, 1441 [NASA ADS] [CrossRef] [Google Scholar]
  150. Zhou, R., Dey, A., Lang, D., et al. 2023, Res. Notes Am. Astron. Soc., 7, 105 [Google Scholar]
  151. Zinn, R., Horowitz, B., Vivas, A. K., et al. 2014, ApJ, 781, 22 [NASA ADS] [CrossRef] [Google Scholar]
  152. Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298 [Google Scholar]

Appendix A: White dwarf contamination

In this section, we investigate the possible white dwarf contamination in identifying BHB stars through our methodology described in Section 2.2. We use the Gaia EDR3 white dwarf catalogue (Gentile Fusillo et al. 2021) to verify whether this stellar population could contaminate the BHB fraction estimation. As we are not interested in the purity of a white dwarf sample, but rather its contamination, we used all their white dwarf candidates to crossmatch with our LSDR9 sample.

The left column in Figure A.1 shows the grz − (g − r) plane, at different g intervals, the white dwarf candidates. The black rectangle highlights the region of interest in our study and in the right column, we show the grz histogram, in light blue, for the objects inside it. We also show, in blue, the histogram for all star objects excluding the white dwarfs, that lie inside the selected region. The solid blue line is the histogram of all objects in the selected grz − (g − r) region. It is clear that for g > 17 there is a low contamination fraction of white dwarf candidates in the expected BHB region.

thumbnail Fig. A.1.

Left column: The grz − (g − r) plane for the white dwarf stars in from the Gaia EDR3 catalogue (Gentile Fusillo et al. 2021). The black rectangle shows the region of interest in our study. Right column: The grz histogram for objects inside the black rectangle. The white dwarf and non-white dwarf objects histograms are shown as the light and dark blue histograms, respectively. The solid blue line corresponds to their sum. The solid red line marks grz = 0, the centre of the BHB distribution.

We verify the effect of not excluding white dwarf candidates on the BHB fraction estimation by fitting the grz distribution for several intervals in 13 < g < 21. We use the same methodology described in Section 3.1 and also include strong prior in the BHB Gaussian parameters. Figure A.2 shows the fractions for the BHB, BS, and background components in blue, orange, and magenta, respectively. The dark/light markers exclude/include white dwarf candidates. We can see that the BHB fraction is independent of the removal or not of white dwarf candidates, at least in the colour-colour regime of our study. The most noticeable effect is on the BS and background contamination. We, therefore, verify that white dwarfs are not significant sources of contamination in our study.

thumbnail Fig. A.2.

Fractions of each component from the MCMC modelling, described in Section 3.1, for objects with −0.3 < (g − r) <  − 0.05 and −0.14 < grz < 0.07 in the LSDR9-Gaia EDR3 crossmatch. The BHB Gaussian component has a hard prior for the mean and standard deviation. The dark/light markers are the fractions for a given g interval when white dwarf candidates are excluded/included. The white dwarf contamination in the BHB fractions is negligible throughout g.

Appendix B: Radial density profiles

In Section 4.2, we divided the Galactocentric sky into 48 patches with the Python package healpy (Górski et al. 2005); healpy using nside = 2. Each region has an equal projected area of 859.44 deg2. We constructed the density profiles, see Section 3 for details, for 10 bins equally space in logarithm in the range 5 < rq/kpc < 120. Figure B.1 shows the profiles of all the patches which have four or more bins in rq.

thumbnail Fig. B.1.

Radial density profiles for several patches in the Galactocentric sky as in Figure 6. We show all the regions with more than four radial bins. See Section 4.2 for details.

Figure B.2 shows the inner and outer slopes, from the broken-power law fitting (see Section 3.4) for each patch of the Galactocentric sky with four or more radial bins. The coordinates of the centre of each patch are given on the x-axis in degrees. The panel highlights the highly anisotropic stellar halo, as observed with BHB stars and discussed throughout in Section 5. We also show the slopes obtained by several studies as compiled by Yu et al. (2024).

thumbnail Fig. B.2.

Inner (outer) slopes of the broken power law model fit of the density profiles shown in Figure B.1 shown as purple (green) circles. The x-axis shows the coordinates in degrees of the centre of the Galactocentric patch using healpy Python package with HEALPIX pixelisation nside = 2. The shaded region corresponds to the slopes measured for the average radial density profile as in Section 4.1. We also included the list compiled by Yu et al. (2024) of several works studying the halo radial density profile: Watkins et al. (2009), Deason et al. (2011), Xue et al. (2015), Pila-Díez et al. (2015), Deason et al. (2018), Fukushima et al. (2018), Hernitschek et al. (2018), Iorio et al. (2018), Medina et al. (2018), Thomas et al. (2018), Starkenburg et al. (2019), Pieres et al. (2020), Stringer et al. (2021), Han et al. (2022), Yu et al. (2024), Medina et al. (2024). These are shown as open diamonds. Han et al. (2022) favours a three-power-law model and thus the extra slope is shown as the black diamond. The dashed vertical line separates data from our work from literature values.

Appendix C: LSDR9 BHB star candidates catalogue

We use objects from Legacy Survey Data Release 9 (Dey et al. 2019) with the colour selection of Equation 2. We then construct a catalogue with 95446 objects containing LSDR9 release, brickid and objid. These three parameters provide a unique identifier hash10. The catalogue also includes R.A., Dec., the extinction corrected g0, r0, z0, grz (Equation 1), the absolute Mg magnitude (Equation 6), and the probability of a stellar object being a BHB star, PBHB, calculated with Equation 3. Table C.1 shows the header and first ten rows of the catalogue available online11.

Table C.1.

First 10 rows of the catalogue constructed with LSDR9 objects. See Appendix C for details.

All Tables

Table 1.

Galactic coordinates and heliocentric distance used in this work for Hecules-Aquila Cloud North and South (HAC-N and HAC-S, Belokurov et al. 2007a; Watkins et al. 2009; Sesar et al. 2010; Simion et al. 2014), Virgo OverDensity (VOD, Jurić et al. 2008; Bonaca et al. 2012), and Pisces (Nie et al. 2015; Belokurov et al. 2019).

Table C.1.

First 10 rows of the catalogue constructed with LSDR9 objects. See Appendix C for details.

All Figures

thumbnail Fig. 1.

(g − r)×(r − z) diagram for the first selected LSDR9 objects. The blue dotted lines are the boundaries for selecting BHB star candidates following Li et al. (2019).

In the text
thumbnail Fig. 2.

grz-colour distribution. The left panels show grz colour, see Equation (1), as a function (g − r) at different g ranges. The right panels show the P(grz), in grey, of the stars within the selection box on the left panel. The corresponding fit to the data following Equation (3) as well as the BHB, BS and background component are represented in green, blue, orange and magenta, respectively. These two examples are part of the exercise illustrated in the left column of Figure 3 which assumes flat priors for the free parameters. Throughout this paper, we do use Gaussian priors for the BHB stellar component when calculating the BHB stellar densities, see Section 3.1 for details.

In the text
thumbnail Fig. 3.

Left column: the variation of the free parameters while fitting Equation (3) as a function of g assuming flat priors. While the horizontal lines correspond to the g interval on which the fit was performed, the vertical lines correspond to the 16th–84th interval of the posterior distribution. The BHB stellar parameters are relatively constant within 15 < g < 19 which motivated the use of Gaussian priors throughout this work, see Section 3.1 for discussion. Right column: now with Gaussian priors on the BHB component and within 14.5 < r < 21.5. The decrease in the BHB fraction as a function of magnitude is expected due to the intrinsic density profile of the MW halo.

In the text
thumbnail Fig. 4.

Top left panel: the Galactic Mollweide projection of LSDR9 stellar objects selected using Equation (2). The highlighted regions represent the location of known substructures: Sagittarius stream (pink), HAC-S (yellow), HAC-N (red), Pisces (green), VOD (blue), and OVO (cyan). Top right panel: the Galactocentric BHB stellar density profile of the whole LSDR9 footprint and without HAC-S, HAC-N, and VOD regions as the solid and dotted markers, respectively. Bottom panels: the Galactocentric radial density profiles for the Galactic regions associated with HAC-N, HAC-S, VOD, and Pisces. The crosses correspond to twenty and ten radial bins equally spaced in logarithm within 5 < rq/kpc < 120 for LSDR9 footprint and overdensities, respectively. We also show the resulting fit for 100 samples drawn from the posterior distribution as the solid lines.

In the text
thumbnail Fig. 5.

Galactocentric radial density profiles for the four quadrants of the Galactocentric projected sky. The crosses correspond to ten radial bins equally spaced in logarithm within 5 < rq/kpc < 120. We also show the resulting fit for 100 samples drawn from the posterior distribution as the brown solid lines. The histogram corresponds to the posterior distribution of the break radius, rbr, and the dashed vertical lines are the corresponding 16th and 84th percentiles. The density profile of the full LSDR9 footprint, HAC-N, HAC-S, VOD, and Pisces are shown as the black, red, yellow, blue, and green lines, respectively.

In the text
thumbnail Fig. 6.

Radial density profile for several patches in the Galactocentric sky covering an area of 859.44 deg2. The purple crosses correspond to the radial densities of the highlighted region shown in the Mollweide map of each panel. The purple solid lines correspond to the double power law fit using 100 samples drawn from the parameters’ posterior distribution, and the dashed vertical lines correspond to the 16th and 84th percentiles of the rbr posterior distribution. The black solid line is the best fit performed on the whole LSDR9 footprint. Each panel also shows the contribution of known substructures in a given region: HAC-N (red), VOD (blue), HAC-S (yellow), Pisces (green), and Sgr stream (pink). The darker the highlighted region is, the stronger the contribution of the substructure on the highlighted Galactocentric region.

In the text
thumbnail Fig. 7.

Mollweide projection of the Galactocentric coordinates (Φgc, Θgc). The projection is divided into a grid of 48 pixels (nside = 2) with equal area. The blue, red, yellow, and orange regions delineate the regions associated with VOD, HAC-N, HAC-S, and Pisces at their typical distances. Top panel: the inner slope, αin, measured at rq = 15 kpc, of the two power law model. Middle panel: the outer slope, αout, measured at rq = 50 kpc, of the two power law model. Bottom panel: the difference between αin and αout. These panels highlight the anisotropy of the MW stellar halo profile, where several patches of the sky have distinct slopes compared to the average halo. We only show patches with more than four density points in the radial binning. The crosses and “x” indicate the regions with |αout|−|αin|< 1.5, and |αout|−|αin|> 1.5 with error in the break radius greater than 25%, respectively.

In the text
thumbnail Fig. 8.

Density contrast of MW BHB stars in the stellar halo at 50 < rq/kpc < 120 seen in Galactocentric coordinates (ϕgc, θgc). The Mollweide projection is divided into a grid of 192 pixels (nside = 4) with equal area. We delineate the approximate Galactocentric regions associated with OVO, Pisces, and Northern OverDensity as the cyan, green, and magenta lines, respectively. The present-day position of the LMC and SMC are shown as the purple scatter points. LMC’s past trajectory (see Vasiliev 2023 for details) and Sgr stream are represented by the solid purple and dashed black lines, respectively.

In the text
thumbnail Fig. 9.

BHB star density of the sky patches crossed by the LMC’s past orbit (Vasiliev 2023) as a function of the (LMS, BMS) coordinates (Nidever et al. 2008). Each patch has an equal projected area of 859.44 deg2 (healpy using nside = 2) and the density is calculated within 50 < r/kpc < 120. The 16th–84th percentiles of the average BHB density within this distance is shown by the shaded green line. Top panel: the density profile along LMC’s orbit from point A to the pixel containing LMC. At B1, the pixel dominated by Pisces overdensity shows a clear signature of the overdensity in contrast to regions further away from it at LMS < 230°. Bottom panel: the density profile as a function of BMS from point B0 to B1. At these large distances, and due to the steepness of the profile, the average density is dominated by the regions with the overdensities as illustrated in both panels. The lighter purple points correspond to the density at smaller patches (healpy using nside = 4) across the LMC’s orbit.

In the text
thumbnail Fig. 10.

BHB density contrast at 50 < rq/kpc < 120 in the Pisces overdensity and the Northen overdensity regions, which are associated with the local wake and the collective response to the LMC orbit, respectively. The purple contours to the shaded region correspond to the 16th, 50th and 84th confidence intervals calculated from the BHB fraction posteriors (see text for details). This result can not confirm the collective response of the MW stellar halo to the LMC passage. For comparison, we show previous measurements using K giant stars (Conroy et al. 2021) and predictions from MW-LMC interaction models (Garavito-Camargo et al. 2019).

In the text
thumbnail Fig. A.1.

Left column: The grz − (g − r) plane for the white dwarf stars in from the Gaia EDR3 catalogue (Gentile Fusillo et al. 2021). The black rectangle shows the region of interest in our study. Right column: The grz histogram for objects inside the black rectangle. The white dwarf and non-white dwarf objects histograms are shown as the light and dark blue histograms, respectively. The solid blue line corresponds to their sum. The solid red line marks grz = 0, the centre of the BHB distribution.

In the text
thumbnail Fig. A.2.

Fractions of each component from the MCMC modelling, described in Section 3.1, for objects with −0.3 < (g − r) <  − 0.05 and −0.14 < grz < 0.07 in the LSDR9-Gaia EDR3 crossmatch. The BHB Gaussian component has a hard prior for the mean and standard deviation. The dark/light markers are the fractions for a given g interval when white dwarf candidates are excluded/included. The white dwarf contamination in the BHB fractions is negligible throughout g.

In the text
thumbnail Fig. B.1.

Radial density profiles for several patches in the Galactocentric sky as in Figure 6. We show all the regions with more than four radial bins. See Section 4.2 for details.

In the text
thumbnail Fig. B.2.

Inner (outer) slopes of the broken power law model fit of the density profiles shown in Figure B.1 shown as purple (green) circles. The x-axis shows the coordinates in degrees of the centre of the Galactocentric patch using healpy Python package with HEALPIX pixelisation nside = 2. The shaded region corresponds to the slopes measured for the average radial density profile as in Section 4.1. We also included the list compiled by Yu et al. (2024) of several works studying the halo radial density profile: Watkins et al. (2009), Deason et al. (2011), Xue et al. (2015), Pila-Díez et al. (2015), Deason et al. (2018), Fukushima et al. (2018), Hernitschek et al. (2018), Iorio et al. (2018), Medina et al. (2018), Thomas et al. (2018), Starkenburg et al. (2019), Pieres et al. (2020), Stringer et al. (2021), Han et al. (2022), Yu et al. (2024), Medina et al. (2024). These are shown as open diamonds. Han et al. (2022) favours a three-power-law model and thus the extra slope is shown as the black diamond. The dashed vertical line separates data from our work from literature values.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.