The stellar mass of the Gaia-Sausage/Enceladus accretion remnant

ABSTRACT

The Gaia-Sausage/Enceladus (GS/E) structure is an accretion remnant that comprises a large fraction of the Milky Way’s stellar halo. We study GS/E using high-purity samples of kinematically selected stars from APOGEE DR16 and Gaia. Employing a novel framework to account for kinematic selection biases using distribution functions, we fit density profiles to these GS/E samples and measure their masses. We find that GS/E has a shallow density profile in the inner Galaxy, with a break between 15 and 25 kpc beyond which the profile steepens. We also find that GS/E is triaxial, with axis ratios 1:0.55:0.45 (nearly prolate), and the major axis is oriented about 80° from the Sun–Galactic centre line and 16° above the plane. We measure a stellar mass for GS/E of |$1.45\, ^{+0.92}_{-0.51}\, \mathrm{(stat.)}\, ^{+0.13}_{-0.37} \mathrm{(sys.)}\ \times 10^{8}$| M_⊙. Our mass estimate is lower than others in the literature, a finding we attribute to the excellent purity of the samples we work with. We also fit a density profile to the entire Milky Way stellar halo, finding a mass in the range of 6.7–8.4 × 10⁸ M_⊙, and implying that GS/E could make up as little as 15–25 per cent of the mass of the Milky Way stellar halo. Our lower stellar mass combined with standard stellar mass-to-halo mass relations implies that GS/E constituted a minor 1:8 mass-ratio merger at the time of its accretion.

1 INTRODUCTION

In our ΛCDM Universe, in which structure is arranged hierarchically, the constituent mass of a galaxy like the Milky Way grows largely due to accretion of smaller nearby structures, which supply gas, stars, and dark matter to the gravitationally dominant host (Searle & Zinn 1978; White & Frenk 1991; Helmi et al. 1999; Bullock & Johnston 2005). One of the best records of these accretion processes lies in the distribution and arrangement of stars in the Milky Way halo, which have comparatively dissapationless kinematics and have therefore retained signatures of the events which brought them into the Galaxy (Freeman & Bland-Hawthorn 2002; Johnston et al. 2008). Mergers play an outsized role in shaping our Galaxy beyond its stellar and dark halo, and are thought to contribute to the formation or evolution of the bulge, bar, disc, as well as the generation of past and ongoing instabilities and perturbations to these stellar components (Bland-Hawthorn & Gerhard 2016; Helmi 2020). An inventory of the stellar constituents of the Milky Way halo is therefore crucial to understanding the role of mergers in the broader context of the formation and evolution of our Galaxy.

The combination of astrometric data from the Gaia mission (Gaia Collaboration et al. 2016) and large ground-based spectroscopic surveys such as SEGUE (Yanny et al. 2009), APOGEE (Majewski et al. 2017), GALAH (Martell et al. 2017), LAMOST (Cui et al. 2012), and H3 (Conroy et al. 2019), has proved a boon for studying the Milky Way stellar halo. One of the most intriguing results of this new era is the detailed characterization using data from the second Gaia data release of an apparent merger remnant, dubbed Gaia-Sausage/Enceladus (GS/E), which dominates the nearby stellar halo (Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018). While GS/E has been the subject of intense recent study, its synthesis can be traced back before the second Gaia data release to the notion of a ‘dual’ galactic halo based on its broken radial density profile and bimodal chemistry (e.g. Carollo et al. 2007; Nissen & Schuster 2010; Deason, Belokurov & Evans 2011; Bonaca et al. 2017; Hayes et al. 2018). GS/E is characterized by a large group of stars on highly radial orbits, with near-zero angular momentum L_z, a wide span of radial velocities, and high eccentricities. The remnant was quickly seen to comprise a significant fraction (≈ 50 per cent) of the density of the nearby stellar halo (Belokurov et al. 2018; Fattahi et al. 2019; Lancaster et al. 2019), and initial estimates of the stellar mass of the progenitor were quite high, ranging from 6 × 10⁸ to nearly 10¹⁰ M_⊙ (Belokurov et al. 2018; Helmi et al. 2018; Fattahi et al. 2019; Myeong et al. 2019; Vincenzo et al. 2019; Das, Hawkins & Jofré 2020; Feuillet et al. 2020). The simplest interpretation of this accretion remnant is that it was deposited during a massive, head-on merger event early in the life of the Milky Way (Belokurov et al. 2018; Helmi et al. 2018; Deason, Belokurov & Sanders 2019; Fattahi et al. 2019; Iorio & Belokurov 2019; Mackereth et al. 2019a). Although other works have argued that the observed chemodynamics are also consistent with multiple smaller accretion events (Donlon Thomas et al. 2022; Donlon & Newberg 2023).

Recent, thorough studies of the GS/E accretion remnant focusing on its density profile have tended to settle on a lower stellar mass for the progenitor in the range 3 × 10⁸ to nearly 10⁹ M_⊙ (Mackereth et al. 2019a; Mackereth & Bovy 2020; Naidu et al. 2021; Han et al. 2022). Additionally, the chemical composition of the remnant has received attention, and specifically its [Fe/H] abundance distribution (which extends to high [Fe/H] ∼ −0.6) as well as its [α/Fe] distribution, both of which appear to reflect those of a massive dwarf galaxy (Myeong et al. 2019; Mackereth et al. 2019a; Feuillet et al. 2020; Monty et al. 2020; Hasselquist et al. 2021; Matsuno et al. 2021; Buder et al. 2022; Gaia Collaboration et al. 2023; Horta et al. 2023a), supporting the earlier findings of an [Fe/H]-rich structure by Belokurov et al. (2018) as well as pre-Gaia second data release studies (Nissen & Schuster 2010; Bonaca et al. 2017; Hayes et al. 2018). The ages of likely GS/E stars were measured by Montalbán et al. (2020) and found to be typically 10 Gyr, supporting other works which have suggested a similar time frame for the occurrence of the accretion episode (z ≈ 2). Our understanding of the kinematics of the remnant have been extended by Lancaster et al. (2019) and Iorio & Belokurov (2021), who measure a strong degree of radial anisotropy in the population, behaviour extending over a wide range of Galactocentric radii. Naidu et al. (2021) and Chandra et al. (2022) compare H3 data with tailored N-body simulations and develop a coherent narrative for the specifics of the infall of the progenitor. They also link several other halo structures such as Arjuna (Naidu et al. 2020), the Hercules-Aquila cloud, and the Virgo Overdensity with GS/E debris. The association of GS/E with the latter two structures was first proposed by Simion, Belokurov & Koposov (2019), with support coming from the analysis of stellar halo density profiles by Iorio & Belokurov (2019), and has recently been affirmed through chemical analysis by Perottoni et al. (2022). Given the menagerie of recently discovered halo substructures (see Helmi 2020; Naidu et al. 2021; Horta et al. 2023a, for censuses of the major ones), there are open questions regarding which are unique and independent of GS/E, and which are associated with, or simply complex dynamical echoes of, GS/E.

Recently, Lane, Bovy & Mackereth (2022, hereafter LBM22) presented kinematic models of the nearby stellar halo based on distribution functions representing the major constituent populations: the metal-rich, radially anisotropic population attributed to GS/E, and the metal-poor, comparatively radially isotropic population attributed to the (in situ) remainder of the stellar halo (e.g. Belokurov et al. 2018; Haywood et al. 2018). LBM22 found that kinematic selection criteria commonly used to identify GS/E stars are subject to bias in the context of their models, and that most selection criteria are at best about 70 per cent, and at worst about 50 per cent pure (i.e. only 5–7 in every 10 stars attributed to GS/E are genuine members of the GS/E remnant). They also use their models to construct high-purity, scale-free (i.e. they are resilient to changes in stellar sample or gravitational potential used for calculation of kinematic quantities) selections for GS/E, which reach as high as 85 per cent purity. LBM22 highlight that while the GS/E remnant is an overdensity occupying a characteristic region of phase-space, it is does not do so without contamination from other stellar populations. They emphasize the need to account for these effects when modelling the GS/E remnant, which is particularly important as studies seek to place tighter constraints on its properties.

In this study we will build on the work of Mackereth & Bovy (2020, hereafter MB20) who studied the density profile of the stellar halo with data from Gaia and the Apache Point Galactic Evolution Experiment (APOGEE Majewski et al. 2017). We integrate the kinematic models of LBM22 into the density modelling framework of MB20 so that we may select a high-purity sample of GS/E stars based on kinematics and directly assess the density profile and stellar mass of the remnant. Our novel method, using distribution function models in concert with density modelling, allows us to study the GS/E remnant unhindered by contaminated kinematic selections, and represents an important step towards future 6D modelling of the entirety of the stellar halo phase space.

This paper is laid out as follows: in Section 2, we describe our APOGEE and Gaia observations and relevant sample selection procedures used to isolate high-purity samples of GS/E stars. In Section 3, we describe the density modelling framework, paying specific attention to the integration of distribution function models. We also thoroughly validate our new approach using mock data. In Section 4, we present the results of our analysis of the GS/E remnant, including fits to the whole stellar halo for comparison. In Section 5, we discuss results in the context of recent literature, and also address the limitations of our new modelling technique. We end in Section 6 by summarizing our results and looking ahead to future applications of distribution functions to studies of the stellar halo.

2 DATA

2.1 Observations and stellar properties

We use data from APOGEE (Majewski et al. 2017), specifically the sixteenth data release (DR16 Ahumada et al. 2020), a component of the Sloan Digital Sky Survey IV (SDSS; Blanton et al. 2017). We do not use the most recent APOGEE data as we prefer to use the same data as LBM22. APOGEE-1 (an SDSS III program; Eisenstein et al. 2011) and its successor APOGEE-2 (an SDSS IV program) obtain high-resolution (⁠|$R \approx 22\, 500$|⁠), near-infrared (⁠|$1.51 \, \mu \mathrm{m} \lt \lambda \lt 1.7 \, \mu \mathrm{m}$|⁠), high signal-to-noise (>100 per pixel) spectra using two near-identical spectrographs (Wilson et al. 2019) at the 2.5 m Sloan Telescope (Gunn et al. 2006) at the Apache Point Observatory in New Mexico (APOGEE and APOGEE-2N) and the 2.5 m Iréné du Pont Telescope (Bowen & Vaughan 1973) at the Las Campanas Observatory in Chile (APOGEE-2S). APOGEE DR16 contains observations of approximately 430 000 unique stars from both hemispheres, which span a range in Galactocentric radii and heights above and below the disc, including all major stellar populations of the Milky Way.

APOGEE targets are selected using 2MASS (Skrutskie et al. 2006) photometry by a two-fold strategy, including stars targeted specifically for focused science goals and ancillary programs as well as the systematic observations of the ‘main’ survey. Stars in the main survey are chosen based on their reddening-corrected (Majewski, Zasowski & Nidever 2011) (J − K_S)₀ colour and H-band magnitude from 2MASS. The targeting procedure is described by Zasowski et al. (2013) for APOGEE-1, and built upon by Zasowski et al. (2017), Beaton et al. (2021), and Santana et al. (2021) for APOGEE-2. Data are reduced using a pipeline originally described by Nidever et al. (2015) but which has been built upon for recent data releases (Holtzman et al. 2018; Jönsson et al. 2020). Throughout this paper we use stellar elemental abundances, atmospheric parameters, and radial velocities from the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP, García Pérez et al. 2016). The pipeline uses a custom line list (Shetrone et al. 2015) which has been updated for DR16 (Smith et al. 2021), and a custom stellar spectral library (Zamora et al. 2015) which has also been updated for DR16 (Jönsson et al. 2020).

To complement the APOGEE data we use astrometry from the third data release (DR3; Gaia Collaboration et al. 2022) of the Gaia space telescope (Gaia Collaboration et al. 2016). We match our APOGEE data to the Gaia DR3 catalogue using the CDS X-match service¹ implemented in the gaia_tools² python package. We specifically rely on Gaia proper motions, as APOGEE provides more accurate radial velocities. We further use spectro-photometric distance estimates determined using the astroNN artificial neural network framework³ (Leung & Bovy 2019a, b). We specifically take those distances from the SDSS DR17 astroNN Value Added Catalogue,⁴ which are obtained from a model trained on newer APOGEE DR17 and Gaia eDR3 data (see Abdurro’uf et al. 2022). The Bayesian artificial neural network predicts stellar luminosity using stellar spectra, and is trained on APOGEE spectra and Gaia parallaxes. It simultaneously predicts distances and accounts for the Gaia parallax zero-point. Distance uncertainties are typically less than 10 per cent, much better than Gaia parallaxes for stars further than a few kpc from the Sun.

We extract a high-quality subsample of APOGEE red giants from the broader data set. First off, we only consider stars that are part of the APOGEE statistical sample, which are stars chosen randomly according to their (J − K_S)₀ colour and H-band magnitude as described above (i.e. they are not part of some ancillary science program, or specifically targeted for some other reason). Stars from this group are chosen which have relative distance uncertainty, σ_d/d, less than 20 per cent, surface gravity between 1 < log g < 3, and surface gravity uncertainty less than 0.1 dex. These giants must also have measured Gaia proper motions and positions, they must have AstroNN distances, and they must have measured APOGEE radial velocities, [Fe/H], and [Al/Fe] abundances (Note that while we will show other abundances, only these are crucial to our analysis). We also remove stars lying in fields near the Galactic bulge, and those containing stars from globular cluster so that our sample contains only stars from the smooth component of the stellar halo. We enforce the first of these criteria by removing fields within −20° < ℓ < 20° and |b| < 20°. The second of these criteria is met by first identifying any fields near the on-sky location of a known globular cluster from the Harris (1996, December 2010 version) catalogue. For each of these fields, we examine the [Fe/H] and radial velocity distributions of the stellar constituents by eye to see if there is a noticeable enhancement of stars with properties similar to those of the cluster, paying specific attention to the region within two tidal radii of the cluster centre. Fields with a distinct group of stars with properties matching the cluster are entirely removed from consideration. Throughout the rest of the paper, we refer to this subset of APOGEE stars from the statistical sample as the ‘quality sample’; it contains 98 245 stars. The distribution of Galactocentric radii of this quality sample is shown in Fig. 1. The distribution is concentrated near the location of the Sun, about 8 kpc, and extends from 2 kpc to nearly 40 kpc.

We calculate kinematic properties of the stars in our sample using galpy⁵ (Bovy 2015). Throughout, we assume that the Sun lies at Galactocentric coordinates (R, z, ϕ)  = (8.275 kpc, 0.028 kpc, 0) (Bennett & Bovy 2019; GRAVITY Collaboration et al. 2021), the peculiar solar motion is (U, V, W) = (11.1, 12.24, 7.25) km s⁻¹ (Schönrich, Binney & Dehnen 2010), and that the circular velocity at the location of the Sun is 220 km s⁻¹. We use a left-handed coordinate system such that angular momenta are positive in the direction of Galactic rotation. We use the MWPotential2014 potential defined in Bovy (2015), and calculate energy, actions, eccentricity, and orbital extrema for each star. Actions are calculated using the ‘Stäckel fudge’ method of Binney (2012) whereby the Milky Way is locally approximated as a Stäckel potential using a focal length estimated with equation (9) of Sanders (2012). Eccentricities are calculated using a method similar to Stäckel fudge presented in Mackereth & Bovy (2018). We principally consider three kinematic spaces of the many which are commonly used in the study of the stellar halo. First is eccentricity as a function of z-axis angular momentum (e − L_z). Second is the square root of the radial action versus z-axis angular momentum (⁠|$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠). Third is the action diamond (ad), which has the difference between the radial and vertical actions as a function of the z-axis angular momentum, but both quantities are normalized by the absolute sum of all three actions which leads to the characteristic diamond boundary in the space. For more information on the action diamond see Section 3 of LBM22 and references therein.

2.2 The halo and GS/E samples

From the quality sample, we identify likely halo stars according to their abundances. Broadly following MB20 and LBM22, we focus on the iron abundance, and specifically those stars with [Fe/H] < −1. The principle issue with this simple selection, however, is that GS/E is known to have [Fe/H] up to ∼−0.6 (e.g. Myeong et al. 2019; Monty et al. 2020; Hasselquist et al. 2021; Horta et al. 2023a). We therefore modify our selection to be more inclusive of GS/E, increasing our maximum [Fe/H] to −0.6. This introduces a different issue though, which is that a sample selected with this [Fe/H] limit would have substantial contamination from thin and thick disc populations. To mitigate this effect, we employ aluminium abundances which have been shown to reliably separate accreted halo and Galactic disc populations (Hawkins et al. 2015). We take specific inspiration from Belokurov & Kravtsov (2022) in constructing our [Al/Fe]-based selection (who also provides an illuminating discussion on the use of [Al/Fe] for this purpose), and note that several other recent works have employed similar selections, or demonstrated the separation of disc and accreted halo components on the basis of [Al/Fe] (Das et al. 2020; Hasselquist et al. 2021; Horta et al. 2023a). Our selection is piecewise-continuous in the [Al/Fe] versus [Fe/H] plane. At [Al/Fe] > −0.1, we define the halo as [Fe/H] < −1. At [Al/Fe] = −0.1, the approximate boundary between accreted and in-situ populations (Hawkins et al. 2015; Das et al. 2020; Hasselquist et al. 2021), our selection increases in [Fe/H] to the locus ([Fe/H],[Al/Fe]) =(− 0.6, −0.3). Our selection then includes [Fe/H] < −0.6 for [Al/Fe] < −0.3.

This selection (dashed line), and the 1523 halo stars from the quality sample chosen with it, are shown in the left panel of Fig. 2. The radial range of this sample is also shown as a dashed line in Fig. 1; it extends about as far as the quality sample. In Fig. 2, we display the stars defined as part of the halo sample in other abundance planes commonly used in the literature. The middle panel shows [Fe/H] versus [Mg/Fe], showing a large group of stars at high [Mg/Fe] and low [Fe/H] traditionally attributed to the in-situ halo and thick disc, accompanied by an extended group at lower [Mg/Fe] and higher [Fe/H] which is emblematic of GS/E and other accreted populations (e.g. Hasselquist et al. 2021; Horta et al. 2023a). The solid lines are fiducials from other studies, the vertical line marks [Fe/H] = −1, which is the cut-off adopted by MB20, and the sloped line is used by Mackereth et al. (2019a) and LBM22 to separate disc and halo populations. It is clear that these cuts, which rely predominantly on [Fe/H], would leave much of GS/E excluded from the sample. The right panel shows [Al/Fe] versus [Mg/Mn], which is used by some authors to separate accreted and in-situ populations (e.g. Das et al. 2020; Horta et al. 2021a; Fernandes et al. 2023). The upper left region of the diagram is where unevolved populations such as those in dwarfs or in the old stellar halo live. The two other regions are where high and low-α evolved, in-situ populations tend to live. Here we can see that the majority of our sample occupies the unevolved region, but some is found in the high-α evolved region, much of which are likely thick disc stars.

From the halo sample, we identify a subset of likely members of the GS/E remnant following the recommendations and approach of LBM22. Those authors built DF-based models of the nearby Milky Way stellar halo using the anisotropy parameter β as a proxy for the two major stellar halo populations: β = 0.9 for the high-[Fe/H], radially anisotropic GS/E remnant, and β = 0.5 for the low-[Fe/H], more isotropic traditional stellar halo. Using these models, LBM22 identified regions of various kinematic spaces where stars from the high-β component were found with low contamination from the low-β component. The ‘best’ kinematic spaces were those where purity was around ∼0.8, and these included e − L_z, the action diamond, and |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠. We use the elliptical selection criteria from table 2 (The Survey selections) of LBM22, defining a separate GS/E subsample for each of these three kinematic spaces. Fig. 3 shows these three spaces, the GS/E selections, as well as the chemically selected halo sample coloured by [Fe/H]. Our GS/E subsamples defined using the e − L_z plane, action diamond, and |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| plane contain 163, 77, and 319 stars respectively.

We consider three kinematic spaces here so that we may examine the differences in outcomes of our model-fitting procedure. This type of selection constitutes a significant bias when modelling GS/E, in so far as we discount a large fraction of stars belonging to the high-β GS/E population where they overlap with the low-β halo. To wit, completeness for these selections is typically between 0.3 and 0.7 and, as we shall show, is spatially non-trivial. We return to these DF-based models in Section 3.2 when we discuss the creation of the kinematic effective selection function which addresses these biases in the density modelling framework.

The final step in cultivating our subsamples of GS/E stars for modelling is to address contamination from thick disc populations. LBM22 include a thin and thick disc component in their models, but acknowledge that this simplifies the complicated underlying nature of the Milky Way disc, which is closer to a continuum of disc populations parametrized by scale height, length, and velocity dispersion (Bovy et al. 2012). The selections of LBM22 find near-perfect separation between thick disc and high-β populations, but their simple model is unable to account for the kinematically hottest thick disc populations. Even though these populations are numerically insignificant in the context of the density of the stellar disc, they are significant in the context of the density of the stellar halo.

Fortunately, these populations should be somewhat distinct in abundance space. Fig. 4 shows the GS/E subsamples in [Fe/H] versus [Al/Fe] (top row) and [Mg/Fe] (bottom row). Since in-situ populations (both disc and stellar halo) should occupy the higher [Al/Fe] portion of our halo selection region, we make a cut at [Al/Fe] = −0.1 and only consider kinematically selected stars with lower [Al/Fe] to be part of our GS/E subsamples. In Fig. 4, the horizontal dashed line is set at [Al/Fe] = −0.1, and stars with greater [Al/Fe] abundances are shown as black crosses in both sets of panels. In [Fe/H] versus [Mg/Fe], many of the stars which are discounted occupy the high-[Mg/Fe], intermediate-[Fe/H] part of the distribution, exactly where thick disc and in-situ halo would expect to be found. We are therefore confident that this cut serves only to enhance the purity of our GS/E subsamples and does not confer any bias. The fraction of stars removed is 5–8 per cent of the sample depending on the kinematic space. The final GS/E subsamples defined by the kinematic selections, including this cut in [Al/Fe] contain 151 (e − L_z), 73 (ad), and 300 (⁠|$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠) stars.

The top and two rightmost panels in Fig. 4 show the marginalized [Fe/H], [Al/Fe] and [Mg/Fe] abundance distributions respectively. The coloured lines are Gaussian kernel density estimates for each of the three kinematic selections and the grey line shows the same for the whole sample. Dashed lines show the medians in the same colours, which are all similar, deviating by no more than 0.02 dex for each abundance. These medians are approximately −1.17, −0.23, and 0.20 for [Fe/H], [Al/Fe], and [Mg/Fe] respectively.

3 METHOD

In this section we describe our approach to modelling the density of both the whole stellar halo and the GS/E subsamples. We describe the manner in which we account for the biases imposed on our samples by the APOGEE selection function and by the high-purity GS/E kinematic selections. We then outline how density profiles are normalized and total masses are derived. Finally, we demonstrate that our technique is robust and yields accurate results by applying it to mock data.

3.1 Density modelling

We fit density profiles to both our halo and GS/E samples, accounting for the survey selection function, dust extinction, and the colour-magnitude density of the tracer population. We generally follow the well-established approach presented by Bovy et al. (2012) (also discussed by Rix & Bovy 2013), and then built upon and used more recently by Bovy et al. (2016a, b); Mackereth et al. (2017); Horta et al. (2021b) and MB20. We describe here the method in broad detail, and encourage the reader to refer specifically to Bovy et al. (2012, 2016a) and MB20 for more information.

We model the observed APOGEE star counts as a Poisson point process, with rate function λ(O|Θ). Here, the data O ≡ [ℓ, b, D, μ_{[ℓ, b]}, v_LOS, H, (J − K_S)₀, [Fe/H]] are the Galactic coordinates ℓ and b, heliocentric distance D, proper motions in Galactic coordinates μ_{[ℓ, b]}, the heliocentric line-of-sight velocity v_LOS, the apparent H-band magnitude H, the dereddened colour (J − K_S)₀, and the iron abundance [Fe/H]. The vector Θ of parameters describes the density profile under consideration and is to be determined by fitting to the data. The rate function can be expressed as

where |$\nu _{\star }(\vec{X} \vert \Theta)$| is the stellar number density as a function of rectangular Galactocentric coordinates |$\vec{X}$| and |$f(\vec{v} \vert \vec{X})$| is the velocity distribution function to which we will return in the next section when we come to the kinematic selection function. While the velocity DF in general depends on parameters Θ fit to the data, here we use a fixed velocity DF |$f(\vec{v}\vert \vec{X})$| to keep the fit computationally tractable. We address the implication of this choice later in the discussion. The function |$\rho (H, (J-K_\mathrm{S})_{0}, [\mathrm{Fe/H}] \vert \vec{X})$| is the apparent density of stars in colour, magnitude, and abundance space for a given position. |$\vert J(\vec{X}, \vec{v}; \ell , b, D, \mu _{[\ell ,b]}, v_\mathrm{LOS}) \vert$| is the Jacobian linking observed coordinates with Galactocentric coordinates, which may be split into a spatial and velocity component as |$\vert J_{1}(\vec{X}; \ell , b, D) J_{2}(\vec{v}; \ell , b, D \mu _{[\ell ,b]}, v_\mathrm{LOS}) \vert$|⁠. Finally, S(ℓ, b, D, μ_{[ℓ, b]}, v_LOS, H, (J − K_S)₀) is the selection function, which depends on sky location, heliocentric distance, velocity, as well as colour and magnitude. This single selection function can be separated into a spatial component, S₁, and a kinematic component, S₂, as S(ℓ, b, D, μ_{[ℓ, b]}, v_LOS, H, (J − K_S0) = S₁(ℓ, b, H, (J − K_S)₀) × S₂(ℓ, b, D, μ_{[ℓ, b]}, v_LOS). The spatial component arises from the APOGEE targeting procedure, which has been briefly discussed in Section 2. This function describes the probability with which a star with a given H-band magnitude and (J − K_S)₀ colour will be included for observation in the APOGEE statistical sample, and varies based on field. We follow the approach of MB20 when constructing S₁, and details may be found in their appendix A. The kinematic component of the selection function arises from our selection of GS/E stars, and its construction will be described below.

As in Bovy et al. (2012), the log likelihood of the model parameters may be written in terms of the rate function as

Derivation of this likelihood uses the fact that the rate function only depends on Θ through |$\nu _{\star }(\vec{X} \vert \Theta)$|⁠, because as we discussed above, the velocity DF does not depend on parameters that we fit. If the velocity DF were to depend on Θ as well, the first term would simply be |$\ln [\nu _{\star }(\vec{X}_{i} \vert \Theta)\, f(\vec{v}_{i} \vert \vec{X}_i,\Theta)]$|⁠. The integral in the likelihood equation is the effective volume, which expresses for a given set of model parameters the expected (non-normalized) number of stars. For APOGEE, the effective volume can be expressed as

Here, the effective volume is a sum over each APOGEE field of area Ω_field under consideration. This is permitted for a survey such as APOGEE which is comprised of a number of small (<2 degree size), non-overlapping fields, since we can assume that the density is constant across the angular extent of each field (an assumption which we test, finding it to be valid at level of about 1 per cent). Recall in Section 2.1 we removed any APOGEE field with noticeable contamination from a globular cluster, which helps to enforce this assumption. In the integral the density is as above, but the positions |$\vec{X}$| are evaluated along the line-of-sight of each field. |$\mathfrak {S}_{1}$| and |$\mathfrak {S}_{2}$| are the spatial and kinematic effective selection functions. Similar to S₁ and S₂, the effective selection functions include additional information required for the modelling procedure, yet which is independent of the model parameters of interest Θ.

The spatial effective selection function |$\mathfrak {S}_{1}$| is determined by marginalizing S₁ over the colour, absolute magnitude, and metallicity distribution of the tracer population, including the effects of dust obscuration. It is given by

Here ρ(M_H, (J − K_S)₀, [Fe/H]) and S₁(field, H, (J − K_S)₀) are as above. We have introduced the absolute magnitude M_H = H − μ(D), which relates to the apparent magnitude by the distance modulus μ(D), and is the preferred quantity for this integration. In the main APOGEE survey each field is divided up into j bins of apparent H-band magnitude and (J − K_S)₀, which are each bounded by a minimum and maximum magnitude H_{[min, max]}. The factor Ω_j/Ω_field is the fraction of observable area of the plate given the local extinction and the limiting H-band magnitudes for each bin. Ω_field is the total area of each plate, and Ω_j is given by

where A_H(ℓ, b, D) is the H-band extinction. Intuitively, as extinction increases stars may be extinguished into or out of a given colour-magnitude bin, and the factor Ω_j expresses this effect. Note that since (J − K_S)₀ is the dereddened colour we assume it is not impacted by extinction and so does not appear in the definition of Ω_j. The H-band extinction is obtained from the dust map of Bovy et al. (2016a), which is a combination of the dust maps from Marshall et al. (2006), Green et al. (2015), and Drimmel, Cabrera-Lavers & López-Corredoira (2003). We query this dust map using the mwdust python package.⁶

In practice, this spatial effective selection function integral is approximated by Monte Carlo sampling the distribution of colours, absolute magnitudes, and iron abundances of the red giant tracer population. These distributions are sourced from a grid of PARSEC v1.2S isochrones (Bressan et al. 2012). The grid is spaced linearly in metal fraction with ΔZ = 0.0001 and a minimum [Fe/H] of −3. The ages of the grid span 10–14 Gyr, with the minimum age being appropriate for GS/E (Montalbán et al. 2020). We draw samples from regions of the isochrone grid bounded by 1 < log g < 3 weighted by a Chabrier (2003) initial mass function. The contributions from each isochrone are weighted such that the samples are spaced linearly in [Fe/H] and linearly in age. When fitting the whole-halo sample we use a [Fe/H] range of [−3, −0.6], and when fitting the GS/E subsamples we use a slightly more metal-rich range: [−2, −0.6] to reflect approximate spread in [Fe/H] observed for GS/E (see our Fig. 4 or refer to e.g. Myeong et al. 2019; Hasselquist et al. 2021; Horta et al. 2023a). Otherwise, the fits to the halo and GS/E samples both use the same isochrone grid. The effective selection function is calculated on a grid for each APOGEE field with 300 bins between 7 < μ(D) < 19 (roughly 0.25 kpc < D < 65 kpc) spaced linearly in distance modulus.

3.2 The kinematic effective selection function

The kinematic effective selection function can be expressed as:

where f, J₂, and S₂ are all as above, but defined at field locations (ℓ, b, D). As discussed by MB20, the ideal approach for handling stellar kinematics is to have a parametrized DF |$f(\vec{v} \vert \vec{X}, \Theta)$| which may be used to kinematically select populations with a set of rules S₂. The downside of this approach is twofold: first, DFs – particularly realistic ones – tend to be computationally expensive to compute, and second, this expands the coordinate space to 6D. These two effects combine such that it is impractical to use DFs with varying parameters which must be integrated on the fly during the optimization of the likelihood.

MB20 bypassed this issue by asserting that the two major halo stellar populations: the radially anisotropic, high-[Fe/H] GS/E remnant, and the more isotropic, low-[Fe/H] traditional stellar halo are well-separated by eccentricity. In this work, we adopt a more sophisticated approach, leveraging the DF-based models of LBM22, which we have already discussed in the context of GS/E sample selection. We will use this model to construct the kinematic effective selection function in such a way that it is independent of any model parameters Θ.

To briefly summarize the model of LBM22, they consider anisotropic distribution functions of the form (see, e.g. Binney & Tremaine 2008)

Here, |$\mathcal {E} = \Psi -\frac{1}{2}v^{2}$| is the binding energy and Ψ = −Φ + Φ₀ is the negative of the gravitational potential normalized such that Ψ(∞) = 0. L is the total angular momentum, v is the velocity magnitude, and β is the anisotropy parameter, defined as

where σ_T is the quadrature sum of the azimuthal and polar velocity dispersions, and σ_r is the velocity dispersion in the radial direction. The function f₁ relates the mass density of the DF to the underlying potential and is, in general, non-trivial to compute, often requiring solving multiple integrals numerically (For more information see LBM22, particularly Appendix A).

We use the same underlying potential (MWPotential2014 of Bovy 2015) and mass density for the DF as LBM22. The mass density is a spherical truncated power law with power-law index α = 3.5 and truncation radius r_c = 25 kpc, taken from the stellar halo fits of MB20. We use β = 0.3 to represent the low-[Fe/H] traditional stellar halo, and β = 0.8 to represent the high-[Fe/H] GS/E halo. We modify our choice of high- and low-β compared with LBM22, who use β = 0.9 and β = 0.5. A smaller value for the low-β component aligns slightly better with the accepted values for the metal-poor halo from the literature (Belokurov et al. 2018; Fattahi et al. 2019; Lancaster et al. 2019; Iorio & Belokurov 2021). With regards to the anisotropy of the high-β component, the literature does support a value of β ∼ 0.9 (Belokurov et al. 2018; Fattahi et al. 2019; Lancaster et al. 2019; Iorio & Belokurov 2021), However we prefer to be conservative with our choice of β = 0.8 since LBM22 did note some minor numerical issues when using DFs with higher values of β. Additionally, as we shall show, by modelling GS/E with β = 0.8 we are actually constructing an upper limit on the mass when compared with a mass derived using a higher value of β.

We set these two components to have equal density in our models in order to enforce the assumption that the two components exist in roughly equal density near the position of the Sun. While it may seem odd that we are fixing the parameters of the mass density distribution of the stellar halo in our DF models before setting out to use them to measure the mass density distribution of the stellar halo, we are confident of a scheme to ensure that this does not significantly bias our results, and will discuss more in Section 5. It is worthwhile to note though that only the high-β component is relevant for calculating the kinematic effective selection function, the low-β component is only used to estimate purity, which is simply informative. Additionally, the overall normalization of this high-β DF is also unimportant, only its shape, which is determined by the value of β as well as the form of the underlying density profile. These two factors contribute to the reasoning why the kinematic effective selection function is resilient to changes in the assumptions of our models, which we will expand upon later.

LBM22 studied samples drawn from these DFs at positions where APOGEE had observed stars. We adopt a similar approach when computing |$\mathfrak {S}_{2}(\mathrm{field},D)$|⁠. At each location where we calculate |$\mathfrak {S}_{2}$|⁠, we draw 10³ velocity samples from our DFs and calculate actions and eccentricities in a manner identical to that used for our stellar samples (Section 2). We then calculate the completeness of the high-β component given each of the kinematic selections used in Section 2.2 (i.e. the fraction of high-β samples lying in the selection compared with the total number of high-β samples), which is equivalent the value of |$\mathfrak {S}_{2}$|⁠. Drawing samples from DFs and calculating kinematic quantities is computationally expensive, and so rather than compute |$\mathfrak {S}_{2}$| for each point on a grid identically to the grid on which we computed |$\mathfrak {S}_{1}$|⁠, we compute it only for 21 points, evenly spanning the same range of distance modulus (7 to 19), for each field. We then linearly interpolate the value of |$\mathfrak {S}_{2}$| between these 21 points when combining it with the more densely sampled |$\mathfrak {S}_{1}$|⁠.

Fig.5 shows the value of |$\mathfrak {S}_{2}$| as a function of distance modulus for each of the three kinematic selections. Each field is shown as a different line, coloured by the absolute value of the Galactic latitude of its on-sky location. Only fields with stars included in the halo sample are shown. In a smaller panel at the top of the figure we show the distribution of distance moduli of our halo sample for reference. The value of |$\mathfrak {S}_{2}$| fluctuates along a field due to the finite number of samples drawn at each location, but this noise clearly averages out in the agregate of many fields. Fig. 6 shows the value of |$\mathfrak {S}_{2}$| marginalized along each line of sight using a spline fit to the distribution of distance moduli in the top panel of Fig. 5 (red line) as a function of Galactic coordinates. Again, only fields which have stars in the halo sample are shown, and the 20° box we use to exclude bulge fields is shown.

These two figures each reveal an interesting perspective of |$\mathfrak {S}_{2}$|⁠. Fig. 5 shows how it can vary along each line of sight, and similarly Fig. 6 shows how it can vary as a function of sky coordinates. These variations are principally due to the change in kinematics as a function of radius, which when mapped into a heliocentric frame and weighted by the distribution of halo star distance moduli yields interesting patterns. A secondary effect is that samples drawn near the Galactic poles will behave differently than those at a similar radius in the plane. Directly above the centre of the Galaxy, the Stäeckel approximation is a poor one, impacting actions and eccentricities, and also the kinematic selections are ill-suited to describe the distributions which arise. While this effect is somewhat pathological, it is at least entirely self-consistent within the modelling framework, such that the data and |$\mathfrak {S}_{2}$| are impacted in the same way. This gives rise to the general trend of worsening |$\mathfrak {S}_{2}$| with larger absolute Galactic latitude seen in both Figs 5 and 6. For both e − L_z and ad, |$\mathfrak {S}_{2}$| is reasonably constant at distances characteristic of our sample (μ ∼ 10–17), with the decreases in |$\mathfrak {S}_{2}$| associated with Galactocentric radius and latitude occurring at slightly lower μ for the ad space. On the other hand, |$\mathfrak {S}_{2}$| has a much more complicated behaviour on a field-by-field basis for |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠.

3.3 Stellar number density models

We consider several models for the stellar number density, which will be fit both to the whole stellar halo sample and the GS/E subsamples. We follow the approach of MB20 in our construction of these density profiles, which are general and flexible. All of our models are triaxial and expressed as functions of an effective radius.

where [X, Y, Z]′ are rectangular coordinates in the rotated frame. p is the ratio of the Y′ to X′ axes, and q is the ratio of the Z′ to X′ axes. Each model has a flexible orientation parametrized by three variables: θ, η, and ϕ. First, the density profile is rotated in the X′–Y′ plane clockwise by an angle ϕ, then the Z′ axis is rotated towards a unit vector |$\hat{z} = \big [\sqrt{ 1-\eta ^{2} } \cos \theta , \sqrt{ 1-\eta ^{2} } \sin \theta , \eta \ \big ]$|⁠. So the parameter θ controls the angle of the transformed Z′-axis of the ellipsoid (i.e. the axis of the ellipsoid originally aligned with the Z′ axis) in the plane of the Galaxy, and η controls the degree to which the transformed Z′-axis is oriented towards the Galactic north pole. Even sampling of these three parameters evenly samples the unit sphere.

The density profiles we consider are as follows. The first is a simple power-law profile

where α₁ is the power-law index. The second is a single power-law profile with an exponential truncation scale r_c, which is expressed as

Finally, we have broken power-law models with two or three break radii which take the form

where r_i are the break radii (we use r to emphasize that these are radii), and α_i are the power-law indices. Together these four models span a range of complexity, and are similar to others used in the literature. Throughout the rest of the paper we refer to the single power-law model as ‘SPL’, the exponentially truncated single power-law model as ‘SC’, and the broken power-law models with one or two breaks as ‘BPL’ and ‘DBPL’, respectively.

Each of these models is considered on its own, but also with a disc contamination model. While we are confident that we have removed most of the disc contamination using abundance cuts, we recognize that some contamination is unavoidable. The disc contamination model is exponential, with radial scale length 2.2 kpc and vertical scale length 0.8 kpc (Mackereth et al. 2017). We parametrize contamination by including a factor f_disc which expresses the fraction of density at the solar position contributed by the disc model. The combined model takes the form

Models including disc contamination are referred to with ‘+D’ after the shorthand introduced above (i.e. ‘SPL+D’).

In fitting these models, power-law indices are allowed to take any value and are nominally positive (indicating decreasing density), but occasionally the exponentially truncated models have negative α. Indices are also required to steepen in broken power-law models such that α₁ < α₂ < α₃. Truncation radii for both exponentially truncated and broken power laws are required to lie in the range [2,55] kpc. The inner range is motivated by the fact that we exclude any fields within a 20° box around the bulge (20° at ∼8 kpc is ∼2 kpc). The outer range is driven by the fact that the effective selection function essentially becomes zero for all fields at this distance. Even though our most distant star lies at nearly 40 kpc, our statistical leverage extends to the point at which the effective selection function is zero, which is about 55 kpc. The shape parameters p and q are defined on the interval (0,1] such that the axis of the ellipsoid aligned with the X′ axis (rotated frame) is always the longest. Of the orientation parameters, θ can take values between [0, 2π], η can take values between [0.5,1], and ϕ is defined between [0, π]. η and ϕ could take wider ranges but these are degenerate with other combinations of orientation parameters as well as p and q. Finally, the disc contamination parameter f_disc is defined on the interval [0,1].

3.4 Fitting density models

To find the best-fitting set of parameters for each model, we use Markov Chain Monte Carlo to sample the posterior of the likelihood (2). We use the affine-invariant ensemble sampler of Goodman & Weare (2010) implemented in emcee (Foreman-Mackey et al. 2013). We first optimize the likelihood function using Powell’s method (Powell 1964) and then initialize one hundred walkers with a small amount of noise around that set of optimized parameters. Each walker moves 10⁴ steps and we remove the first 10³ steps as burn-in. We define the best-fitting parameters as the median of this set of posterior samples, with the uncertainty being the 16th and 84th percentiles. Our posteriors are complicated, however, and while we explored other more sophisticated methods of determining the best-fitting parameters and uncertainties we defer back to this approach for its simplicity. For this reason, we make our chains available online.⁷

3.5 The mass from normalization of the amplitude of the density function

To determine the mass of each density profile and its array of parameter samples drawn from the posterior, we follow the approach of MB20 which we describe here briefly. First, we normalize the density profiles from Section 3.3 and then integrate the density profiles over a range of Galactocentric radii. With regards to normalization, each density profile is fit such that the density is 1 at the solar position: |$\nu _{\star }(\vec{X}_{\odot }) = 1$|⁠. The number of stars expected for a given density profile and parameter vector Θ is therefore the rate function of equation (1), which is the product of the density and the effective selection functions, integrated over the volume of the survey. This quantity is identical to the effective volume expressed in equation (3). Dividing the total number of observed APOGEE stars to which the density profile is fit, N_obs, by this expected number of stars given |$\nu _{\star }(\vec{X}_{\odot }) = 1$| gives the proper normalization of the density profile in units of observed APOGEE red giants per volume, and is expressed as

We then convert this number density to a stellar mass-density using the same grid of PARSEC v1.2s isochrones (Bressan et al. 2012) used to calculate |$\mathfrak {S}_{1}$|⁠, again weighted by a Chabrier (2003) initial mass function. We calculate the average mass of giants 〈M_giants〉 observed in APOGEE by applying the equivalent of the observational cuts (from the APOGEE targeting procedure) in (J − K_S)₀ and our imposed cut between 1 < log (g) < 3 to a subset of isochrone grid defined by the minimum and maximum [Fe/H] for the sample of interest ([−3, −0.6] for the halo sample and [−2, −0.6] for the GS/E samples) and finding the initial mass function-weighted mean mass of the remaining points. Similarly, we find the fraction of stellar mass in giants, ω, by taking the ratio between the initial mass function-weighted sum of isochrone points within these boundaries and those outside, again for isochrones defined by the [Fe/H] range of the sample at hand. The conversion factor between giant number counts and total stellar mass is then simply calculated as

We compute this factor for each field, adjusting the limit in (J − K_S)₀ to reflect the minimum (J − K_S)₀ of the bluest bin adopted in that field. The edges in colour binning for each field are then accounted for by our integration under ρ(M_H, (J − K_S)₀, [Fe/H]) for the effective selection function. This factor ranges between 220 and 370 M_⊙ star⁻¹ depending on the colour and [Fe/H] limits. See appendix B of MB20 for a detailed validation of this approach using Hubble Space Telescope photometry. Combining these factors with the number density normalization, we attain the appropriate halo mass normalization at the Sun, |$\nu _{0} = A(\Theta)~ \chi \left(\mathrm{[Fe/H]_\mathrm{min,max}} \right)$| for a given sample defined by a range in [Fe/H] and a set of parameters from the posterior Θ. The now properly normalized density profiles are then integrated between Galactocentric radius 2 < r < 70 kpc. We choose this radial range to avoid over-extrapolating our fits, and to match the range used by MB20 for efficient comparison.

3.6 Validation of the method using mock data

Since we are introducing a new layer to the modelling approach outlined in the previous section by incorporating the kinematic effective selection function, it is prudent to test the method. We have created a pipeline to generate mock APOGEE data, including kinematics, for this purpose. We begin by considering any spherical potential expressible in galpy to act as the density profile representing the spatial distribution of the data. We draw mass samples from a Chabrier (2003) initial mass function such that the total mass of the samples is equal to the target total mass of the profile. We then assign these samples a radial position by sampling from the inverse cumulative mass function for the (spherical) density profile under consideration. Azimuthal and polar angles are assigned to the samples by using random sphere point picking. We scale the Y and Z positions by factors p and q respectively to make the spherical profile triaxial and rotate it by the same method described in Section 3.3, using a rotation by ϕ followed by a tilt towards |$\hat{z}$| defined by η and θ. In this way we generate a set of mass samples which have positions sourced from a density profile which can be made triaxial and rotated in the same manner as our candidate density profiles.

We assign stellar parameters and magnitudes to these samples with a single 10 Gyr old, Z = 0.001 ([Fe/H] ∼ −1.3) isochrone from the same PARSEC v1.2 grid described above using mass as the correspondence. We then remove all samples lying outside the APOGEE footprint, compute H-band extinction using the dust maps described above, and remove all stars with H-band magnitude and (J − K_S)₀ colours outside of the range considered by APOGEE on a field-by-field basis. For the remaining stars we apply the same selection procedure outlined in Zasowski et al. (2013) and Zasowski et al. (2017), randomly selecting stars in colour-magnitude bins with the appropriate probabilities to remain in the sample. We are then left with a sample of stars that would be seen by APOGEE as if it had ‘observed’ the underlying stellar population. The repository containing the code to create these mock data is available online.⁸

To test our model, we generate a mock from a stellar population with a total mass of 2 × 10⁸ M_⊙ and a single power-law density profile with the following parameters: α = 3.5, p = 0.8, q = 0.5, θ = π/4, |$\eta =1/\sqrt{2}$|⁠, and ϕ = π/6. We choose a low mass such that the generated mock has a small number of stars (645 for the whole mock, and 89–185 for the kinematically selected GS/E analogue samples). Then we may test whether our modelling approach can recover the input parameters when the number of observations is comparable to the number we see in our GS/E subsamples (∼100). To create analogue GS/E subsamples the whole mock is split in half, and each half is assigned kinematics from a DF with β = 0.8 and β = 0.3 respectively. We apply the kinematic cuts for each of the kinematic spaces considered in this work and then perform the fitting procedure described in Section 3.4 to the kinematically selected subsamples. We also fit the whole mock without considering kinematics, a fit which obviously does not include the kinematic effective selection function.

The results to these fits, including the total inferred mass, are shown in Fig. A1. We are able to recover most of the correct parameters of the underlying density profile within reasonable confidence intervals, and the inferred mass of the β = 0.8 population (which should be 10⁸ M_⊙, half the value of the total population) is also within the confidence interval. This not only demonstrates that our method works correctly, but that the biases induced by the kinematic selections we employ to isolate GS/E stars can be accounted for with the correct kinematic effective selection function. The masses for the GS/E subsamples are overestimated, with typical median values being about 10^8.1 M_⊙. The reason for this bias is that purity for the kinematic selections is not 100 per cent, but typically close to 80 per cent as per LBM22 (also see e.g. Limberg et al. 2022). In that work purity is determined for each kinematic selection criterion using kinematics for APOGEE DR16 stars assigned by the same distribution function models that we use in this work, and so the statistics should be applicable here. The impact of an impure kinematic selection is that stars from the mock low-β component are included in the number counts for the kinematically selected population, inflating the mass. There is a potential for a smaller effect from this contamination on the density profile parameters, since the two populations are sourced from the same density profile, but the spatially complicated nature of the kinematic effective selection function (Figs 5 and 6) may imprint slight spatial irregularities on this contamination. Indeed, we do see some slight mismatches between input parameters and inferred parameters for the kinematically selected populations, but they are minor. The mass on the other hand should be purely overestimated by roughly a factor of the inverse purity, which is what we see: 10^8.1 ≈ 10⁸/0.8. We address the implications of this effect on the fittings to real data later in the discussion when we assess systematic uncertainties.

To ensure that the disc contamination model works as expected we independently analyse the same mock, but without kinematics (i.e. no kinematic cuts on the fitting sample and no kinematic effective selection function). To the mock we add stars sourced from an identical stellar isochrone and initial mass function but tracing the exponential disc density profile used in the contamination model outlined in Section 3.3, which has radial scale h_R = 2.2 kpc and vertical scale h_z = 0.8 kpc. We add stars from this disc model sufficient to set f_disc = 0.4 (76 stars for this mock). The process to determine the correct number of stars to add is nearly identical to the process of normalizing the density profile described in Section 3.5, requiring calculation of the effective volume for both the halo and disc density profile, from which the relevant density at the position of Sun can be determined. We then fit the appropriate model to these mock data, and show the results also in Fig. A1. Again, we confirm that all model parameters, including f_disc are recovered within the confidence intervals. Note that here we did not consider contamination from a disc component in the context of our kinematically selected mock subsamples, since we do not have a good kinematic representation of these stars (see LBM22, and discussion above).

4 RESULTS

4.1 Fits to GS/E subsamples

We first perform the fitting procedure introduced above, including the kinematic effective selection functions, to each of our kinematically defined GS/E samples, using each of the density profiles introduced in Section 3.3. We consider each profile with and without disc contamination. The best-fitting parameters and derived masses are shown for each density profile and each kinematic selection in Table 1. In order to gauge how well each model fits the data, we also determine the maximum likelihood by choosing the set of parameters with the highest likelihood from the MCMC chain and optimizing our likelihood function from that set of parameters, again using Powell’s method (Powell 1964). The optimized set of parameters is never far from the highest likelihood set of parameters from the posterior, and the routine always converges quickly. To compare models with different numbers of parameters, and select the best one among them, we use this optimized maximum likelihood, |$\mathrm{max}(\mathcal {L})$|⁠, to calculate the standard Bayesian information criterion |$\mathrm{BIC} = k\ln (n)-2\mathrm{max}(\mathcal {L})$| (Schwarz 1978) where k is the number of parameters, and n is the number of stars in the specific kinematically selected sample. We also calculate the similar Akaike information criterion |$\mathrm{AIC} = 2k-2\mathrm{max}(\mathcal {L})$| (Akaike 1974). Both quantities penalize fits which use more parameters, but the BIC places stronger emphasis on the number of stars used for fitting, while the AIC places more weight on the value of the likelihood function. Maximum likelihoods, BICs, and AICs are shown in Table 2. The value of each parameter is normalized by the value derived for the single power-law model for the corresponding kinematic selection. Different kinematic selections cannot be compared in the manner we have outlined since the samples are different.

4.1.1 The e − L_z and ad samples

The fits to the e − L_z and ad samples are similar and so we describe them together first. Beginning with the power-law models: for the SPL model the power-law indices are modest, lying in the range 2 < α₁ < 3. BPL models have very shallow inner indices, which tend to be less than ∼1.2. These BPL models have outer indices which are steeper than the single index found in the SPL case, however, with values tending to be greater than 3.5. This difference in index is large, indicating a steep drop-off in density at the break radius, which tends to lie between 15 and 30 kpc. The e − L_z sample is best-fit by smaller break radii, only 15–20 kpc, while the ad sample is best-fit by larger break radii between 20 and 30 kpc. The uncertainties, particularly on the fits to the ad sample, are large enough that the two sets of values are broadly consistent with one another. The exponentially truncated SC models are interesting in that their power-law indices are negative, indicating that density rises with radius. However they have incredibly short truncation radii, less than 10 kpc. So the density only increases slightly in the very inner halo, and then drops off quite quickly with increasing radius as the exponential term overtakes. Thus, the overall density profile looks similar to that obtained for the broken power-law models (Fig. 9, which will be introduced in Section 4.3, shows an example of these sorts of density profiles). Adding a second break to the profile increases the steepness even more, with DBPL models showing outer power-law indices of α₃ > 6, and inner and outer break radii between 15 and 20 kpc, and 30 and 40 kpc, respectively. As in the case for the BPL break radii, when factoring in uncertainties, the sets of break radii for the DBPL models are broadly consistent between the e − L_z and ad samples.

Disc contamination appears substantial at a glance for both samples, with f_disc taking values between 0.3 and 0.6. However, it is important to keep in mind that f_disc expresses the fraction of density at the solar position contributed by the contamination model, not the fraction of contaminating disc stars in each sample. Given that the observable volume of the contaminant disc model is vastly smaller than the observable volume of any halo model for equivalent densities at the solar position (i.e. f_disc = 0.5), even moderate values of f_disc do not indicate substantial contamination. The relationship between f_disc and the fraction of contaminating stars is non-linear, and the number of contaminating stars only becomes substantial once f_disc becomes very close to 1. To this point is the fact that during the construction of mocks in Section 3.6 we only had to add to stars from the disc mock at the level of about 10 per cent of the total (76 mock disc stars added to 645 mock halo stars) to set f_disc to 0.4. We convert f_disc from an expression of density to the number of contaminating stars for our e − L_z and ad models, finding that the range 0.3 < f_disc < 0.6 corresponds roughly to a range of disc star contamination fraction of 0.05–0.1. This is in line with our expectations from the construction of mocks and indicates that disc contamination is not a major concern from a numerical perspective. A final note on this point is the fact that best-fitting density profile parameters and masses are broadly similar between profiles with and without disc contamination.

Shape parameters vary somewhat between e − L_z and ad models, but are broadly consistent across non-SPL fits. In general, we find p of order 0.5–0.7 and q of order 0.4–0.6, which suggests significant triaxiality. For ad models the values of p and q are very similar, while for the e − L_z models p tends to be greater than q. For e − L_z models, the value of η is close to unity, indicating near alignment with the Galactocentric coordinate system. For density profiles such as these, where η is close to 1, the value of θ is not particularly significant (refer to the definition of |$\hat{z}$| in Section 3.3). ad models, on the other hand, tend to have a large spread in η, with median values being around 0.8. The range of values is allowed because p and q are similar, causing degeneracy in rotation about the principal axis. Despite the flexibility in η, the value associated with the maximum likelihood set of parameters is typically close to 1. The final shape parameter, the azimuthal rotation ϕ, is consistent between e − L_z and ad samples between 90° and 100°.

Regarding mass, the e − L_z sample fits are found to have typically lower masses than those to the ad sample, ranging between 10^7.9–10^8.2 (∼0.8–1.6 × 10⁸) M_⊙ and 10^8.1–10^8.5 (∼1.3–3 × 10⁸) M_⊙, respectively. The wide range in masses for each set of samples is mostly driven by the difference between SPL and non-SPL density profiles. SPL density profiles will carry substantial mass at large radii due to the large integrable volume, while the other profiles are all broken or truncated in a way that substantially decreases the density at large radii. Discounting SPL profiles, the range of masses for the e − L_z and ad sample fits are much narrower, being about 10^7.9–10⁸ (∼0.8–1 × 10⁸) M_⊙ and 10^8.1–10^8.3 (∼1.3–2 × 10⁸) M_⊙, respectively. The mass is not consistent between these two samples, however, differing by nearly a factor of two, although the uncertainties, of order a quarter of a decade, are large enough that the mass ranges can be reconciled with one another.

4.1.2 The |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample

While the best-fitting parameters are in reasonable agreement for the e − L_z and ad samples, the results of the fits to the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample tell a different story. Here, the power-law indices tend to be somewhat steeper, and the fits using the exponentially truncated SC models do not return the type of fit described above where the power-law index is negative and the truncation radius short. Instead the inner power-law index is very reasonable (α₁ ≈ 1.5–2.5) and the truncation radius is at about 35 kpc. BPL and DBPL fits do steepen, however, at break radii which are not dissimilar from those found for e − L_z and ad models, but the difference between the indices is not as stark. This difference in power-law index and break radii are reflected in the derived masses, which are substantially larger for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample than the e − L_z or ad samples. In contrast to e − L_z and ad models, now q is greater than p with values ranging from 0.65 to 0.75 and 0.5 to 0.65 respectively. η, θ, and ϕ are similar to the values derived for the ad and e − L_z models, with η being near unity, indicating only modest tilting of the ellipsoid out of the Galactocentric X-Y plane. The angle ϕ is between 100° and 110°, somewhat larger than the e − L_z and ad values, but still orienting the major axis of the ellipsoid towards the Galactocentric Y axis. Finally, disc contamination is relatively constant at about f_disc ∼ 0.65, which corresponds to a contamination fraction of ∼0.08.

The masses of fits to the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|samples are substantially larger than those for the e − L_z and ad fits, roughly 10^8.6–10^8.8 (∼4–6 × 10⁸) M_⊙. This is largely due to the fact that the density profiles tend to be unbroken, and therefore contribute substantial mass at both large radii where the density otherwise drops off, and small radii where the density profile otherwise flattens. The masses for SPL fits are larger than broken fits, but the difference is not as marked as in the case for e − L_z and ad fits, again since the density profiles tend to be relatively unbroken. These masses obviously cannot be reconciled with the findings for the e − L_z and ad samples, even when factoring in uncertainties. We will discuss later in this section as well as the next why we favour the results for the e − L_z and ad samples over these obtained for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample, and also provide explanation for why the fits to the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample resemble unbroken power laws.

4.1.3 Choice of the best-fitting density models

We select the best-fitting model for each kinematically selected sample guided by the values of |$\mathrm{max}(\mathcal {L})$|⁠, AIC, and BIC from Table 2. Larger values of |$\mathrm{max}(\mathcal {L})$| indicate a better fit, while smaller values of AIC or BIC do the same. In general, we see that maximum likelihood increases or is constant as the number of model parameters increases. The values of AIC and BIC, however, reveal that increasing the parametrization beyond SC models to BPL and DBPL models is not statistically justified. For the e − L_z and ad samples the values of AIC and BIC peak for the SC+D and SC model respectively. For |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| AIC and BIC peak for SPL+D model. In each of these cases the maximum likelihood does not increase substantially beyond the chosen model as model complexity increases. We highlight these best-fitting models in grey in both Tables 1 and 2. We also show the posteriors for these best-fitting models in Fig. 7.

We do note that even though these goodness-of-fit criteria favour simpler models over more complex models, that the more complex models tend to match the simpler models in the case of each sample. For instance, the SC and BPL models are very similar in regards to shape parameters and derived masses in the case of the e − L_z and ad samples. While the power law and break radii appear at a glance dissimilar, as mentioned above and as we shall show, these profiles are not particularly different. With this in mind, even though the simpler models are favoured from a statistical perspective, we will discuss the results of the BPL and DBPL models, particularly the break radii and different power laws.

4.2 Fits to the whole halo

We also performed fits using the same set of density profiles to the whole halo sample, without using any kinematic selection functions but otherwise following the exact procedure laid down in Section 3. We also calculate |$\mathrm{max}(\mathcal {L})$|⁠, AIC, and BIC with the results. The best-fitting parameters and likelihood-based parameters are shown at the bottom of Tables 1 and 2, respectively. Best-fitting structural parameters are broadly homogeneous for each of the density profiles. The inner power-law index, α₁, is steeper than for any of the kinematically selected samples, tending to lie in the range 2.5–2.9. Break radii are somewhat trivial, with the posteriors hugging the upper domain limit of 55 kpc. Shape parameters are not dissimilar from the e − L_z and ad fits, although the halo is much closer to a simple oblate ellipsoid, with p between 0.8 and 0.9 and q taking values between 0.5 and 0.6. Orientation parameters suggest strict alignment with the Galactocentric frame, with η very close to 1, which renders θ uninformative. ϕ takes values of about 100°, consistent with the values found for the kinematic samples, suggesting alignment of the major axis of the density ellipsoid with the Galactocentric Y-axis. Finally, f_disc is higher than for the kinematic subsample fits, with best-fits occupying a narrow range 0.73–0.75, which are equivalent to a numerical contamination fraction of 0.18–0.2.

Masses for the whole halo are much larger than for the GS/E subsamples, lying between 10^8.8 and 10^8.9 (∼6.5–8.5 × 10⁸) M_⊙ depending on the density profile. Since the break radii lie at large Galactocentric distances, the fractional difference between SPL and non-SPL models is much less substantial than for the other kinematic subsamples. Taking the same approach to gauging goodness-of-fit as above, we can say that the best-fitting density profile is the SPL+D model.

One interesting trend is that fits to the halo sample are very similar to those for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample. Power-law indices of the best-fitting models are both about 2.5, and the masses are comparable. The shape parameters are dissimilar between the two best-fits, But the orientation parameters are very similar. Both f_disc and the mass of the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| best-fitting profile is more similar to those of the halo best-fitting profile than either the e − L_z or ad best-fits. All this suggests that perhaps the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| selection is not as high in purity as suggested by LBM22, and the resulting sample may be closer to an intrinsic whole-halo mix of high and low β populations, therefore yielding similar fits.

4.3 Interpretation of results

Fig. 8 shows histograms of the distance modulus distribution for each of the samples considered here: the halo sample as well as the three kinematically selected GS/E samples. Overlaid on each of the histograms are characteristic distributions of distance modulus from the posteriors of the best-fitting density profiles, calculated by determining the effective volume as a function of distance modulus (marginalized over all fields). These characteristic distributions are derived by taking the median of one hundred different samples from the posterior parameter distribution, selected randomly. While calculating the median of the posterior parameter samples, we also record the 16th and 84th percentile of the distribution at each distance modulus. We then take an average of these percentiles across each of the eight density profiles (SPL, SC, BPL, DBPL with and without disc contamination), which produces the characteristic uncertainty in the model distributions shown in grey at the top of each panel. Each of the density profiles are clearly able to provide reasonable fits to each of the stellar samples, with posteriors lying well within the counting uncertainty (⁠|$\sqrt{N}$| error on the counts in each distance modulus bin, shown as grey error bars) of the distance modulus distributions when also factoring in model uncertainty. There are small variations between density profiles for each sample, and particularly for density models with (dashed lines) and without (solid lines) disc contamination.

Fig. 9 shows the mass density as a function of triaxial radius m (top panel) and Galactocentric X (bottom panel) for the best-fitting models for each of the samples. The solid lines and the surrounding fill show the median and central 68th percentile confidence interval for the density profile generated using one hundred different sets of parameters, drawn randomly from the posterior distributions. We see emphasized here the similarities between the fits to e − L_z and ad samples, as well as those between the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| and whole-halo samples outlined above. The e − L_z and ad fits are extremely flat in the inner Galaxy, even rising in the case of the e − L_z fit, and drop-off precipitously at about log₁₀(m/kpc) ∼ 1.2–1.5 = 15–30 kpc. The |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample best-fit is very similar to the whole-halo fit, a simple power law with nearly identical indices, but offset to lower density.

The dashed grey lines in Fig. 9 shows density profiles which are half and twice as dense as the median profile for the halo sample. We show these lines because one important consistency check of our results is the observed fact that GS/E constitutes roughly 50 per cent of stars in the solar vicinity (Belokurov et al. 2018; Fattahi et al. 2019; Lancaster et al. 2019). Our GS/E density profile fits do not quite rise to this level, but are not far off, sitting between one quarter and one half of the density of the best-fitting halo model within uncertainties. Interestingly, at slightly larger radii between log₁₀(m) ∼ 1–1.5 or m ∼ 10–30 kpc the e − L_z and ad fits do rise to – and exceed for a small range of radii – this level. It is useful here to compare with the density profiles as a function of Galactocentric X, which both factor in the effect of p and q and represent the density in a manner more focused on the solar neighbourhood. This panel reveals that indeed the e − L_z and ad density profiles do peak near the solar vicinity and not beyond, but do not dominate the density at any distance.

Fig.10 shows density contours for cross-sectional views of the best-fitting ad profile. We show only one set of density contours, since the approximate shape of the best-fitting profiles is similar among each of the samples (i.e. triaxial, slightly tilted, and with the major axis oriented near the Galactocentric Y-axis). We can now see the density profile described by the best-fitting parameters outlined in in the previous sections. The principal axis of the ellipsoid is rotated towards the Galactocentric Y-axis (ϕ ∼ 90°), and is tilted slightly out of the Galactic X-Y plane (η near 1). Note that profiles where η is equal to 1 will simply have the principal axis in the X-Y plane, but oriented still towards the Y-axis. The inclination of the principal axis with respect to the Galactocentric X-Y plane for the ad best-fitting profile shown here is −16.4° (The central 68 percentile interval for this value is [ −35.5°, 6.6°]).

Moving to the mass, in general we find evidence for a light GS/E remnant, roughly 0.8–1.6 × 10⁸ M_⊙ in two out of our three kinematic samples under consideration: those created with the e − L_z and ad selections. The |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample fits suggests the opposite, that GS/E is heavy, roughly 5 × 10⁸ M_⊙. We posit that this is an issue of contamination, arising from the lower purity of the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| selection. Indeed, LBM22 find that the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| is the least pure among the kinematic spaces that they studied which we also use here. To further explore this angle, we determine the purity of each of our selections in a manner similar to how we calculated the kinematic effective selection function in Section 3.2. The purity here is defined as the fraction of samples within each kinematic selection from the high-β component with respect to the total number of samples in the selection, from both high- and low-β components. We marginalize along each line of sight and over each field, weighted by the effective selection function (which generally describes the number of stars expected at a given location, modulo the density). We find the characteristic purity for the three spaces is as follows: e − L_z: 0.84, ad: 0.90, and |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠: 0.78. This result is similar to the purity quoted by LBM22 for each space, which was derived using similar models, albeit in a different manner.

A lower purity would impact the results to the fits in two ways. First, more stars would be counted in the fitting sample, which would drive up the normalization factor for the density profiles in equation (14), which in turns drives up the calculated masses (see discussion of this effect in Section 3.6 in the context of our mock data). A second, more subtle, effect is that we would expect the fit to tend more closely to the fit to the whole halo sample. The best-fitting density profile for the halo sample is steep in the inner part of the Galaxy, and lacks a cut-off at 20 to 40 kpc such as the e − L_z and ad best-fitting density profiles have. Both of these effects contribute to a heavier overall mass, and if the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| profile were to trend towards that sort of fit due to lower purity in the underlying sample then we would expect the mass to come out large for that reason as well. Although, it may be argued that the difference in purity (of order 6 to 12 per cent) between the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| and e − L_z/ad spaces could not possibly account for this difference. We also acknowledge that the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| space is ‘absolute’ in the sense that the dependant variable is the radial action. This is contrasted to the ad and e − L_z spaces which are much more regular, the ad space being explicitly normalized by the total action, while the e − L_z space being normalized in the sense that eccentricity may only vary between 0 and 1, and our selection includes eccentricity of 1, which orbits approach as they become increasingly radial. The kinematic effective selection functions are derived assuming a potential, specifically MWPotential2014 of Bovy (2015), but it is applied to real data. Now the actions and conserved quantities of the real data are also derived assuming the MWPotential2014, however their underlying kinematics are dictated by the true potential of the Galaxy, whereas the kinematics of the DF samples used to compute the kinematic effective selection function are dictated by MWPotential2014. This mismatch could cause problems for spaces such as |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠, which will be much more sensitive to the underlying potential than the ad or e − L_z spaces.

Interestingly, the works of Feuillet et al. (2020) and (Carrillo et al. 2023) suggest the opposite of the conclusion at which we have arrived here, and they posit that |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| is the superior space for selecting GS/E debris. Feuillet et al. (2020) conclude this by studying the metallicity distribution functions of kinematically selected halo stars, although they do not consider eccentricity or a space similar to the action diamond, finding that |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| is better than energy versus L_z, the Toomre diagram, and radial versus tangential velocities (this is actually consistent with the conclusions of LBM22). So while they determine that |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| is better than energies and velocities we should not take that to mean that eccentricities and the action diamond may not be better than |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠. Carrillo et al. (2023) come to their conclusions while studying simulated GS/E remnants in a cosmological context. They consider all of the spaces studied by LBM22 but use different specific kinematic selections in some instances. For example they consider eccentricities greater than 0.7 (our eccentricity selection which is very close to 1), and a slightly larger selection in the action diamond. Their findings regarding eccentricity are broadly consistent with LBM22, who also find that e > 0.7 yields poor purity (again, note the eccentricity selection used in this work is much more restrictive). With respect to the action diamond, their results are tougher to reconcile with those of LBM22 since the selections are similar, although their selection is again broader. In general, the work of Carrillo et al. (2023) benefits from knowledge of the ground truth and the realistic dynamics innate to simulations, however they lack the necessary abundance information in their simulations with which to make cuts which are typical in observations. So while perhaps their work does not discount our conclusions about the purity of kinematic spaces it certainly highlights the need to take care when kinematically selecting GS/E debris. Finally, the work of Donlon & Newberg (2023) suggests that |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| may be contaminated at high |$\sqrt{J_\mathrm{R}}$| at the level of ∼1/3 or greater by non-GS/E accretion remnants, which is in rough agreement with the findings of LBM22.

4.4 Altering the kinematic effective selection

When constructing the kinematic effective selection function, we noted the inconsistency between the fact that we had to assert a density profile for the stellar halo when defining the DFs, yet we are employing the DFs to determine the density profile of the stellar halo. To investigate the impact of this assumption, we compute new high and low-β DFs based on the density profiles that we have determined for the GS/E and whole-halo samples. For the low-β component we construct a density profile inspired by the best-fitting SPL model to the whole-halo sample. We use a single power-law density profile with index α = 2.5, normalized such that it has a mass of 6 × 10⁸ M_⊙ between 2 and 70 kpc (i.e. roughly the mass of our whole-halo fit less our mass derived for the GS/E remnant). For the high-β component we take our cue for the density profile from the BPL fits, which for the ad and e − L_z samples have approximate inner and outer power laws of 0.8–1.2 and 3.8–5.2, and break radii between 17 and 27 kpc. We therefore conveniently approximate the high-β component using a Hernquist density profile (Hernquist 1990), which has an inner power-law index of 1 and an outer index of 4, with scale radius 20 kpc and a mass between 2 and 70 kpc of 1.5 × 10⁸ M_⊙. We choose to model the high-β density profile in this way for ease of computation, and the BPL models which inspire it provide fits to the data with equivalent likelihood compared with the SC model which is our claimed preference. We keep the values of β fixed at 0.3 and 0.8.

Samples are drawn and kinematic quantities are computed as outlined in Section 3.2. Now, however, we compute purity by counting the number of high-β samples in each selection with respect to the total number of samples, but we weight the contributions from the low and high-β components by the local density of the respective density profile (note the completeness of the high-β component is independent of the local densities of the two components). We replicate Figs 5 and 6, which show the kinematic effective selection functions, as Figs B1 and B2, which are shown in Appendix  B. Examining the new kinematic effective selection function (which we denote at |$\mathfrak {S}_{2}^{\prime }$| in these figures) reveals remarkable consistency between them and the original kinematic effective selection functions shown in Figs 5 and 6. On Fig. B1 we overlay the median trends, calculated across all fields at each distance modulus, of each of the kinematic effective selection functions to summarize the differences. The overall typical value of |$\mathfrak {S}_{2}$| decreases slightly for each selection, from ∼0.5 to ∼0.45 for the e − L_z space for example, but the overall form remains very similar. This implies that we could expect very similar fits to the data using these new kinematic effective selection functions, except that the derived mass should increase by a small amount, since the typical value of |$\mathfrak {S}_{2}$| is somewhat smaller. We compute the ratio between these two median trends, and calculate an average over the range of distance modulus weighting by the distribution of distance moduli in the sample (top panel in Fig. B1 which is the same as the top panel of Fig. 5). We find that the resulting typical fractional change in the kinematic effective selection function is −0.09 (e − L_z), −0.07 (ad), and −0.13 (⁠|$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$|⁠). We would therefore expect the masses to increase by a small amount, roughly 10 per cent, but overall this effect is negligible. In general, these results suggest that the DF-based corrections we derive and apply in our modelling framework are resilient to modest changes, and therefore that as long as the assumptions underlying the models used to generate the corrections are realistic that their specific form is generally unimportant.

The most interesting aspect of this analysis is actually in examining purity. Above we noted that the purity is lower for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| selection than the e − L_z or ad selections, but not by a significant amount (0.78 versus 0.84 and 0.9 respectively). We recalculate the typical purity using the new DFs here, again by summing the purity over all fields and distances weighted by the value of the effective selection function. We find that the new purity values are 0.78 for the ad selection, 0.65 for the e − L_z selection, and 0.52 for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| selection, so all values have been lowered. Note that the overall lowering of the purities here will be largely driven by the fact that the low-β density profile is quadruple the density at the solar position compared with the high-β component. This should not be taken as indication that these selections are now inferior, and they are likely still optimal (i.e. were one to repeat the analysis of LBM22). The differences in purity among the kinematic spaces are even more stark now than in the case of the original DF models, and suggest that |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| is truly an inferior kinematic space for which to select GS/E samples. Fractionally, the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| space experiences the largest reduction in purity, decreasing 33 per cent. The fractional decrease for the e − L_z and ad spaces are smaller in magnitude, 22 and 14 per cent respectively. So not only does the ad space maintain the highest purity, but its purity is changed by the smallest amount when pivoting from the original DF models to those based on our results. These findings suggest that the ad space, in particular, should be preferred given its high purity and resilience of its purity to changes in DF models.

4.5 Our preferred results for the density profile of GS/E

We close out this section by summarizing our main findings for the density profile and mass of GS/E in light of our primary results and secondary analyses. We have already noted that there appears to be broad consensus between the fits to the e − L_z and ad samples, a consensus which differs significantly from the results of the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample. We also noted that the results for the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample are similar in many ways to the results for the fits to the whole halo. When factoring in the assessment of purity, for both the original kinematic effective selection functions and the second set of kinematic effective selection functions derived based on our initial results, we arrive at a neat conclusion. The e − L_z and ad results more faithfully trace the underlying reality of the GS/E remnant, while the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| results are biased by excessive contamination from the low-β halo such that the density profile parameters are similar to those for the whole halo and the mass is overestimated. For this reason we defer to the former results, and specifically those for the ad sample since it has the highest purity, when crafting our conclusions regarding the density profile and mass of GS/E.

With this in mind we can say that the best-fitting density profile for the GS/E remnant is the SC model, a single power law with exponential truncation, with |$\alpha _{1} = -0.57^{+1.58}_{1.95}$| and |$r_{1} = 8.57^{+12.58}_{-4.36}$| kpc. This type of density profile, which gently rises in the inner Galaxy and drops rapidly at modest Galactocentric radii, is non-traditional. For this reason we also report the power-law indices and break radius of our BPL fit (no disc contamination) which has inner power law |$\alpha _{1} = 1.05^{+0.86}_{-0.72}$|⁠, outer power law |$\alpha _{2} = 3.79^{+2.44}_{-0.98}$|⁠, and break radius |$r_{1} = 23.59^{+15.55}_{-9.30}$| kpc. These two profiles are, in practice, very similar over the range of radii probed and are equivalent in terms of maximum likelihood, yet the SC model is favoured by BIC/AIC since it has one fewer parameter. We find that GS/E is moderately triaxial, with the ratios of the Y–X and Z–X axes in the rotated (not Galactocentric) frame being respectively |$p = 0.54^{+0.22}_{-0.2}$| and |$q = 0.46^{+0.27}_{-0.2}$|⁠. The principal axis is rotated towards the direction of Galactic rotation with |$\phi = 99.5^{+13.9}_{-14.9}$| degrees, and is tilted slightly out of the Galactic plane with |$\eta =0.84^{+0.12}_{-0.22}$| and |$\theta = 147.0^{+43.5}_{-95.3}$| degrees, which equates to an inclination of the principal axis of |$-16.4^{+23.0}_{-19.1}$| degrees. Finally, the mass is |$1.45^{+0.92}_{-0.51}\times 10^{8}$| M_⊙.

We also calculate a few supplementary parameters that derive from our best-fitting density profile and associated mass. First, we consider the triaxiality parameter, T, commonly given as

where b = max([p, q]) and c = min([p, q]) in the context of our model parameters, and a = 1. We calculate this parameter for each of the samples in the posterior for the SC fit to the ad sample, arriving at a value of |$T=0.79^{+0.13}_{-0.35}$|⁠. This indicates that the GS/E density profile is preferrably prolate, but only to a modest degree. Second, we take our stellar mass and determine the corresponding dark matter halo mass using a stellar mass–halo mass (M₂₀₀) relation. According to the relation derived by Behroozi, Wechsler & Conroy (2013) (and validated by Read et al. 2017 for stellar masses ∼10⁸ M_⊙) a stellar mass of 1.5 × 10⁸ M_⊙ implies a halo mass of about 5 × 10¹⁰ M_⊙ at redshift z = 0, with halo masses larger by a few tenths of a decade expected for the same stellar mass at redshift z = 2 (the approximate GS/E merger epoch). The Milky Way-oriented stellar mass–halo mass relation of Nadler et al. (2020) predicts a comparable halo mass of 6 × 10¹⁰ M_⊙. Given that the Milky Way at the time of the merger had a halo mass that is approximately half of that today (so ≈5 × 10¹¹ M_⊙; e.g. Mackereth et al. 2018), this means that the GS/E merger was approximately 1:8 to 1:10, constituting a minor merger.

5 DISCUSSION

5.1 The distribution function models and the kinematic effective selection function

Our results presented here are perched upon the DF-based models introduced by LBM22. While we have gone to great lengths in this paper to demonstrate that the corrections we derive for kinematic selections allow recovery of the properties of the underlying population, we must address the limits of these types of models. First, the functional form of our DFs is a simple, but flexible, ansatz (7) which is designed to incorporate many of the basic features we expect from stellar halo populations: simple radial dependency and a variable anisotropy which is assumed to be radially constant. It lacks, however, many of the ingredients that we know stellar halo populations to possess. For example, we know that GS/E is likely to have modest net rotation (Helmi et al. 2018; Lancaster et al. 2019; Iorio & Belokurov 2021), and to have asymmetric or complicated phase-space features when comparing the debris at high and low L_z (Helmi et al. 2018; Naidu et al. 2020), or at large Galactocentric radii (Chandra et al. 2022). Additionally, as we have shown here and as has been shown in nearly every other work that has measured the density profile of GS/E, the remnant is triaxial (Iorio & Belokurov 2019; Han et al. 2022). None of these features are included in our DF models (but of course our density models are triaxial), which may result in discrepancies between our corrected results and the true GS/E population which are challenging to predict. Finally, the mathematical form of the DF implies divergence when β > 0 as L → 0, which may compromise its suitability for modelling purposes (as discussed by Binney 2014). These effects may be particularly exacerbated when modelling the highly anisotropic GS/E remnant.

Despite these admitted deficits, our models do an excellent job of reproducing the basic assumptions and features of GS/E, namely that it is a stellar population with extreme radial anisotropy and near-zero L_z (see fig. 4 in LBM22, for example). It is likely that our models are closer to a true description of the phase-space density of GS/E than not, and therefore that our results should be taken at face value. LBM22 tested their models against changes in underlying potential, and they also tested the resiliency of their selections to changes in radius, finding that neither change produced a noticeable effect on the interpretation of their results. Those findings may be carried over to our work to suggest that our results would be resilient to changes in underlying potential as well, while variations the quality of the selection with radial position are explicitly handled by our spatially dependent kinematic effective selection function. Additionally, we explored variation in the kinematic effective selection function and found that at least one of our kinematic selections, ad, is resilient to iterative improvements in the selection function based on our initial results. These tests all suggest that our preferred findings of a light, triaxial GS/E are accurate and would not change substantially given more sophisticated corrections based on more realistic kinematic models.

Despite our confidence in these results, we do plan follow-up work which will both help to reinforce our findings, but also to expand these sorts of techniques for use with other remnants, which tend to have more complicated kinematic properties, requiring more complicated kinematic selections and corrections during the modelling procedure. A tangential, yet perhaps more realistic, approach would be to model the stellar halo as a mixture of self-consistent DFs (see, e.g. Lancaster et al. 2019; Iorio & Belokurov 2021) for similar studies using Gaussian mixtures), however the large number of parameters and the intractability of computing such DFs on the fly mean that it would be somewhat challenging. On a similar note, exploring the use of more realistic DFs should be of top priority. For example action-based DFs (Binney 2014; Posti et al. 2015; Williams & Evans 2015) including flattening and rotation (Binney 2014). The construction of self-consistent triaxial DFs is an area of ongoing effort (Sanders & Binney 2015; Binney 2018), yet these would obviously prove ideal given the observed geometry of the GS/E remnant.

Another fruitful avenue for future work would to test these sorts of modelling techniques on merger remnants in simulated Milky Way analogues sourced from cosmological simulations. The EAGLE (Schaye et al. 2015), AURIGA (Grand et al. 2017), and FIRE (Wetzel et al. 2016) sets of simulations, among others, all have rich veins of recent work on GS/E-like remnants (e.g. Fattahi et al. 2019; Mackereth et al. 2019a; Orkney et al. 2023; Horta et al. 2023b). In particular, examining the phase space distributions of debris in detail would be useful. There have been recent works that have looked at the phase space distribution of simulated debris (e.g. Amorisco 2017; Jean-Baptiste et al. 2017; Naidu et al. 2021; Amarante et al. 2022; Carrillo et al. 2023; Orkney et al. 2023), and we advocate a more in-depth look specifically at how well various families of DFs can be fit to merger debris. This type of study could place necessary constraints or priors on the vast volume of parameter space which would undoubtedly accompany a multi-DF assessment of the Milky Way stellar halo.

5.2 The measurement of the mass of the GS/E remnant

The core measurement made in this paper is the stellar mass of GS/E, which we estimate to be about 1.5 × 10⁸ M_⊙ in our best-fitting models. We will discuss our results in the context of the literature in the next section, but so as to not bury the lead: our derived mass for GS/E is much lower than other findings. It therefore behooves us to assess our measurement of the mass of GS/E and investigate possible sources of systematic error.

First, we examine the impact of a few of our sample-selection and model-fitting decisions on our final results. We repeat our fits without the log g uncertainty and relative distance uncertainty cuts presented in Section 2.1. The Halo sample and the three GS/E samples increase in size by between 2 and 8 per cent when these cuts are removed, but the changes to the results – the values of the best-fitting parameters, the derived masses, and the overall narrative presented in Section 4 – are negligible. According to the BIC/AIC the best-fitting density profiles are the same for the Halo sample and each of the three GS/E samples, and the masses change by only a few per cent (the fit to the ad sample increases by 0.01 dex), far smaller than the statistical uncertainty on each measurement. We also check the assumptions which go into the construction of the effective selection function and the density profile normalization factors, namely the choice of [Fe/H] and age ranges. While Montalbán et al. (2020) show that typical GS/E stars ages are about 10 Gyr old, the distribution extends down to ∼7 Gyr old. Additionally, with regards to our choice of [Fe/H] range, whereas we elected to use a minimum [Fe/H] of −3 for the halo sample, the data only barely extend lower than [Fe/H] of −2. We replicate our fits but using a new isochrone grid which extends from 7–14 Gyr old (therefore rebuilding our effective selection function) and using a minimum [Fe/H] of −2 for the halo population (same as for the GS/E samples), but holding all other aspects of our analysis constant. We again find that this has a negligible impact on our results or overall narrative. Masses decrease by between 0.02 and 0.07 dex (the largest decrease is for the ad GS/E sample), best-fitting density profile parameters change little, and the AIC/BIC prefer the same density profiles. We therefore assign a systematic uncertainty which lowers the mass of GS/E by 0.07 dex (15 per cent), or about 0.22 × 10⁸ M_⊙, due to choices in parameters for sample selection and effective selection function grid construction.

Second, as discussed when validating our method in Section 3.6, we expect the mass of our fits to be slightly overestimated due to contamination from the low-β component. Now in the case of the mock data, both the low and high-β components were drawn from the same density profile, and therefore one would only expect bias due to the excess in number counts (modulo some higher order effects from the spatial dependence of the selection). Secondarily, and as discussed above in the context of the disparity between e − L_z and ad fits versus |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| and whole-halo fits, when the underlying density profiles differ then contamination from the low-β component can lead to a best-fitting density profile which more closely resembles that of the low-β component than the target density profile of the high-β component. Given that we assume the low-β component has a ‘more massive’ density profile (i.e. the whole halo has a steeper power law in the inner Galaxy and a less marked break at intermediate Galactocentric radii) then we might also expect that contamination leads to an overestimation of the mass.

The first of these effects is simple to gauge: we expect to overestimate the mass by a factor equal to the inverse of the purity of the selection (this is neatly observed when determining the mass of our mocks). In the case of the ad selection the purity is estimated to be 0.9 for our initial kinematic effective selection function, and drops to 0.78 for our more realistic kinematic effective selection function created in Section 4.4. The purity found by LBM22 for the ad selection was 0.86. The second effect, the biasing of the shape of the density profile itself is challenging to gauge. We have posited above that contamination from the low-β halo causes the increased mass found with the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| fit, however this is higher than that for the ad fits by a factor of 3–4. Without any good way to gauge the impact of the contamination on the shape of the GS/E halo density profile we therefore fall back on the simple, conservative estimate being of order the inverse purity. We therefore adopt as our first systematic uncertainty one which lowers the mass to 0.78 (the lowest purity estimate for the ad selection) times the best-fitting mass, which is about 0.3 × 10⁸ M_⊙. We do note though that we have attempted to mitigate this source of systematic uncertainty by cutting stars with high [Al/Fe] from our kinematically selected GS/E samples (see Fig. 2), since these stars are more likely to be part of the in-situ, low-β stellar halo.

Building on this, while we believe the trimming of stars with high [Al/Fe] is a well-motivated way to remove low-β and thick disc contaminants from our kinematic selections, it is possible that these are genuine GS/E stars which should be included in our calculations. For the ad selection there are 4 stars which are removed on this basis, compared with 73 stars in the sample (there are similar fractions removed for the e − L_z and |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| spaces, of order 5–8 per cent). Assuming that inclusion of these stars would not change the density profile, only increase the normalization factor, this would result in an increase in the mass by 0.08 × 10⁸ M_⊙.

Finally, we consider that in Section 4.4 we investigated iterative improvements to the kinematic effective selection function, finding that more realistic kinematic effective selection functions constructed based on our fits generally have smaller values for a fixed kinematic selection. For example, referring to Fig. B2 we see a typical value of 0.5 for the value of |$\mathfrak {S}_{2}$| for the ad selection, while a typical value for the more realistic |$\mathfrak {S}_{2}^{\prime }$| is more like 0.45. In Section 4.4 we specifically calculated a fractional decrease for the ad space of −0.07. Again, since this is for a fixed kinematic selection, and therefore a constant number of stars in the sample, we expect this decrease in the value of the kinematic effective selection function to manifest as an increase in the derived mass by an approximate factor 0.07, or about 0.1 × 10⁸ M_⊙.

Before summarizing our systematics we discuss three related points. First, we must contend with the fact that we find a much lower fractional density for the high-β GS/E component at the position of the Sun compared with other studies, which tend to attribute about 50 per cent of the density near the Sun to GS/E (Belokurov et al. 2018; Lancaster et al. 2019; Iorio & Belokurov 2021). We find that the density at the position of the Sun for the ad sample is |$\nu _{0} = A\chi = 5.2^{+0.7}_{-0.6}\times 10^{3}$| M_⊙kpc⁻³, compared with a density of |$2.7^{+0.2}_{-0.2}\times 10^{4}$| M_⊙kpc⁻³ for the halo sample. This is a difference by about a factor of 5, about 2.5 less than we would expect for an even mixture of high- and low-β stars at the solar position. Since in our analysis mass scales directly with the density at the solar position, if we demand that the density of the GS/E component rises to the level of 50 per cent of the whole-halo sample, then we would find the mass increases to about 3.8 × 10⁸ M_⊙. Now, we do not consider this to be a genuine source of systematic uncertainty, unlike those noted above, since we have no reason based on a sober assessment of our analysis to think that we have substantially underestimated the density at the solar position. For this reason we simply note that this is what our mass would become in light of these constraints from prior studies.

A second point is that we elected to model the high-β portion of the stellar halo using an anisotropy of β = 0.8 to avoid minor numerical instabilities noted by LBM22 when using DFs with β = 0.9. This is slightly lower than the value of β ≈ 0.9 which is preferred by recent literature (Belokurov et al. 2018; Fattahi et al. 2019; Lancaster et al. 2019; Iorio & Belokurov 2021). Our choice of β = 0.8 actually implies that our derived mass is an upper limit, assuming that the geometry and form of the best-fitting GS/E density profile would be similar if we fit with a kinematic effective selection function derived using β = 0.9. This is because the kinematic effective selection function should typically approach 1 as β increases, which in turn decreases the derived mass. Intuitively, the selections from LBM22 contain within them a larger fraction of a β = 0.9 population than a β = 0.8 population since they probe regions of kinematic space where extremely radial orbits lie, and therefore the required correction for stars missed by the restrictive kinematic cuts is smaller and the inferred mass is smaller. To check this explicitly, we compute a kinematic effective selection function assuming the same parameters as presented in Section 3.2 except that we increase the high-β anisotropy to 0.9. We find that the form of the correction is very similar to the β = 0.8 correction both as a function of distance along each line of sight, as well as on-sky position, except that values are as expected systematically higher. For the ad space, we find typical values increase from ∼0.4 to ∼0.6, and similarly for the e − L_z space we find typical values increase from ∼0.5 to ∼0.7. Assuming that were we to fit our samples with these kinematic effective selection functions that resulting density profiles would be similar to our current findings, we can say that the derived mass would be about 60–70 per cent our current derived mass, or about 1 × 10⁸ M_⊙.

A third and final auxiliary point is to mention that while we are confident in our thick disc outlier model, it is possible that it does not fully capture all possible contaminants. It is well known that the thick disc has a range of geometries depending on the age and abundance of the specific subpopulation (Bovy et al. 2012; Mackereth et al. 2019b), and the oldest, hottest, thickest parts may not be well-described by our simple single-exponential model. This is less likely to impact the measurement of the mass of GS/E, since these contaminants likely have higher [α/Fe] and [Al/Fe] such that they are excluded from our chemical cuts on the GS/E subsamples. The whole-halo sample, on the other hand, may be plagued by such contaminants since it includes higher [Al/Fe] stars. It may be the case then, that our measurements of the total mass of the stellar halo are impacted, and by extension our estimates of the mass ratio of the GS/E merger. However since our stellar halo mass measurements are less than 10⁹ M_⊙ and therefore on the low end of modern estimates (e.g. Deason et al. 2019; Mackereth & Bovy 2020) it is likely that the degree of contamination is minimal.

Now by combining the sources of systematic uncertainty together, and adding uncertainties in quadrature, we arrive at an estimate for the mass including systematics of |$1.45\ ^{+0.92}_{-0.51} \mathrm{(stat.)}\ ^{+0.13}_{-0.37} \mathrm{(sys.)}\ \times 10^{8}$| M_⊙. And to be clear we have not included the overestimations due to the discrepancy between our measured fractional solar densities for GS/E and the literature estimates, or the modification of the mass due to a change in β in this quoted value.

5.3 Comparison of our findings with other work

It is challenging to compare our results with those of other studies due to the unique nature by which we select our sample. Although the spirit of our kinematic and chemical selections is not that different from similar studies (e.g. Mackereth & Bovy 2020; Han et al. 2022), our kinematic selections are far more restrictive and we correct for missing stars using a novel kinematic effective selection function. None the less, we attempt here to provide a context for the density profile parameters and mass of GS/E reported here.

5.3.1 The mass of GS/E

When it was first studied after its discovery, the stellar mass of the GS/E progenitor was thought to be quite high. Helmi et al. (2018) posited the stellar mass was 6 × 10⁸ M_⊙ based on chemical abundance trends, and Belokurov et al. (2018) suggested that the virial mass of the progenitor was >10¹⁰ M_⊙. Early simulation-based studies suggested that the stellar mass of the progenitor could by in the range 10^8.5–10¹⁰ M_⊙ (Fattahi et al. 2019; Mackereth et al. 2019a). In particular, Fattahi et al. (2019) find that within 20 kpc the mass of the progenitor could be between about 2 × 10⁸ and 2 × 10⁹ M_⊙. Das et al. (2020) estimate a total virial mass of about 10¹¹ M_⊙ based on the dispersion of the remnant in action space, which equates to a stellar mass of about 10^9.5 M_⊙ using canonical abundance matching relations. GS/E was also studied using chemistry, with Vincenzo et al. (2019) estimating a mass of 5 × 10⁹M_⊙ and Feuillet et al. (2020) estimating a mass of 7 × 10⁸–7 × 10⁹ M_⊙. Naidu et al. (2020) and Naidu et al. (2021) use H3 data – and N-body simulations in the latter study – to infer a mass in the range of 4–7 × 10⁸ M_⊙. The work of Naidu et al. (2021), where the mass was found to be 5 × 10⁸ M_⊙, in particular is consistent with a number of auxiliary constraints on GS/E, such as the modestly retrograde tendency for the debris (Belokurov et al. 2018; Helmi et al. 2018), the existence of a retrograde tail of debris at large L_z (Helmi et al. 2018; Naidu et al. 2020), as well as the concurrence and alignment with the major axes of the density ellipsoid of the debris with the Hercules-Aquila Cloud and the Virgo Overdensity (first noted by Simion et al. 2019). More recently, Rey et al. (2023) use the VINTERGATEN suite of cosmological zoom simulations to study GS/E-like mergers of varying mass ratios using the ‘genetic modification’ technique. Their results favour a similar mass ratio for the GS/E merger as our finding, but most interestingly they demonstrate that a wide range of merger mass ratios (∼1: 24–1: 2) can lead to very similar halo chemodynamics at z = 0. This may help to explain why our results are at odds with many of the abundance-based measurements which predict a higher mass for GS/E.

The work of MB20 is the most closely analogous to ours, since they use a very similar technique, similar density profiles, and similar sample (they use APOGEE DR14). They find the mass of GS/E to be approximately 3 × 10⁸ M_⊙, although they base this estimate purely on a chemical selection of GS/E inspired by Mackereth et al. (2019a). Another recent work which is very similar to ours is that of Han et al. (2022), who fit triaxial broken power-law density profiles to H3 data, finding a mass of 5.8–7.6 × 10⁸ M_⊙, with the lower mass end of that range corresponding to more restrictive eccentricity cuts on the GS/E sample. Compared to all of these works, we find a much lower mass, with a value of 1.45 × 10⁸ M_⊙ derived for our most favoured kinematic selection and density profile. While it does bring itself into closer alignment with other findings, we do not favour the results of our |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| sample due to contamination issues as discussed in previous sections.

Another set of constraints on the mass of GS/E comes from the study of globular clusters systems associated with the remnant (as tabulated by e.g. Myeong et al. 2018; Massari, Koppelman & Helmi 2019), for which the age-metallicity relation may be linked to the total mass of the host. Kruijssen et al. (2020) use the E-MOSAICS simulations to train a neural network to predict the accretion redshift and mass of the major Milky Way accretion remnants. For GS/E they find a mass of 2.7 × 10⁸ M_⊙ with error bars that are large enough to be consistent with our preferred measurement. Forbes (2020) link the total number of GS/E-attributed globular clusters to its halo mass, which can then be used to infer the stellar mass using canonical relations, which yields an estimate of approximately 8 × 10⁸ M_⊙. Callingham et al. (2022) perform a similar analysis, relying on simulations and kinematics for additional context when grouping GS/E globular clusters, and arrive at a mass for the remnant of about 3.2 × 10⁸ M_⊙. Limberg et al. (2022) carry out a joint analysis of stars and globular clusters associated with GS/E, additionally relying on the hypothesis that the cluster ω Centauri was a nuclear star cluster of GS/E, and estimate a stellar mass of 1.3 × 10⁹ M_⊙ for the progenitor system.

Many lines of evidence for the mass of GS/E rely on the stellar mass-metallicity relation. It is therefore useful to ask whether the mass we derive for GS/E is unphysical given that its mean [Fe/H] remains unchallenged (cf. our Fig. 4 with Myeong et al. 2019; Monty et al. 2020; Hasselquist et al. 2021; Horta et al. 2023a). Naidu et al. (2022) collect the stellar mass and metallicity information for Milky Way stellar halo remnants from the H3 survey. Comparing with their fig. 2, we see that were GS/E to have a mass ∼1.5 × 10⁸ M_⊙ that it would still be consistent with the observed mass-[Fe/H] and mass-[α/Fe] trends. In fact, GS/E would more closely resemble the Helmi streams (Helmi et al. 1999) in both stellar mass and metallicity. Under the assumption that the masses of the progenitors of these two structures are similar, the mismatch in their [α/Fe] abundances (of order 0.1 dex) may require explanation, although there is intrinsic scatter in the relation that may permit the deviation. It could perhaps be the case that the Helmi streams has lower [α/Fe] because it experienced lower star formation efficiency for its mass (e.g. Matsuno et al. 2022). Mackereth et al. (2019a) study the progenitors of high-eccentricity stellar halo debris in Milky Way analogues the EAGLE simulations in order to derive a correspondence between mass and metallicity for comparison with APOGEE. Their observed trends of progenitor mass with [Fe/H] and [α/Fe] (see their fig. 11) would by no means exclude a progenitor with stellar mass ∼1.5 × 10⁸ M_⊙ given a typical [Fe/H] for GS/E of −1.2. Horta et al. (2023b) use the FIRE simulations to investigate the properties of merger remnants in the haloes of Milky Way analogues, and their findings would indeed be challenging to reconcile with our low stellar mass for GS/E. They show that fully phase-mixed remnants with [Fe/H] comparable to GS/E tend to have higher stellar masses. These comparisons suggest that while there may be some tension in placing such a light GS/E in a cosmological context, it would be by no means an outlier when compared with other remnants in the Milky Way stellar halo.

When comparing our results both with all other mass constraints for GS/E, but also specifically those of Han et al. (2022) and Mackereth & Bovy (2020) who also study the GS/E stellar density profile directly, our mass is certainly among the lowest ever claimed. Our findings of a low mass require explanation when compared with these two studies since the techniques are so similar. We posit that the reason for this is that our kinematic selections are far more restrictive. LBM22 compared literature kinematic selections for GS/E with their derived high-purity selections, finding that some commonly used selections for GS/E can achieve purity as low as 0.5–0.6, and many achieve purity in the range of 0.6–0.7. This adds a substantial degree of contamination to samples selected in this way, artificially increasing inferred masses. This effect is compounded, however, by the fact that a contaminated sample will also yield a different fit to the GS/E sample, one which is more likely to resemble the low-β stellar halo. We specifically attribute the large derived mass of the |$\sqrt{J_\mathrm{R}}-L_\mathrm{z}$| selection sample to this exact behaviour, since the purity using our original DF models is 0.78 and it drops to 0.52 when using our improved kinematic effective selection function.

5.3.2 The density profile of GS/E

We find that the density profile of GS/E is described by a shallow power law in the inner parts of the Galaxy, and steepens between 20 and 30 kpc. These findings are in good agreement with recent studies such as Han et al. (2022), who detect breaks at 12 and 28 kpc. Our power-law models with two breaks find those radii to be at about 15–20 kpc and 30–40 kpc, which is broadly consistent with the findings of Han et al. (2022). These authors find inner, middle, and outer power laws of 1.7, 3.1, and 4.6 respectively. These power laws are roughly consistent with both our BPL and DBPL findings for the ad selection, where for the BPL profile we have inner power laws of about 1.2 steepening to about 3.8, and for the DBPL profile the inner, middle, and outer power laws are about 0.9, 3.0, and 6.2. We are therefore shallower in the inner Galaxy and steeper in the outer Galaxy when compared with the results of Han et al. (2022). Our break radius is also consistent with the findings of MB20 who see a break radius of 25 kpc. They do find a much steeper inner power law (α ∼ 3.5), but they were not specifically targeting GS/E in their fits. in general, the finding of a break radius between 15–30 kpc is in broad agreement with other literature findings (Deason et al. 2011; Sesar, Jurić & Ivezić 2011; Xue et al. 2015; Deason et al. 2019). While we do detect power law breaks at reasonable radii, the broken power-law models with one or two breaks are never favoured compared with simpler models. It is likely that this is caused by the extreme paucity of stars in our GS/E samples and the relative lack of stars in our halo sample compared with, for example, the H3 data set (Conroy et al. 2019). Upcoming surveys such as SDSS-V (Kollmeier et al. 2017) will have vastly more stars with which to probe more complicated density profiles using techniques with restrictive samples such as we use here.

Interestingly, we find that our best-fitting density profiles for GS/E are similar in orientation to those found by Iorio & Belokurov (2019). Specifically, the major axis of our density profiles is rotated by 90°–100° (specifically |$99.5^{+13.9}_{-14.9}$| for our best-fitting ad model) in the Galactocentric X–Y plane towards the direction of Galactic rotation (positive Y), and inclined about 16° below the plane. Iorio & Belokurov (2019) find that the principal axis is rotated about 20° from the Galactocentric Y-axis (about 70° towards anti-rotation from the X-axis) and inclination of about 20° (oriented the same as our model). When we cast our rotation angles in the same frame as theirs we would find rotation of about 9.5° from the Y-axis and 80.5° towards anti-rotation from the X-axis. These results are very similar (cf. our Fig. 10 with their figs 2 and 3), reinforcing the idea that the Hercules-Aquila Cloud and Virgo Overdensity are related to the GS/E remnant. Similar findings were also recently reported by Han et al. (2022), although they find the rotation of the principal axis with respect to the Sun-Galactic centre axis to be only 24°, but a similar inclination of the principal axis of 25°.

Regarding axis scale ratios, we find that the GS/E remnant is significantly more triaxial than other studies have reported. Han et al. (2022) find axis ratios of 1:0.8:0.7 for the GS/E remnant, and Iorio et al. (2018) finds 1:0.8:0.5–0.6, while we find axis ratios 1:0.55:0.45. The cause of this mismatch could be attributed to our selections, and perhaps the GS/E is truly more triaxial than these other studies have noted. A more likely interpretation of the difference in these findings is that GS/E changes in shape with increasing Galactocentric radius. The model of Iorio et al. (2018) has variable flattening with radius, while the shape parameters of the models used by Han et al. (2022), like ours, are fixed with radius. It may be that surveys which are sensitive to different ranges of Galactocentric radius (due to the differences in observational footprint, choice of tracer, etc.) find different shape parameters for GS/E. It is curious though that the results of Han et al. (2022), using data from the H3 survey which should probe larger Galactocentric radii compared with APOGEE which probe smaller Galactocentric radii, should find less triaxiality since one might generally expect large degrees of triaxiality at larger Galactocentric radii. Nevertheless, when these results are collated with recent findings of complex structure associated with the GS/E remnant at large Galactocentric radii by Chandra et al. (2022), we can likely say that GS/E has a complicated shape-radius dependency.

6 SUMMARY AND CONCLUSIONS

In this work we have measured the mass of the GS/E remnant using high purity samples of APOGEE DR16 red giants selected using kinematics and chemistry. We employ a density modelling approach which accounts for the APOGEE survey selection function, dust obscuration, and the density of the tracer population in colour–magnitude space. To this method we add a novel technique for correcting the biases induced by kinematic selections which is based the creation of a kinematic effective selection function using fiducial distribution function models of the Milky Way stellar halo. We identify GS/E in multiple kinematic spaces using high-purity kinematic selections based on the work of LBM22, and we fit triaxial rotated density profiles of varying radial complexity to these samples. We validate our density modelling approach, and specifically our new kinematic effective selection functions, using mock APOGEE data. Our main results can be summarized as follows:

• We find that our kinematically selected GS/E samples are well fit by density profiles which have relatively shallow or flat power-law indices in the inner Galaxy (0.5 < α₁ < 1.5) and drop sharply at ∼ 20–30 kpc to α > 3.5 at larger radii.

• The best-fitting GS/E density profiles are triaxial, with axis ratios 1:0.55:0.45. The principal axis is rotated by about 100° towards the direction of Galactic rotation and is inclined 16° below the plane, which is in good agreement with previous studies.

• We determine the mass of GS/E to be |$1.45\ ^{+0.92}_{-0.51}\, 10^8\, M_\odot$|⁠, which is lower than found in previous studies. An assessment of sources of systematic error suggests they are smaller in magnitude than the statistical uncertainty quoted above. We attribute our finding of a lower mass to our use of restrictive kinematic cuts which produces a high-purity sample of GS/E stars and removes contamination from stellar halo populations with lower velocity anisotropy.

• A stellar mass of |$1.5\times 10^8\, \mathrm{M}_{\odot }$| corresponds to a total halo mass of GS/E of |$\approx 5-6\times 10^{10}\, \mathrm{M}_{\odot }$| using standard stellar mass-to-halo mass relations. Given that the Milky Way’s halo mass at the time of the merger was likely |$\approx 5\times 10^{11}\, \mathrm{M}_{\odot }$|⁠, GS/E constituted a minor 1:8–1:10 merger.

• We also fit density profiles to the whole stellar halo, finding that the best-fitting profiles have power-law index around α = 2.5, with no evidence for breaks in the density profile. We find the mass of the stellar halo to be between 6.7 × 10⁸  and  8.4 × 10⁸ M_⊙ out to 70 kpc. GS/E therefore represents only a modest fraction of the mass of the stellar halo, about 15 per cent at the low end or 25 per cent at the high end.

We have shown evidence that the GS/E remnant is lighter than many previous studies have suggested. We argue that our precise kinematic selections for GS/E driven by fiducial stellar halo distribution function models leads to less contamination by stellar halo populations with low-velocity anisotropy. While we are confident in our results, they represent the first application of a new technique for integrating kinematic selections into Galactic density modelling through the use of distribution functions. Additional work needs to be done to explore distribution functions which more accurately represent the stellar populations known to exist in the stellar halo, including triaxial, rotating distribution functions for example. We hope that simulations may be illuminating for this purpose, and plan to study this technique in the context of cosmological N-body simulations. We envision a future where stellar halo density modelling does not rely on kinematic selections to define samples, but involves modelling of the full six dimensions of phase space.

ACKNOWLEDGEMENTS

We first thank the referee for their comments, which have certainly improved the quality of the manuscript. JMML and JB acknowledge financial support from NSERC (funding reference number RGPIN-2020-04712) and an Ontario Early Researcher Award (ER16-12-061). We thank Josh Speagle and Ting Li for helpful discussions and comments on the manuscript. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for the Sloan Digital Sky Survey IV (SDSS IV) has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

DATA AVAILABILITY

The APOGEE DR16 data used in this article are available at: https://www.sdss.org/dr16. The Gaia data used in this article are available at: https://gea.esac.esa.int/archive/.

Footnotes

References

APPENDIX A: VALIDATION OF THE METHOD WITH MOCK DATA

Here we show the results of applying our density fitting technique to mock data, as described in Section 3.6. Fig. A1 shows the posterior distributions for fits to mock data with and without disc contamination, as well as mock data created using kinematic cuts that mirror those used to select the GS/E samples. The contours of the posterior distributions are set at 1σ and 2σ. The black lines show the input parameters for the models. The posterior distributions always contain the input model parameters, yet are broad in general. This is due to the fact that we chose to fit a lightweight mock halo (only 2 × 10⁸ M_⊙ total mass) to test whether our modelling approach would work with a small number of stars. Since each stellar halo realization is subject to a certain degree of randomness, we generate multiple stellar haloes with similar parameters and a variety of masses. We find that the input parameters are always contained in the posterior, and that there are no major systematic over- or underestimations of input parameters beyond those noted in Section 3.6.

APPENDIX B: COMPLETENESS OF THE MODIFIED KINEMATIC EFFECTIVE SELECTION FUNCTION

Here we show the kinematic effective selection function derived based on our fits to the GS/E and whole-halo populations described in Section 4.4. Figs B1 and B2 are identical to Figs 5 and 6, except they use this new kinematic effective selection function (denoted |$\mathfrak {S}_{2}^{\prime }$|⁠). On to Fig. B1 we add the median of |$\mathfrak {S}_{2}^{\prime }$| and |$\mathfrak {S}_{2}$|⁠, calculated across all fields at each distance modulus, to aid efficient comparison of typical values.