The chemical characterization of halo substructure in the Milky Way based on APOGEE Free

ABSTRACT

Galactic haloes in a Λ-CDM universe are predicted to host today a swarm of debris resulting from cannibalized dwarf galaxies. The chemodynamical information recorded in their stellar populations helps elucidate their nature, constraining the assembly history of the Galaxy. Using data from APOGEE and Gaia, we examine the chemical properties of various halo substructures, considering elements that sample various nucleosynthetic pathways. The systems studied are Heracles, Gaia-Enceladus/Sausage (GES), the Helmi stream, Sequoia, Thamnos, Aleph, LMS-1, Arjuna, I’itoi, Nyx, Icarus, and Pontus. Abundance patterns of all substructures are cross-compared in a statistically robust fashion. Our main findings include: (i) the chemical properties of most substructures studied match qualitatively those of dwarf Milky Way satellites, such as the Sagittarius dSph. Exceptions are Nyx and Aleph, which are chemically similar to disc stars, implying that these substructures were likely formed in situ; (ii) Heracles differs chemically from in situ populations such as Aurora and its inner halo counterparts in a statistically significant way. The differences suggest that the star formation rate was lower in Heracles than in the early Milky Way; (iii) the chemistry of Arjuna, LMS-1, and I’itoi is indistinguishable from that of GES, suggesting a possible common origin; (iv) all three Sequoia samples studied are qualitatively similar. However, only two of those samples present chemistry that is consistent with GES in a statistically significant fashion; (v) the abundance patterns of the Helmi stream and Thamnos are different from all other halo substructures.

1 INTRODUCTION

‘How did the Milky Way form?’ is likely the most fundamental question facing the Field of Galactic archaeology. When posed in a cosmological context, the Λ-CDM model predicts that the Galaxy formed in great measure via the process of hierarchical mass assembly. In this scenario, nearby satellite galaxies are consumed by the Milky Way due to them being attracted to its deeper gravitational potential, and as a result merge with the Galaxy. In such cases, these merger events shape the formation and evolution of the Milky Way. Therefore, an understanding of the assembly history of the Milky Way in the context of Λ-CDM depends critically on the determination of the properties of the systems accreted during the Galaxy’s history, including their masses and chemical compositions. Moreover, the merger history of the Galaxy has a direct impact on its resulting stellar populations at present time, and plays a vital role in shaping its components.

Since the seminal work by Searle & Zinn (1978), many studies have aimed at characterizing the stellar populations of the Milky Way, linking them to either an ‘in situ’ or accreted origin. Albeit detection of substructure in phase space has worked extremely well for the identification of on-going and/or recent accretion events (e.g. Sagittarius dSph, Ibata, Gilmore & Irwin 1994; Helmi stream, Helmi et al. 1999), the identification of accretion events early in the life of the Milky Way has proven difficult due to phase-mixing. A possible solution to this conundrum resides in the use of additional information, typically in the form of detailed chemistry and/or ages (e.g. Nissen & Schuster 2010; Hawkins et al. 2015; Hayes et al. 2018; Haywood et al. 2018; Mackereth et al. 2019; Das, Hawkins & Jofré 2020; Montalbán et al. 2020; Hasselquist et al. 2021; Horta et al. 2021; Buder et al. 2022; Carrillo et al. 2022).

The advent of large spectroscopic surveys such as APOGEE (Majewski et al. 2017), GALAH (De Silva et al. 2015), SEGUE (Yanny et al. 2009a), RAVE (Steinmetz et al. 2020), LAMOST (Zhao et al. 2012), H3 (Conroy et al. 2019), amongst others, in combination with the outstanding astrometric data supplied by the Gaia satellite (Gaia Collaboration 2018, 2021), revolutionized the field of Galactic archaeology, shedding new light into the mass assembly history of the Galaxy.

The core of the Sagittarius dSph system and its still forming tidal stream (Ibata et al. 1994) have long served as an archetype for dwarf galaxy mergers in the Milky Way. Moreover, in the past few years, several phase-space substructures have been identified in the Field of the Galactic stellar halo that are believed to be the debris of satellite accretion events, including the Gaia-Enceladus/Sausage (GE/S; Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Mackereth et al. 2019), Heracles (Horta et al. 2021), Sequoia (Barbá et al. 2019; Matsuno, Aoki & Suda 2019; Myeong et al. 2019), Thamnos 1 and 2 (Koppelman et al. 2019a), Nyx (Necib et al. 2020), LMS-1 (Yuan et al. 2020),¹ the substructures identified using the H3 survey: namely Aleph, Arjuna, and I’itoi (Naidu et al. 2020), Icarus (Re Fiorentin et al. 2021), Cetus (Newberg, Yanny & Willett 2009), and Pontus (Malhan et al. 2022). While the identification of these substructures is helping constrain our understanding of the mass assembly history of the Milky Way, their association with any particular accretion event still needs to be clarified. Along those lines, predictions from numerical simulations suggest that a single accretion event can lead us to multiple substructures in phase space (e.g. Jean-Baptiste et al. 2017; Koppelman, Bos & Helmi 2020). Therefore, in order to ascertain the reality and/or distinction of these accretion events, one must combine phase-space information with detailed chemical compositions for large samples.

Previous studies dedicated to characterizing the chemical properties of halo substructures have been based on either large samples from spectroscopic surveys, focusing on a relatively small number of elemental abundances (e.g. Hayes et al. 2020; Feuillet et al. 2021; Hasselquist et al. 2021; Buder et al. 2022), or more detailed studies of smaller samples from follow-up programs at an even higher resolution (e.g. Monty et al. 2020; Limberg et al. 2021; Matsuno et al. 2022a, b; Naidu et al. 2022). This is the first attempt at mapping the detailed abundance patterns of all halo substructures reported thus far in the literature.

In this work, we set out to combine the latest data releases from the APOGEE and Gaia surveys in order to dynamically determine and chemically characterize previously identified halo substructures in the Milky Way. We attempt, where possible, to define the halo substructures using kinematic information only, so that the distributions of stellar populations in various chemical planes can be studied in an unbiased fashion. This allows us to understand in more detail the reality and nature of these identified halo substructures, as chemical abundances encode more pristine fossilized records of the formation environment of stellar populations in the Galaxy.

For the convenience of the reader, we hope to facilitate navigation of the paper by supplying a detailed list of its contents, as follows:

• Section 2 describes the data and the selection criteria adopted to select the parent sample upon which this work is based. This is followed in Section 2.1 by a detailed account of the criteria adopted to define each substructure, building on techniques from previous work. For the keen reader interested solely on the study of the substructures’ chemical compositions, we suggest jumping straight into Section 4.

• Section 3 very briefly presents the resulting distributions of the various structures in the orbital energy versus angular momentum plane.

• Section 4 presents a qualitative examination of the halo substructures in different chemical planes. Specifically, we show the α elements in Section 4.1, the iron-peak elements in Section 4.2, the odd-Z elements in Section 4.3, the carbon and nitrogen abundances in Section 4.4, a neutron capture element (namely, Ce) in Section 4.5, and the [Mg/Mn]-[Al/Fe] chemical composition plane in Section 4.6. Readers interested in the quantitative results may skip this section and go straight to Section 5.

• Section 5 describes a statistical technique we developed to calculate a quantitative estimate of the chemical similitude between any two substructures. We run comparisons between all possible pairs of substructures and the results are encapsulated in the form of a confusion matrix in Fig. 18.

• We then discuss our results in the context of previous work in Section 6, and summarize our conclusions in Section 7.

2 DATA AND SAMPLE

This paper combines the latest data release (DR17; Abdurro’uf et al. 2021) of the SDSS-III/IV (Eisenstein et al. 2011; Blanton et al. 2017) and APOGEE survey (Majewski et al. 2017) with distances and astrometry determined from the early third data release from the Gaia survey (EDR3; Gaia Collaboration 2021). The celestial coordinates and radial velocities supplied by APOGEE (Nidever et al. 2015; Holtzman et al, in preparation), when combined with the proper motions and inferred distances (Leung & Bovy 2019b) based on Gaia data, provide complete 6D phase space information for over ∼700 000 stars in the Milky Way, for most of which exquisite abundances for up to ∼20 different elements have been determined.

All data supplied by APOGEE are based on observations collected by (almost) twin high-resolution multifibre spectrographs (Wilson et al. 2019) attached to the 2.5m Sloan telescope at Apache Point Observatory (Gunn et al. 2006) and the du Pont 2.5 m telescope at Las Campanas Observatory (Bowen & Vaughan 1973). Elemental abundances are derived from automatic analysis of stellar spectra using the ASPCAP pipeline (García Pérez et al. 2016) based on the ferre² code (Allende Prieto et al. 2006) and the line lists from Cunha et al. (2017) and Smith et al. (2021). The spectra themselves were reduced by a customized pipeline (Nidever et al. 2015). For details on target selection criteria, see Zasowski et al. (2013) for APOGEE, Zasowski et al. (2017) for APOGEE-2, Beaton et al. (2021) for APOGEE north, and Santana et al. (2021) for APOGEE south.

We make use of the distances for the APOGEE DR17 catalogue generated by Leung & Bovy (2019a), using the astroNN python package (for a full description, see Leung & Bovy 2019b). These distances are determined using a re-trained astroNN neural-network software, which predicts stellar luminosity from spectra using a training set comprised of stars with both APOGEE DR17 spectra and Gaia EDR3 parallax measurements (Gaia Collaboration 2021). The model is able to simultaneously predict distances and account for the parallax offset present in Gaia-EDR3, producing high precision, accurate distance estimates for APOGEE stars, which match well with external catalogues (Hogg, Eilers & Rix 2019) and standard candles like red clump stars (Bovy et al. 2014). We note that the systematic bias in distance measurements at large distances for APOGEE DR16 as described in Bovy et al. (2019) have been reduced drastically in APOGEE DR17. Therefore, we are confident that this bias will not lead to unforeseen issues during the calculation of the orbital parameters. Our samples are contained within a distance range of ∼20 kpc and have a mean d_err/d ∼0.13 (except for the Sagittarius dSph, which extends up to ∼30 kpc and has a mean d_err/d ∼0.16).

We use the 6D phase space information³ and convert between astrometric parameters and Galactocentric cylindrical coordinates, assuming a solar velocity combining the proper motion from Sgr A* (Reid & Brunthaler 2020) with the determination of the local standard of rest of Schönrich, Binney & Dehnen (2010). This adjustment leads to a 3D velocity of the Sun equal to [U_⊙, V_⊙, W_⊙]  = [–11.1, 248.0, 8.5] km s⁻¹. We assume the distance between the Sun and the Galactic Centre to be R₀  = 8.178 kpc (Gravity Collaboration 2019), and the vertical height of the Sun above the midplane z₀  = 0.02 kpc (Bennett & Bovy 2019). Orbital parameters were then determined using the publicly available code galpy⁴ (Bovy 2015; Mackereth & Bovy 2018), adopting a McMillan (2017) potential and using the Stäckel approximation of Binney (2012).

The parent sample employed in this work is comprised of stars that satisfy the following selection criteria:

• APOGEE-determined atmospheric parameters: 3500 < T_eff < 5500 K and log g < 3.6,

• APOGEE spectral S/N > 70,

• APOGEE STARFLAG  = 0,

• astroNN distance accuracy of |$d_{\odot , \rm err}/d_{\odot }\lt 0.2$|⁠,

where d_⊙ and |$d_{\odot \rm err}$| are heliocentric distance and its uncertainty, respectively. The S/N criterion was implemented to maximize the quality of the elemental abundances. The T_eff and log g criteria aimed to minimize systematic effects at high/low temperatures, and to minimize contamination by dwarf stars. We also removed stars with STARFLAG flags set, in order to not include any stars with issues in their stellar parameters. A further 7750 globular cluster stars were also removed from consideration using the APOGEE Value Added Catalogue of globular cluster candidate members from Schiavon et al. (in preparation), (building on the method from Horta et al. 2020, using primarily radial velocity and proper motion information). Finally, stars belonging to the Large and Small Magellanic clouds were also excluded using the sample from Hasselquist et al. (2021) (removing 3748 and 1002 stars, respectively). The resulting parent sample contains 199 030 stars.

In the following subsection, we describe the motivation behind the selection criteria employed to select each substructure in the stellar halo of the Milky Way. The criteria are largely built on selections employed in previous works and are summarized in Table 1 and in Fig. 1.

2.1 Identification of substructures in the stellar halo

We now describe the method employed for identifying known substructures in the stellar halo of the Milky Way. We set out to select star members belonging to the various halo substructures.

We strive to identify substructures in the stellar halo by employing solely orbital parameter and phase-space information where possible, with the aim of obtaining star candidates for each substructure population that are unbiased by any chemical composition selection.

We take a handcrafted approach and select substructures based on simple and reproducible selection criteria that are physically motivated by the data and/or are used in previous works, instead of resorting to clustering software algorithms, which we find cluster the n-dimensional space into too many fragments and vary wildly depending on the input parameters to the clustering model. In the following subsections, we describe the selection procedure for identifying each substructure independently.

We note that our samples for the various substructures are defined by a strict application to the APOGEE survey data of the criteria defined by other groups, often on the basis of different data sets. The latter were per force collected as part of a different observational effort, based on specific target selection criteria. It is not immediately clear whether or how differences between the APOGEE selection function and those of other catalogues may imprint dissimilarities between our samples and those of the original studies. We nevertheless do not expect such effects to influence our conclusions, regarding the chemical compositions of these structures, in an important way.

2.1.1 Sagittarius

Since its discovery (Ibata et al. 1994), many studies have sought to characterize the nature of the Sagittarius dwarf spheroidal (hereafter Sgr dSph; e.g. Ibata et al. 2001; Majewski et al. 2003; Johnston, Law & Majewski 2005; Belokurov et al. 2006; Yanny et al. 2009b; Koposov et al. 2012; Carlin et al. 2018; Antoja et al. 2020; Ibata et al. 2020; Vasiliev & Belokurov 2020), as well as interpret its effect on the Galaxy using numerical simulations (e.g. Johnston, Spergel & Hernquist 1995; Ibata et al. 1997; Law, Johnston & Majewski 2005; Law & Majewski 2010; Purcell et al. 2011; Gómez et al. 2013). More recently, the Sgr dSph has been the subject of comprehensive studies on the basis of APOGEE data. This has enabled a detailed examination of its chemical compositions, both in the satellite’s core and in its tidal tails (e.g. Hasselquist et al. 2017, 2019; Hayes et al. 2020). Moreover, in a more recent study, Hasselquist et al. (2021) adopted chemical evolution models to infer the history of star formation and chemical evolution of the Sgr dSph. Therefore, in this paper the Sgr dSph is considered simply as a template massive satellite whose chemical properties can be contrasted to those of the halo substructures that are the focus of our study.

While it is possible to select high confidence Sgr dSph star candidates using Galactocentric positions and velocities (Majewski et al. 2003), Hayes et al. (2020) showed it is possible to make an even more careful selection by considering the motion of stars in a well-defined Sgr orbital plane. We identify Sgr star members by following the method from Hayes et al. (2020). Although the method is fully described in their work, we summarize the key steps for clarity and completeness. We take the Galactocentric positions and velocities of stars in our parent sample and rotate them into the Sgr orbital plane according to the transformations described in Majewski et al. (2003). This yields a set of position and velocity coordinates relative to the Sgr orbital plane, but still centred on the Galactic Centre. As pointed out in Hayes et al. (2020), Sgr star members should stand out with respect to other halo populations in different Sgr orbital planes. Using this orbital plane transformation, we select from the parent sample Sgr star members if they satisfy the following selection criteria:

• |β_GC| < 30 (°),

• 18 < L_{z, Sgr} < 14 (× 10³ kpc km s⁻¹),

• −150 < V_{z, Sgr} < 80 (km s⁻¹),

• X_Sgr > 0 or X_Sgr < −15 (kpc),

• Y_Sgr > −5 (kpc) or Y_Sgr < −20 (kpc),

• Z_Sgr > −10 (kpc),

• pm_α > −4 (mas),

• d_⊙ > 10 (kpc),

where β_GC is the angle subtended between the Galactic Centre and the Sgr dSph, L_{z, Sgr} is the azimuthal component of the angular momentum in the Sgr plane, V_{z, Sgr} is the vertical component of the velocity in the Sgr plane, (X_Sgr, Y_Sgr, Z_Sgr) are the Cartesian coordinates centred on the Sgr dSph plane, pm_α is the right-ascension proper motion, and d_⊙ is the heliocentric distance, which for Sgr has been shown to be ∼ 23 kpc (Vasiliev & Belokurov 2020) (although we follow the distance cut from Hayes et al. (2020) for this work, to prevent the exclusion of more nearby parts of the Sgr stream). Our selection yields a sample of 266 Sgr star members, illustrated in the L_{z, Sgr} versus V_{z, Sgr} plane in Fig. 2.

2.1.2 Heracles

Heracles is a recently discovered metal-poor substructure located in the heart of the Galaxy (Horta et al. 2021). It is characterized by stars on eccentric and low energy orbits. Due to its position in the inner few kpc of the Galaxy, it is highly obscured by dust extinction and vastly outnumbered by its more abundant metal-rich (in situ) co-spatial counterpart populations. Only with the aid of chemical compositions has it been possible to unveil this metal-poor substructure, which is discernible in the [Mg/Mn]-[Al/Fe] plane. It is important at this stage that we mention a couple of recent studies which, based chiefly on the properties of the Galactic globular cluster system, proposed the occurrence of an early accretion event whose remnants should have similar properties to those of Heracles (named Kraken and Koala, by Kruijssen et al. 2020; Forbes 2020, respectively). In the absence of a detailed comparison of the dynamical properties and detailed chemical compositions of Heracles with these putative systems, and recent results by Pagnini et al. (2022), a definitive association is impossible at the current time.

In this work, we define Heracles candidate star members following the work by Horta et al. (2021), and select stars from our parent sample that satisfy the following selection criteria:

• e > 0.6,

• −2.6 < E < −2 (× 10⁵ km² s⁻²),

• [Al/Fe] < −0.07 |$\&$| [Mg/Mn] ≥ 0.25,

• [Al/Fe] ≥ −0.07 |$\&$| [Mg/Mn] ≥ 4.25 × [Al/Fe]  + 0.5475.

Moreover, we impose a [Fe/H] > −1.7 cut to select Heracles candidate star members in order to select stars from our parent sample that have reliable Mn abundances in APOGEE DR17. Our selection yields a resulting sample of 300 Heracles star members.

2.1.3 Gaia-Enceladus/Sausage

Recent studies have shown that there is an abundant population of stars in the nearby stellar halo (namely, R_GC ≲ 20–25 kpc) belonging to the remnant of an accretion event dubbed the Gaia-Enceladus/Sausage (GES; e.g. Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Mackereth et al. 2019). This population is characterized by stars on highly radial/eccentric orbits, which also appear to follow a lower distribution in the α-Fe plane, presenting lower [α/Fe] values for fixed metallicity than in situ populations.

For this paper, we select GES candidate star members by employing a set of orbital information cuts. Specifically, GES members were selected adopting the following criteria:

• |L_z| < 0.5 (× 10³ kpc km⁻¹),

• −1.6 < E < −1.1 (× 10⁵ km² s⁻²).

This selection is employed in order to select the clump that becomes apparent in the E–L_z plane at higher orbital energies and roughly L_z ∼ 0 (see Fig. 4), and to minimize the contamination from high-α disc stars on eccentric orbits (namely, the ‘Splash’ Bonaca et al. 2017; Belokurov et al. 2020), which sit approximately at E∼–1.8 × 10⁵ km² s⁻² (see Kisku et al, in preparation). The angular momentum restriction ensures we are not including stars on more prograde/retrograde orbits. We find that by selecting the GES substructure in this manner, we obtain a sample of stars with highly radial (J_R ∼1×10³ kpc km s⁻¹) and therefore highly eccentric (e ∼ 0.9) orbits, in agreement with selections employed in previous studies to identify this halo substructure (e.g. Mackereth & Bovy 2020; Naidu et al. 2020; Feuillet et al. 2021; Buder et al. 2022). The final GES sample is comprised of 2353 stars.

We note that in a recent paper, Hasselquist et al. (2021) undertook a thorough investigation into the chemical properties of this halo substructure and compared it to other massive satellites of the Milky Way (namely, the Magellanic Clouds, Sagittarius dSph, and Fornax). Although their selection criteria differs slightly from the one employed in this study, we find that their sample is largely similar to the one employed here, as both studies employed APOGEE DR17 data.

2.1.4 Retrograde halo: Sequoia, Thamnos, Arjuna, and I’Itoi

A number of substructures have been identified in the retrograde halo. The first to be discovered was Sequoia (Barbá et al. 2019; Matsuno et al. 2019; Myeong et al. 2019), which was suggested to be the remnant of an accreted dwarf galaxy. The Sequoia was identified given the retrograde nature of the orbits of its stars, which appear to form an arch in the retrograde wing of the Toomre diagram. Separately, an interesting study by Koppelman et al. (2019a) showed that the retrograde halo can be further divided into two components, separated by their orbital energy values in the E–L_z plane. They suggest that the high energy component corresponds to Sequoia, whilst the lower energy population would be linked to a separate accretion event, dubbed Thamnos. In addition, Naidu et al. (2020) proposed the existence of additional retrograde substructure overlapping with Sequoia, characterized by different metallicities, which they named Arjuna and I’Itoi.

As the aim of this paper is to perform a comprehensive study of the chemical abundances of substructures identified in the halo, we utilize all the selection methods used in prior work and select the same postulated substructures in multiple ways, in order to compare their abundances later. Specifically, we build on previous works (e.g. Koppelman et al. 2019a; Myeong et al. 2019; Naidu et al. 2020) that have aimed to characterize the retrograde halo and select the substructures following a similar selection criteria.

For reference, throughout this work we will refer to the different selection criteria of substructures in the retrograde halo as the ‘M19’, ‘K19’, and ‘N20’ selections, in reference to the Myeong et al. (2019), Koppelman et al. (2019c), and Naidu et al. (2020) studies that determined the Sequoia substructure, respectively. We will now go through the details of each selection method independently.

The M19 method (used in Myeong et al. 2019) selects Sequoia star candidates by identifying stars that satisfy the following conditions:

• E > −1.5 (× 10⁵ km² s⁻²),

• J_ϕ/J_tot < −0.5,

• J|$_{\mathrm{(J_{z}-J_{R})}}$|/J_tot <0.1.

Here, J_ϕ, J_R, and J_z are the azimuthal, radial, and vertical actions, and J_tot is the quadrature sum of those components. This selection yields a total of 116 Sequoia star candidates.

The K19 method (used in Koppelman et al. 2019c) identifies Sequoia star members based on the following selection criteria:

• −0.65 < η < −0.4,

• −1.35 < E < −1 (× 10⁵ km² s⁻²),

• L_z < 0 (kpc km s⁻¹),

where η is the circularity. These selection criteria yield a total of 45 Sequoia stars.

Lastly, we select the Sequoia based on the N20 selection criteria (used in Naidu et al. 2020) as follows:

• η < −0.15,

• E > −1.6 (× 10⁵ km² s⁻²),

• L_z < −0.7 (10³ kpc km s⁻¹).

This selection produces a preliminary sample comprised of 478 stars. However, we note that by purely selecting the high-energy retrograde structures in APOGEE DR17-Gaia-DR3 as done so in Naidu et al. (2020) leads to some contamination from Magellanic Cloud stars not included in the sample from Hasselquist et al. (2021) that we removed earlier. We further exclude these stars by performing an additional [Fe/H]>–0.7 and d_⊙ < 20 kpc cut. This leads to a high-energy retrograde sample comprising of 227 Sequoia stars. Furthermore, Naidu et al. (2020) used this selection to define not only Sequoia, but all the substructures in the high-energy retrograde halo (including the Arjuna and I’itoi substructures). In order to distinguish Sequoia, I’itoi, and Arjuna, Naidu et al. (2020) suggest performing a metallicity cut, motivated by the observed peaks in the metallicity distribution function (MDF) of their retrograde sample. Thus, we follow this procedure and further refine our Sequoia, Arjuna, and I’itoi samples by requiring −1.6 < [Fe/H] cut for Arjuna, −2 < [Fe/H] < −1.6 for Sequoia, and [Fe/H] < −2 for I’itoi, based on the distribution of our initial sample in the MDF (see Fig. 3). This further division yields an N20 Sequoia sample comprised of 72 stars, an Arjuna sample constituted by 132 stars, and an I’itoi sample comprised of 25 stars.

Following our selection of substructures in the high-energy retrograde halo, we set out to identify stars belonging to the intermediate-energy and retrograde Thamnos 1 and 2 substructures. Koppelman et al. (2019c) state that Thamnos 1 and 2 are separate debris from the same progenitor galaxies. For this work, we consider Thamnos as one overall structure, given the similarity noted by Koppelman et al. (2019c) between the two smaller individual populations in chemistry and kinematic planes. Stars from our parent sample were considered as Thamnos candidate members if they satisfied the following selection criteria:

• −1.8 < E < −1.6 (× 10⁵ km² s⁻²),

• L_z < 0 (kpc km s⁻¹),

• e < 0.7,

These selection cuts are performed in order to select stars in our parent sample with intermediate orbital energies and retrograde orbits (see Fig. 4 for the position of Thamnos in the E-L_z), motivated by the distribution of the Thamnos substructure in the E–L_z plane illustrated by Koppelman et al. (2019c). This selection yields a Thamnos sample comprised of 121 stars.

2.1.5 Helmi stream

The Helmi stream was initially identified as a substructure in orbital space due to its high V_z velocities (Helmi et al. 1999). More recent work by Koppelman et al. (2019b) characterized the Helmi stream in Gaia DR2, and found that this stellar population is best defined by adopting a combination of cuts in different angular momentum planes. Specifically, by picking stellar halo stars based on the azimuthal component of the angular momentum (L_z), and its perpendicular counterpart (L_⊥  = |$\sqrt{L_{x}^{2}+L_{y}^{2}}$|⁠), the authors were able to select a better defined sample of Helmi stream star candidates. We build on the selection criteria from Koppelman et al. (2019b) and define our Helmi stream sample by including stars from our parent population that satisfy the following selection criteria:

• 0.75 < L_z < 1.7 (× 10³ kpc km s⁻¹),

• 1.6 <L_⊥ < 3.2 (× 10³ kpc km s⁻¹).

Our final sample is comprised of 85 Helmi stream stars members.

2.1.6 Aleph

Aleph is a newly discovered substructure presented in a detailed study of the Galactic stellar halo by Naidu et al. (2020) on the basis of the H3 survey (Conroy et al. 2019). It was initially identified as a sequence below the high α-disc in the α-Fe plane. It is comprised by stars on very circular prograde orbits. For this paper, we follow the method described in Naidu et al. (2020) and define Aleph star candidates as any star in our parent sample which satisfies the following selection criteria:

• 175 < V_ϕ < 300 (km s⁻¹),

• |V_R| < 75 (km s⁻¹),

• [Fe/H] > −0.8,

• [Mg/Fe] < 0.27,

where V_ϕ and V_R are the azimuthal and radial components of the velocity vector (in Galactocentric cylindrical coordinates), and we use Mg as our α tracer element. The selection criteria yield a sample comprised of 128 578 stars. We find that the initial selection criteria determine a preliminary Aleph sample that is dominated by in situ disc stars, likely obtained due to the prograde nature of the velocity cuts employed as well as the chemical cuts. Thus, in order to remove disc contamination and select true Aleph star members, we employ two further cuts in vertical height above the plane (namely, |z| > 3 kpc) and in vertical action (i.e. 170 < J_z < 210 kpc km s⁻¹), which are motivated by the distribution of Aleph in these coordinates (see Section 3.2.2 from Naidu et al. 2020). We also note that the vertical height cut was employed in order to mimic the H3 survey selection function. After including these two further cuts, we obtain a final sample of Aleph stars comprised of 28 star members.

2.1.7 LMS-1

LMS-1 is a newly identified substructure discovered by Yuan et al. (2020). It is characterized by metal-poor stars that form an overdensity at the foot of the omnipresent GES in the E–L_z plane. This substructure was later also studied by Naidu et al. (2020), who referred to it as Wukong. We identify stars belonging to this substructure adopting a similar selection as Naidu et al. (2020), however adopting different orbital energy criteria to adjust for the fact that we adopt the McMillan (2017) Galactic potential (see fig. 23 from appendix B in that study). Stars from our parent sample were deemed LMS-1 members if they satisfied the following selection criteria:

• 0.2 < L_z < 1 (× 10³ kpc km s⁻¹),

• −1.7 < E < −1.2 (× 10⁵ km⁻² s⁻²),

• [Fe/H] < −1.45,

• 0.4 < e < 0.7,

• |z| > 3 (kpc).

We note that the e and z cuts were added to the selection criteria listed by Naidu et al. (2020). This is because we conjectured that instead of eliminating GES star members from our selection (as Naidu et al. (2020) do), it is more natural in principle to find additional criteria that distinguishes these two overlapping substructures. Thus, we select stars on less eccentric orbits than those belonging to GES, but still more eccentric than most of the Galactic disc (i.e. 0.4 < e < 0.7). Furthermore, in order to ensure that we are observing stars at the same distances above the Galactic plane as in Naidu et al. (2020) (defined by the selection function of the H3 survey), we add a vertical height cut of |z| > 3 (kpc). Our selection identifies 31 stars belonging to the LMS-1 substructure.

2.1.8 Nyx

Nyx has recently been put forward by Necib et al. (2020), who identified a stellar stream in the solar neighbourhood, that they suggest to be the remnant of an accreted dwarf galaxy (Necib et al. 2020). Similar to Aleph, it is characterized by stars on very prograde orbits, at relatively small mid-plane distances (|Z| < 2 kpc) and close to the solar neighbourhood (i.e. |Y| < 2 kpc and |X| < 3 kpc). The Nyx structure is also particularly metal-rich (i.e. [Fe/H] ∼ −0.5). Based on the selection criteria used in Necib et al. (2020), we select Nyx star candidates employing the following selection criteria:

• 110 < V_r < 205 (km s⁻¹),

• 90 < V_ϕ < 195 (km s⁻¹),

• |X| < 3 (kpc), |Y| < 2 (kpc), |Z| < 2 (kpc).

The above selection criteria yield a sample comprising of 589 Nyx stars.

2.1.9 Icarus

Icarus is a substructure identified in the solar vicinity, comprised by stars that are significantly metal-poor ([Fe/H] ∼ −1.45) with circular (disc-like) orbits (Re Fiorentin et al. 2021). In this work, we select Icarus star members using the mean values reported by those authors and adopting a two sigma uncertainty cut around the mean. The selection used is listed as follows:

• [Fe/H] < −1.05,

• [Mg/Fe] < 0.2,

• 1.54 < L_z < 2.21 (× 10³ kpc km s⁻¹),

• L_⊥ < 450 (kpc km s⁻¹),

• e < 0.2,

• z_max < 1.5 (kpc).

These selection criteria yield an Icarus sample comprised of one star. As we have only been able to identify one star associated with this substructure, we remove it from the main body of this work and focus on discerning why our selection method only identifies 1 star in Appendix  A. Furthermore, we combine the one Icarus star found in APOGEE DR17 with 41 stars found by Re Fiorentin et al. (2021) in APOGEE DR16, in order to study the nature of this substructure in further detail. Our results are discussed in Appendix  A.

2.1.10 Pontus

Pontus is a halo substructure recently proposed by Malhan et al. (2022), on the basis of an analysis of Gaia EDR3 data for a large sample of Galactic globular clusters and stellar streams. These authors identified a large number of groupings in action space, associated with known substructures. Malhan et al. (2022) propose the existence of a previously unknown substructure they call Pontus, characterized by retrograde orbits and intermediate orbital energy. Pontus is located just below Gaia-Enceladus/Sausage in theE–L_z plane, but displays less radial orbits (Pontus has an average radial action of J_R ∼ 500 kpc km s⁻¹, whereas Gaia-Enceladus/Sausage displays a mean value of J_R ∼ 1250 kpc km s⁻¹). In this work, we utilize the values listed in Section 4.6 from Malhan et al. (2022) to identify Pontus candidate members in our sample. We note that because both that study and ours use the McMillan (2017) potential to compute the IoM, the orbital energy values will be on the same scale. Our selection criteria for Pontus are the following:

• −1.72 < E < −1.56 (× 10⁵ km⁻² s⁻²),

• −470 < L_z < 5 (× 10³ kpc km s⁻¹),

• 245 < J_R < 725 (kpc km s⁻¹),

• 115 < J_z < 545 (kpc km s⁻¹),

• 390 < L_⊥ < 865 (kpc km s⁻¹),

• 0.5 < e < 0.8,

• 1 < R_peri < 3 (kpc),

• 8 < R_apo < 13 (kpc),

• [Fe/H] < −1.3,

where R_peri and R_apo are the perigalacticon and apogalacticon radii, respectively. Using these selection criteria, we identify two Pontus candidate members in our parent sample. As two stars comprise a sample too small to perform any statistical comparison, we refrain from comparing the Pontus stars in the main body of this work. Instead, we display and discuss the chemistry of these two Pontus stars in Appendix  B, for completeness.

2.1.11 Cetus

As a closing remark, we note that we attempted to identify candidate members belonging to the Cetus (Newberg et al. 2009) stream. Using the selection criteria defined in table 3 from Malhan et al. (2022), we found no stars associated with this halo substructure that satisfied the selection criteria of our parent sample. This is likely due to a combination of two factors: (i) Cetus is a diffuse stream orbiting at large heliocentric distances (d_⊙ ≳ 30 kpc; Newberg et al. 2009), which APOGEE does not cover well; (ii) it occupies a region of the sky around the southern polar cap, where APOGEE does not have many field pointings, at approximately l ∼ 143° and b ∼ −70° (Newberg et al. 2009).

3 KINEMATIC PROPERTIES

In this Section, we present the resulting distributions of the identified halo substructures in the orbital energy (E) versus the azimuthal component of the angular momentum (L_z) plane in Fig. 4. The parent sample is illustrated as a density distribution and the halo substructures are shown using the same colour markers as in Fig. 1. By construction, each substructure occupies a different locus in this plane. However, we do notice some small overlap between some of the substructures (for example, between GES and Sequoia, or GES and LMS-1), given their similar selection criteria. More specifically, we find that Heracles dominates the low energy region (E < −2 × 10⁵ km⁻² s⁻²), whereas all the other substructures are characterized by higher energies. As shown before (e.g. Koppelman et al. 2019c; Horta et al. 2021), we find that GES occupies a locus at low L_z and relatively high E, which corresponds to very radial/eccentric orbits. We find the retrograde region (i.e. L_z < 0 × 10³ kpc km s⁻¹) to be dominated by Thamnos at intermediate energies (E ∼ −1.7 × 10⁵ km⁻² s⁻²), and by Sequoia, Arjuna, and I’itoi at higher energies (E > −1.4 × 10⁵ km⁻² s⁻²); on the other hand, in the prograde region (L_z > 0 × 10³ kpc km s⁻¹), we find the LMS-1 and Helmi stream structures, which occupy a locus at approximately E ∼ −1.5 × 10⁵ km⁻² s⁻² and L_z ∼ 500 kpc km s⁻¹, and E ∼ −1.4 × 10⁵ km⁻² s⁻² and L_z ∼ 1000 kpc km s⁻¹, respectively. Furthermore, the loci occupied by the Aleph and Nyx substructures closely mimic the region defined by disc orbits. Lastly, sitting above all other structures we find the Sagittarius dSph, which occupies a position at high energies and spans a range of angular momentum between 0 < L_z < 2000 kpc km s⁻¹.

4 CHEMICAL COMPOSITIONS

In this Section, we turn our attention to the main focus of this work: a chemical abundance study of substructures in the stellar halo of the Milky Way. In this Section, we seek to first characterize these substructures qualitatively in multiple chemical abundance planes that probe different nucleosynthetic pathways. In Section 5 we then compare mean chemical compositions across various substructures in a quantitative fashion. Our aim is to utilize the chemistry to further unravel the nature and properties of these halo substructures, and in turn place constraints on their star formation and chemical enrichment histories. We also aim to compare their chemical properties with those from in situ populations (see Fig. 5 for how we determine in situ populations). By studying the halo substructures using different elemental species we aim to develop an understanding of their chemical evolution contributed by different nucleosynthetic pathways, contributed either by core-collapse supernovae (SNII), supernovae type Ia (SNIa), and/or Asymptotic Giant Stars (AGBs). Furthermore, as our method for identifying these substructures relies mainly on phase space and orbital information, our analysis is not affected by chemical composition biases (except for the case of particular elements in the Heracles, Aleph, LMS-1, Arjuna, I’itoi, and (H3) Sequoia sample).

Our results are presented as follows. In Section 4.1 we present the distribution of the halo substructures in the α-Fe plane, using Mg as our α element tracer. In Section 4.2, we show the distribution of these substructures in the Ni-Fe abundance plane, which provides insight into the chemical evolution of the iron-peak elements. Section 4.3 displays the distribution of the halo substructures in an odd-Z-Fe plane, where we use Al as our tracer element. Furthermore, we also show the C and N abundance distributions in Section 4.4, the Ce abundances (namely, an s-process neutron capture element) in Section 4.5, and the distribution of the halo substructures in the [Mg/Mn]-[Al/Fe] plane in Section 4.6.⁵ This last chemical composition plane is interesting to study as it has recently been shown to help us distinguish stellar populations from ‘in situ’ and accreted origins (e.g. Hawkins et al. 2015; Das et al. 2020; Horta et al. 2021). Upon studying the distribution of the substructures in different chemical composition planes, we finalize our chemical composition study in Section 5 by performing a quantitative comparison between the substructures studied in this work for all the (reliable) elemental abundances available in APOGEE. For simplicity, we henceforth refer to the [X/Fe]-[Fe/H] plane as just the X-Fe plane.

As mentioned in Section 2, we exclude the Pontus and Icarus substructures from our quantitative chemical comparisons as the number of candidate members of these substructures in the APOGEE catalogue is too small. The properties of Icarus and Pontus are briefly discussed in Appendices  A and  B, respectively.

4.1 α-elements

We first turn our attention to the distribution of the substructures in the α-Fe plane. This is possibly the most interesting chemical composition plane to study, as it can provide great insight into the star formation history and chemical enrichment processes of each substructure (e.g. Matteucci & Greggio 1986; Wheeler, Sneden & Truran 1989; McWilliam 1997; Tolstoy, Hill & Tosi 2009; Nissen & Schuster 2010; Bensby, Feltzing & Oey 2014). Specifically, we seek to identify the presence of the α-Fe knee. For this work, we resort to magnesium as our primary α element, as this has been shown to be the most reliable α element in previous APOGEE data releases (e.g. DR16; Jönsson et al. 2020). For the distributions of the remaining α elements determined by ASPCAP (namely, O, Si, Ca, S, and Ti), we refer the reader to Figs D1–D5 in Appendix  D.

Fig. 6 shows the distribution of each substructure in the Mg–Fe plane (coloured markers) compared to the parent sample (2D density histogram). We find that all the substructures − except for Aleph and Nyx – occupy a locus in this plane which is typical of low mass satellite galaxies and accreted populations of the Milky Way (e.g. Tolstoy et al. 2009; Hayes et al. 2018; Mackereth et al. 2019), characterized by low metallicity and lower [Mg/Fe] at fixed [Fe/H] than in situ disc populations. Moreover, we find that different substructures display distinct [Mg/Fe] values, implying certain differences despite their overlap in [Fe/H]. However, we also note that at the lowest metallicities ([Fe/H] < −1.8), the overlap between different halo substructures increases.

Next we discuss the distribution on the Mg–Fe plane of stars in our sample belonging to each halo substructure.

• Sagittarius: As shown in previous studies (e.g. Monaco et al. 2005; Sbordone et al. 2007; Carretta et al. 2010; McWilliam, Wallerstein & Mottini 2013; Hasselquist et al. 2017, 2019; Hayes et al. 2020; Minelli et al. 2021), the stellar populations of the Sgr dSph galaxy are characterized by substantially lower [Mg/Fe] than even the low-α disc at the same [Fe/H] and traces a tail towards higher [Mg/Fe] with decreasing [Fe/H]. Conversely, at the higher [Fe/H] end, we find that Sgr stops decreasing in [Mg/Fe] and shows an upside-down ‘knee’, likely caused by a burst in SF at late times that may be due to its interaction with the Milky Way (e.g. Hasselquist et al. 2017, 2021). Such a burst of star formation intensifies the incidence of SNe II, with a consequent boost in the ISM enrichment in its nucleosynthetic by-products, such as Mg, over a short time-scale (∼10⁷ yr). Because Fe is produced predominantly by SN Ia over a considerably longer time-scale (∼10⁸–10⁹ yr), Fe enrichment lags behind, causing a sudden increase in [Mg/Fe].

• Heracles: The [Mg/Fe] abundances of the Heracles structure occupy a higher locus than that of other halo substructures of similar metallicity, with the exception of Thamnos. As discussed in Horta et al. (2021), the distribution of Heracles in the α-Fe plane is peculiar, differing from that of GES and other systems by the absence of the above-mentioned α-knee. As conjectured in Horta et al. (2021), we suggest that this distribution results from an early quenching of star formation, taking place before SN Ia could contribute substantially to the enrichment of the interstellar medium (ISM).

• Gaia-Enceladus/Sausage: Taking into account the small yet clear contamination from the high-α disc at higher [Fe/H] (see Fig. 7 for details), the distribution of GES dominates the metal-poor and α-poor populations of the Mg-Fe plane (as pointed out in previous studies e.g. Hayes et al. 2018; Helmi et al. 2018; Mackereth et al. 2019), making it easily distinguishable from the high-/low-α discs. We find that GES reaches almost solar metallicities, displaying the standard distribution in the Mg–Fe plane with a change of slope − the so-called ‘α-knee’ – occurring at approximately [Fe/H] ∼ −1.2 (Mackereth et al. 2019). The metallicity of the ‘knee’ has long been thought to be an indicator of the mass of the system (e.g. Tolstoy et al. 2009), and indeed it occurs at [Fe/H]>–1 for both the high- and low-α discs. As a result, GES stars in the shin part of the α-knee are characterized by lower [Mg/Fe] at constant metallicity than disc stars. Interestingly, even at the plateau ([Fe/H]<–1.2), GES seems to present lower [Mg/Fe] than the high-α disc, although this needs to be better quantified. Furthermore, we note the presence of a minor population of [Mg/Fe] < 0 stars at −1.8 < [Fe/H] < −1.2, which could be contamination from a separate halo substructure, possibly even a satellite of the GES progenitor (see Fig. 7 for details).

  Based purely on the distribution of its stellar populations on the Mg–Fe plane, one would expect the progenitor of GES to be a relatively massive system (see Mackereth et al. 2019, for details). The fact that the distribution of its stellar populations in the Mg–Fe plane covers a wide range in metallicity, bracketing the knee and extending from [Fe/H]<–2 all the way to [Fe/H] ∼ −0.5 suggests a substantially prolonged history of star formation.

• Sequoia: The distribution of the M19, K19, and N20 selected samples (selected on the criteria described in Myeong et al. 2019; Koppelman et al. 2019c; and Naidu et al. 2020, respectively) occupy similar [Mg/Fe] values ([Mg/Fe] ∼ 0.1) at lower metallicities ([Fe/H] ≲–1). More specifically, we find that all three Sequoia samples occupy a similar position in the Mg–Fe plane, one that overlaps with that of GES and other substructures at similar [Fe/H] values. Along those lines, we find that the K19 and M19 Sequoia samples seem to connect with the N20 Sequoia sample, where the N20 sample comprises the lower metallicity component of the K19/M19 samples. We examine in more detail the distribution of the Sequoia stars in the α-Fe plane in Section 6.1.3.

• Helmi stream: Despite the lower number of members associated to this substructure, its chemical composition in the Mg–Fe plane appears to follow a single sequence, with slope and scatter similar to that of GES, with the addition of a small number of likely contaminants overlapping with the low- and high-alpha disc. Furthermore, it is confined to low metallicities ([Fe/H]<–1.2) and intermediate magnesium values ([Mg/Fe] ∼ 0.2). However, we do note that the stars identified for this substructure appear to be scattered across a wide range of [Fe/H] values. Specifically, the Helmi stream occupies a locus that overlaps with the GES for fixed [Fe/H]. In Fig. 20 we show that the best-fitting piece-wise linear model prefers a knee that is ‘inverted’ similar to, although less extreme than, the case of Sgr dSph. We discuss this in more detail in Section 6.1.4.

• Thamnos: The magnesium abundances of Thamnos suggest that this structure is clearly different from other substructures in the retrograde halo (namely, Sequoia, Arjuna, and I’itoi). It presents a much higher mean [Mg/Fe] for fixed metallicity than the other retrograde substructures, and appears to follow the Mg–Fe relation of the high-α disc. We find that Thamnos presents no α-knee feature, and occupies a similar locus in the Mg–Fe plane to that of Heracles. The distribution of this substructure in this plane with the absence of an α-Fe ‘knee’ suggests that this substructure likely quenched star formation before the onset of SN Ia.

• Aleph: By construction, Aleph occupies a locus in the Mg–Fe plane that overlaps with the metal-poor component of the low-α disc. Given the distribution of this substructure in this chemical composition plane, and its very disc-like orbits, we suggest it is possible that Aleph is constituted by warped/flared disc populations. Because the data upon which our work and that by Naidu et al. (2020) are based come from different surveys, it is important however to ascertain that selection function differences between APOGEE and H3 are not responsible for our samples to have very different properties, even though they are selected adopting the same kinematic criteria. In an attempt to rule out that hypothesis we cross-matched the Naidu et al. (2020) catalogue with that of APOGEE DR17 to look for Aleph stars in common to the two surveys. We find only two such stars,⁶ which we highlight in Fig. 6 with purple edges. While this is a very small number, the two stars seem to be representative of the chemical composition of the APOGEE sample of Aleph stars, which is encouraging.

• LMS − 1: Our results from Fig. 6 show that the LMS-1 occupies a locus in the Mg–Fe plane which appears to form a single sequence with the GES at the lower metallicity end. Based on its Mg and Fe abundances, and the overlap in kinematic planes of LMS-1 and GES, we suggest it is possible that these two substructures could be linked, where LMS-1 constitutes the more metal-poor component of the GES. We investigate this possible association further in Section 5.

• Arjuna: This substructure occupies a distribution in the Mg–Fe plane that follows that of the GES. Despite the sample being lower in numbers than that for GES, we still find that across −1.5 < [Fe/H] < −0.8 this halo substructure overlaps in the Mg–Fe plane with that of the GES substructure. Given the strong overlap between these two systems, as well as their proximity in the E–L_z plane (see Fig. 4), we suggest it is possible that Arjuna could be part of the GES substructure, and further investigate this possible association in Section 5.

• I′itoi: Despite the small sample size, we find I’itoi presents high [Mg/Fe] and low [Fe/H] values (the latter by construction), and occupies a locus in the Mg–Fe plane that appears to follow a single sequence with the Sequoia (all three samples) and the GES sample.

• Nyx: The position of this substructure in the Mg–Fe strongly overlaps with that of the high-α disc. Given this result, and the disc-like orbits of stars comprising this substructure, we conjecture that Nyx is constituted by high-α disc populations, and further investigate this association in Section 5. Furthermore, in a similar fashion as done for the Aleph substructure, we highlight in Fig. 6 with purple edges those stars in APOGEE DR17 that are also contained in the Nyx sample from Necib et al. (2020), in order to ensure that our results are not biased by the APOGEE selection function. We find that the overlapping stars occupy a locus in this plane that overlaps with the Nyx sample determined in this study, and the high-α disc.

4.2 Iron-peak elements

Following our analysis of the various substructures in the Mg–Fe plane, we now focus on studying their distributions in chemical abundance planes that probe nucleosynthetic pathways contributed importantly by Type Ia supernovae. We focus on nickel (Ni), which is the Fe-peak element that is determined the most reliably by ASPCAP, besides Fe itself. For the distribution of the structures in other iron-peak element planes traced by ASPCAP (e.g. Mn, Co, and Cr), we refer the reader to Figs F1–F3 in Appendix  F.

The distributions of the halo substructures in the Ni–Fe plane are shown in Fig. 8. We find that the distributions of GES, Sgr dSph, the Helmi stream, Arjuna, LMS-1, and the three Sequoia samples occupy a locus in this plane that is characteristic of low mass satellite galaxies and/or accreted populations of the Milky Way, displaying lower [Ni/Fe] abundances than the low- and high-α disc populations (e.g. Shetrone et al. 2003; Mackereth et al. 2019; Horta et al. 2021; Shetrone et al., in preparation). In contrast, the data for Heracles and Thamnos display a slight correlation between [Ni/Fe] and [Fe/H], connecting with the high-α disc at [Fe/H] ∼ −1 (despite the differences of these substructures in the other chemical composition planes with in situ disc populations). Conversely, we find that the Aleph and Nyx structures clearly overlap with in situ disc populations at higher [Fe/H] values, agreeing with our result for these substructures on the Mg–Fe plane. The distribution of I’itoi shows a spread in [Ni/Fe] for a small range in [Fe/H], that is likely due to observational error at such low metallicities.

Interpretation of these results depends crucially on an understanding of the sources of nickel enrichment. Like other Fe-peak elements, nickel is contributed by a combination of SNIa and SNe II (e.g. Weinberg et al. 2019; Kobayashi, Karakas & Lugaro 2020). The disc populations display a bimodal distribution in Fig. 8, which is far less pronounced than in the case of Mg. This result suggests that the contribution by SNe II to nickel enrichment may be more important than previously thought (but see below). It is thus possible that the relatively low [Ni/Fe] observed in MW satellites and halo substructures has the same physical reason as their low [α/Fe] ratio, namely, a low star formation rate (e.g. Hasselquist et al. 2021). This hypothesis can be checked by examining the locus occupied by halo substructures in a chemical plane involving an Fe-peak element with a smaller contribution by SNe II, such as manganese (e.g. Kobayashi et al. 2020). If indeed the [Ni/Fe] depression is caused by a decreased contribution by SNe II, one would expect [Mn/Fe] to display a different behaviour. Fig. F1 confirms that expectation, with substructures falling on the same locus as disc populations on the Mn–Fe plane.

Another possible interpretation of the reduced [Ni/Fe] towards the low metallicity characteristic of halo substructures is a metallicity dependence of nickel yields (Weinberg et al. 2021). We may need to entertain this hypothesis since, in contrast to the results presented in Fig. 8, no [Ni/Fe] bimodality is present in the solar neighbourhood disc sample studied by Bensby et al. (2014), which may call into question our conclusion that SNe II contribute relevantly to nickel enrichment. It is not clear whether the apparent discrepancy between the data for nickel in Bensby et al. (2014) and this work is due to lower precision in the former, sample differences, or systematics in the APOGEE data.

Given the distribution of the substructures in the Ni–Fe plane, our results suggest that: (i) Sgr dSph, GES, Sequoia (all three samples), and the Helmi stream substructures show a slightly lower mean [Ni/Fe] than in situ populations at fixed [Fe/H], as expected for accreted populations in the Milky Way on the basis of previous work (e.g. Nissen & Schuster 1997; Shetrone et al. 2003); (ii) Heracles and Thamnos fall on the same locus on the Ni–Fe plane, presenting a slight correlation between [Ni/Fe] and [Fe/H]; (iii) as in the case of the Mg–Fe plane, Arjuna and LMS-1 occupy a similar locus in the Ni–Fe plane to that of GE/S, further supporting the suggestion that these substructures may be associated; (iv) Aleph/Nyx mimic the behaviour of in situ low-/high-α disc populations, respectively.

4.3 Odd-Z elements

Aside from α and iron-peak elements, other chemical abundances provided by ASPCAP/APOGEE that are interesting to study are the odd-Z elements. These elements have been shown in recent work to be depleted in satellite galaxies of the MW and accreted systems relative to populations formed in situ (e.g. Hawkins et al. 2015; Das et al. 2020; Hasselquist et al. 2021; Horta et al. 2021). For this paper, we primarily focus on the most reliable odd-Z element delivered by ASPCAP: aluminium. For the distribution of the structures in other odd-Z chemical abundance planes yielded by APOGEE (namely, Na, and K), we refer the reader to Figs G1 and G2 in Appendix  G.

Fig. 9 displays the distribution of the substructures and parent sample in the Al–Fe plane, using the same symbol convention as adopted in Fig. 6. We note that the parent sample shows a high density region at higher metallicities, displaying a bimodality at approximately [Fe/H] ∼ −0.5, where the high-/low-[Al/Fe] sequences correspond to the high-/low-α discs, respectively. In addition, there is a sizeable population of aluminium-poor stars with −0.5 ≲ [Al/Fe] ≲ 0 ranging from the most metal-poor stars in the sample all the way to [Fe/H] ∼ −0.5.⁷ This is the locus occupied by MW satellites and most accreted substructures, with the exception of Aleph and Nyx. Note also that the upper limit of the distribution of the Heracles population on this plane is determined by the definition of our sample (see Horta et al. 2021, for details).

The majority of the substructures studied occupy a similar locus in this plane, which agrees qualitatively with the region where the populations from MW satellites are usually found (e.g. Hasselquist et al. 2021). There is strong overlap between stars associated with the GES, the Helmi stream, Arjuna, Sequoia (all three samples), and LMS-1 substructures. More specifically, we find that GES dominates the parent population sample at [Fe/H] < −1, being located at approximately [Al/Fe] ∼ −0.3. At a slightly higher value of [Al/Fe] ∼ −0.15 and similar metallicities, we find Heracles and Thamnos. In contrast, Sgr dSph is characterized by an overall lower [Al/Fe] ∼ −0.5 value, which extends below the parent disc population towards higher [Fe/H], reaching almost solar metallicity. Within the −2 ≲ [Fe/H] ≲ −1 interval, Heracles, GES, Helmi streams, Thamnos, Nyx, and, to a lesser extent, Sequoia, show some degree of correlation between [Al/Fe] and [Fe/H]. Towards the metal-poor end, we find the LMS-1 located at [Al/Fe] ∼ −0.3, which is consistent with the value found for I’itoi, although the sample of aluminium abundances for this latter structure is very small and close to the detection limit. As in the case of magnesium and nickel, all three Sequoia samples occupy the same locus in the Al-Fe plane as Arjuna, which strongly overlap with GES. Again in the case of the Al-Fe plane, we find that the case for Nyx and Aleph follow closely the trend established by in situ disc populations.

4.4 Carbon and nitrogen

In this subsection, we examine the distribution of stars belonging to various substructures in the C–Fe and N–Fe abundance planes, shown in Figs 10 and  11, respectively. We note that in these chemical planes, we impose an additional surface gravity constraint of 1 < logg < 2 in order to minimize the effect of internal mixing along the giant branch.

In the C–Fe plane, most substructures are characterized by subsolar [C/Fe], displaying a clear correlation between that abundance ratio and metallicity. The exceptions, as in all previous cases, are Aleph and Nyx, which again follow the same trends as in situ populations. Interestingly, the Sgr dSph presents the lowest values of [C/Fe] at fixed [Fe/H], tracing a tight sequence at approximately [C/Fe] ∼ −0.5, spanning from −1.4 < [Fe/H] < −0.2, approximately ∼0.5 dex below that of the Galactic disc. In the case of I’itoi, due to the low numbers of stars in this sample we are unable to draw any conclusions.

The distribution of substructures in the N–Fe plane follows a different behaviour than seen in all other chemical planes. Again, except for Aleph and Nyx, all systems display a trend of increasing [N/Fe] towards lower metallicities, starting at [Fe/H] ≲ −1. This trend cannot be ascribed to systematics in the ASPCAP abundances or evolutionary effects, as the abundances are corrected for variations with log g. Nitrogen abundances are notoriously uncertain, particularly in the low metallicity regime. In their analysis of DR16 elemental abundances, Jönsson et al. (2020) report the presence of a trend of nitrogen abundances with T_eff. They suggest that such trend is likely caused by systematics in the comparison samples, which are based on rather uncertain analyses of optical spectra. They find that the relation between [C/N] and log g is in good agreement with theoretical expectations. The compilation by Kobayashi et al. (2020) shows that the [N/Fe] trend at low metallicity is strongly dependent on the analysis methods. Discerning the source of systematics in the ASPCAP abundances at [Fe/H] ≲ −1 is beyond the scope of this paper, which focuses on a strictly differential analysis of the data, within a metallicity regime where ASPCAP elemental abundances attain exceedingly high precision and accuracy (Section 5).

For completeness, data covering the whole range of log g for all substructures are displayed in the (C + N)–Fe plane in Fig. E1 in Appendix  E. By combining carbon and nitrogen abundances, we minimize the effect of CNO mixing along the giant branch. In this plane, MW satellites and accreted populations typically display a lower [(C + N)/Fe] chemical composition than their in situ counterparts (e.g. Hasselquist et al. 2021; Horta et al. 2021). This is in fact what we observe for all the structures identified, again with the exception of Aleph and Nyx, whose locus overlaps with that of in situ disc populations.

4.5 Cerium

Cerium is a neutron capture element of the s-process family, with a large enrichment contribution from AGB stars (Sneden, Cowan & Gallino 2008; Jönsson et al. 2020; Kobayashi et al. 2020). In Fig. 12, the disc sample at [Fe/H]>−1 has a roughly horizontal locus at [Ce/Fe]≈−0.1 dominated by stars in the high-α population and an upward-pointing triangular locus dominated by stars in the low-α sequence, reaching [C/Fe]≈+0.4 at [Fe/H]≈−0.2. The scatter within each of these components is large and may be partly observational. The presence of substantial Ce in high-α stars suggests that massive stars with short lifetimes make a significant, prompt contribution. The rising-then-falling trend in the low-α population is expected from the metallicity-dependent yield of intermediate mass AGB nucleosynthesis: at low [Fe/H] the number of seeds available for neutron capture increases with increasing metallicity, but at high [Fe/H], the number of neutrons per seed becomes to low to produce the heavier s-process elements (Gallino et al. 1998). See Weinberg et al. (2021) for plots of [Ce/Mg] versus [Mg/H] and further discussion of the disc trends.

In this chemical composition plane, we find that all the identified substructures, with the exception of Aleph and Nyx, present [Ce/Fe] abundances that follow the mean trend with [Fe/H] of the parent population until [Fe/H] ∼−1. Aleph and Nyx have higher [Fe/H] stars that generally lie within the broad disc locus. Interestingly, the Sgr stars with [Fe/H] >−1 show a rising [Ce/Fe] trend that tracks the behaviour of the low-α disc population. This trend is not obvious in the other substructures, though with the exception of Aleph and Nyx they have few stars at [Fe/H] >−1. We interpret this upturn in both Sgr and the low-α disc as the signature of an AGB contribution with a metallicity-dependent yield.

At [Fe/H] <−1, the parent halo population and most substructure stars exhibit mildly subsolar [Ce/Fe] with substantial scatter, which may have a significant observational component. The [Ce/Fe] is similar to that of typical high-α disc stars at [Fe/H] >−1. However, these stars have lower [α/Fe] than the high-α disc, the signature of Fe enrichment from Type Ia SNe, so if the Ce in these populations is a prompt contribution from massive stars one might have naively expected them to have depressed [Ce/Fe]. It is difficult to disentangle the effects of metallicity-dependent Ce yields, differences in the relative contributions of high-mass and intermediate-mass stars, and the impact of Type Ia SN enrichment on the Fe abundance; further observational investigation and theoretical modelling will be needed to do so. The [Ce/Fe] locus of substructure stars at −2<[Fe/H]<−1 is similar to that in the dwarf satellites studied by Hasselquist et al. (2021).

4.6 The [Al/Fe] versus [Mg/Mn] plane

Having studied the distribution of the identified structures in chemical abundance planes that aimed to give us an insight into the different nucleosynthetic pathways, contributed either by core-collapse, type Ia supernovae, and AGB stars, we now focus our attention on analysing the distribution of substructures in the stellar halo in an abundance plane that lends insights into the accreted or in situ nature of Galactic stellar populations: namely, the [Mg/Mn]–[Al/Fe] plane.

Fig. 13 shows the resulting distribution of the various structures in the [Mg/Mn]–[Al/Fe] plane. This chemical plane has been proposed by Das et al. (2020) as a means to distinguishing accreted populations from those formed in situ. Horta et al. (2021) showed that in situ stellar populations with a small degree of chemical evolution occupy the same locus in that plane as accreted populations. By construction, the Heracles substructure falls in the accreted locus of the diagram (see fig. 1 in Horta et al. (2021) for reference). However, our results show that all the other structures, except for Aleph and Nyx, also occupy the accreted locus of this plane. Interestingly, we find that although the GES, Sgr dSph, the Helmi stream, Sequoia (all three samples), Thamnos, LMS-1, Arjuna, and I’itoi substructures occupy the same locus, they appear to show some small differences. Specifically, we find that Sgr dSph occupies a locus in this plane positioned at lower mean [Mg/Mn] than the other structures. This is likely due to Sgr being more recently accreted by the Milky Way, and thus had more time to develop stellar populations with enriched Mn abundances that have been contributed on a longer time-scale by type Ia supernovae. This feature is also seen to a lesser extent for I’itoi. Conversely, at higher [Mg/Mn] values (but still low [Al/Fe]) we find GES, Heracles (by construction), Thamnos, LMS-1, the Helmi stream, Sequoia (all three samples), and Arjuna. The distribution of these substructures in this chemical plane reinforces the hypothesis of these halo substructures arising from an accreted origin.

In a similar fashion to the other chemical composition planes, we find that Aleph and Nyx overlap with in situ (disc) populations at higher [Al/Fe], suggestive that Aleph and Nyx are likely substructures comprised of in situ disc populations.

5 A QUANTITATIVE COMPARISON BETWEEN HALO SUBSTRUCTURE ABUNDANCES

After qualitatively examining the chemical compositions of the previously identified halo substructures in a range of chemical abundance planes, we now focus on comparing the abundances in a quantitative fashion using a χ² method. To do so, we compare the mean value of thirteen different elemental abundances, manufactured across different nucleosynthetic channels, for each substructure at a fixed metallicity that is well covered by the data. The set of elemental abundances chosen to run this quantitative test was determined based on the distribution of the parent sample in the respective chemical composition plane, where we only chose those elements that did not display a large scatter towards low metallicity due to increased abundance uncertainties (on the order of σ ∼0.15 dex). Out of the initial 20 elemental abundances available in ASPCAP (excluding Fe), we utilize the following thirteen elements: C, N, O, Mg, Al, Si, S, K, Ca, Ti, Mn, Ni, and Ce. We note that Na, P, V, Cr, and Co were removed due to the large scatter at low metallicity, whereas Cu and Nd were not considered due to ASPCAP not being able to determine abundances for these elements in APOGEE DR17.

We employ in our quantitative comparison a methodology that is immune to, or very minimally affected by, selection effects. For that reason, the metallicity distribution function does not enter our analysis. Our quantitative comparisons between various substructures are instead based on the abundance ratios at a reference [Fe/H], and the slope of the relation between abundance ratios and [Fe/H]. The method proceeds as follows:

• We select a high- and low-α (disc/halo) population based in Fig. 5 for reference, and utilize these as representative samples for any comparison between in situ populations and halo substructures. In order to account for any distance selection function effects, we restrict our high-/low-α samples to stars within d<2 kpc, and also determine an ‘inner high-α’ sample (restricted to R_GC<4 kpc) which we use to compare to Heracles (which has a spatial distribution that is largely contained within ∼4 kpc from the Galactic Centre).

• Before performing any chemical composition comparisons, we correct the abundances for systematic trends with surface gravity. Systematic abundance variations trends with log g can be caused, on one hand, by real physical chemical composition variations as a function of evolutionary stage and/or, on the other, by systematic errors in elemental abundances as a function of stellar parameters. The former chiefly impact elements such as C, N, and O, whose atmospheric abundances are altered by mixing during evolution along the giant branch. The latter impact various elements in distinct, though more subtle, ways (see discussion in Weinberg et al. 2021). Surface gravity distributions of various substructures differ in important ways (Fig. 1), so that such systematic abundance trends with log g can induce spurious artificial chemical composition differences between substructures. We thus follow a procedure similar to that outlined by Weinberg et al. (2021) to correct each elemental abundance using the full parent sample. As systematic trends with log g are more important towards the low and high ends of the log g distribution, we restrict our sample to stars within the 1<log g<2 range. We then fit a second order polynomial to the [X/H]–log g relation, and calculate the difference between that fit and the overall [X/H] median. The difference between these two quantities for any given log g is then added to the original [X/H] values so as to produce a flat relation between [X/H]_corrected and log g. We then use these values to determine corrected [X/Fe] abundances. In a recent study, Eilers et al. (2021) pointed out that simple corrections for abundance trends as a function of log g could erase real differences associated with abundance gradients within the Galaxy. That is because a magnitude limited survey may cause an artificial dependence of log g on distance. Our study aims at contrasting the chemical compositions of substructures that are in principle associated with spatially self-contained progenitors. Thus, systematic differences linked to spatial abundance variations within each structure are irrelevant for our purposes, so a straightforward correction for abundance variations with log g are perfectly acceptable for our goals.

• Upon obtaining corrected abundances for every halo substructure, we determined the uncertainties in the abundances using a bootstrapping resampling with replacement method (utilizing the astropy.stats.bootstrap routine by Price-Whelan et al. 2018). We generated 1000 realizations of the X–Fe chemical composition planes for every element and every halo substructure sample in order to assess the scatter in the abundance distribution. For example, we generated 1000 realizations of the C–Fe distribution for GES by drawing 2353 values from the observed distribution, with replacement (where 2353 is the size of our GES sample).

• For each one of the 1000 bootstrapped realizations of a chemical composition plane of a halo substructure, we determine the [X/Fe] value at a given metallicity ([Fe/H]_comp) that is covered by both halo substructures being compared by taking a 0.05 dex slice in [Fe/H] around [Fe/H]_comp and determining the median value for stars in that [Fe/H] interval. This yields 1000 [X/Fe] median values for each of the thirteen elements studied for every halo substructure (and disc sample) compared. We note that the [Fe/H]_comp values selected to perform this comparison are safely within the metallicity (and stellar parameter regime) where ASPCAP analysis is very reliable.

  We take the mean and standard deviation of the medians distribution for every [X/Fe] as our representative [X/Fe] and uncertainty value, respectively, and use these to quantitatively compare the chemical abundances between two populations. This sample median from the mean of the bootstrap medians is always close to the full sample median itself.

• Upon obtaining the mean and uncertainty chemical abundance values for every halo substructure at [Fe/H]_comp (for a range of thirteen reliable elemental abundances in APOGEE), we quantitatively compare the chemical compositions across halo substructures adopting a χ² statistic to assess the chemical similarities between different substructures. This quantity was computed by using the following relation:

  $$\begin{eqnarray} \chi ^{2} = \sum _{i} \frac{\Big (\mathrm{[X/Fe]}_{i,\mathrm{sub}} - \mathrm{[X/Fe]}_{i,\mathrm{ref}}\Big)^{2}}{\Big (\sigma _{\mathrm{[X/Fe]}_{i,\mathrm{sub}}}^{2} + \sigma _{\mathrm{[X/Fe]}_{i,\mathrm{ref}}}^{2}\Big)}, \end{eqnarray}$$

  (1)

  where [X/Fe]_sub and [X/Fe]_ref are the abundances of the halo substructure and the compared reference stellar population, respectively, and σ_{[X/Fe], sub} and σ_{[X/Fe], ref} are the corresponding uncertainties to those abundance values. Since GES is the halo substructure for which our sample is the largest, we use this substructure as our main reference against which all other substructures, as well as in situ disc populations are contrasted.

• Lastly, in order to infer if two stellar populations present consistent chemical abundances in a statistical manner, we determine the probability value of the χ² result for twelve degrees of freedom using the scipy (Virtanen et al. 2020) stats.chi2.cdf routine. In addition to the χ² value, we also compute a metric of separation (defined as Σ_[X/Fe]), that is calculated by setting the denominator of equation (1) equal to 1. This separation metric provides an additional way to quantify how similar the chemical compositions of two halo substructures are that is unaffected by the sample size (as smaller halo substructure samples will have larger uncertainties on their mean abundance values). This separation metric provides an additional way to quantify how similar the chemical compositions of two halo substructures are (see Fig. 14 for an example) that is unaffected by the sample size.

• For the case of Heracles and Aleph, as the selection of these substructures relies heavily on the use of [Al/Fe] and [Mg/Fe], respectively, we remove these elements when comparing these substructures, and reduce the numbers of degrees of freedom to eleven when calculating the χ² probability value.

In order to develop a more clear notion of the meaning of the resulting χ² values resulting from the above comparisons, we perform an additional exercise aimed at gauging the expected χ² values for the cases where two samples are identical to each other, or very different. To accomplish this, we draw, for each substructure, three N_sub-sized random samples, two from the high-α and one from the low-α disc samples, where N_sub is the size of the sample of that substructure. We then calculate, for each substructure, two χ² values, one resulting from the comparison of the high-α disc against itself, and the other from the comparison of the high-α disc against the low-α disc samples. As the high- and low-α disc both cover a similar range in [Fe/H], and their abundances vary with [Fe/H], we select a narrow bin in [Fe/H] from which to draw our high- and low-α disc samples (namely, between −0.45 < [Fe/H] < −0.35). This enables us to obtain random samples of high- and low-α disc populations at the same [Fe/H], and allows us to directly compare the mean and scatter values of the chemical abundances using the χ² method.

The reader may inspect the resulting χ² and associated probability (⁠|$p_{\chi ^{2}}$|⁠) values obtained in Table 2. Here, a |$p_{\chi ^{2}}$| ∼ 1 signifies that two populations are statistically equal, and |$p_{\chi ^{2}}$| ∼ 0 means that they are statistically different. For the quantitative comparisons, we employ a threshold of |$p_{\chi ^{2}}$|  = 0.1 as our benchmark, where we will deem two substructures to be statistically similar if their associated probability value is higher than |$p_{\chi ^{2}}$| > 0.1, and different if below |$p_{\chi ^{2}}$| < 0.1. We note that this definition is arbitrary, and chosen based on testing the similarity between two substructures in chemical space.

The resulting quantitative comparison between the halo substructures indicates that approximately half of the halo substructures are statistically equal with regards to their chemical abundances (see Figs 15−17). For example, the χ² comparison between the GES and two of the Sequoia samples implies that these three substructures are statistically the same, as we find that the M19 and N20 Sequoia samples all have a high probability of being statistically similar to the GES substructure (⁠|$p_{\chi ^{2}}$| > 0.4). Conversely, for the K19 Sequoia shows clear differences with the GES substructure, given its |$p_{\chi ^{2}}$| = 0 and χ² = 47.9 values. Moreover, we also find striking similarities between GES and other halo substructures. For example, when comparing LMS-1 to GES we obtain a probability value of 0.46. Similarly, an ever closer match is found for the GES-Arjuna comparison, yielding a probability value of 0.92. Along similar lines, we find that when comparing Nyx to the high-α disc, we obtain a probability value of 0.78, which confirms our initial hypothesis that the Nyx is not an accreted substructure, but instead is a stellar population constituted of high-α disc stars. Moreover, we find that when comparing Heracles, Sgr dSph, and Thamnos with the GES that the probability of these substructures being statistically equal is ∼0. This result is not entirely surprising, as all these substructures are postulated to be debris from separate accretion events, and thus their chemical compositions are likely to differ. Interestingly, we find two surprising results: (i) despite the Aleph substructure presenting qualitatively the same chemistry as the low-α disc, its χ² value yields a probability of 0.04, suggesting that these are not as similar as initially hypothesized. These differences are manifested most prominently in [Ti/Fe] and [Ce/Fe]; (ii) although Heracles occupies a position in several chemical composition planes that appears to follow a single sequence with the inner high-α populations, the χ² yielded by that comparison indicates that these are statistically different (with a probability value of |$p_{\chi ^{2}}$| = 0). For this comparison, we find that [O/Fe], [Mg/Fe], and [Si/Fe] (i.e. the lighter α elements) are the main culprits for this difference, and to a lesser degree [S/Fe].

For the complete resulting probability values obtained when comparing the chemistry between every halo substructure, as well as the high- and low-α disc, we refer the reader to Fig. 18.

6 DISCUSSION

6.1 Summary of substructure in the stellar halo

Having qualitatively and quantitatively compared the chemical abundances of all the halo substructures under study, in this Section we discuss the results obtained for each substructure in the context of previous work.

6.1.1 Heracles

Stars from this substructure follow low energy, often eccentric orbits with low L_z, being largely phase-mixed in velocity and action space. All these features are to be expected in the scenario where Heracles was a massive system that merged early in the history of the Milky Way. Under this hypothesis, both the fact that the young Milky Way system was smaller in size, and that dynamical friction would have driven the system quickly into low energy orbits, it would lead to this stellar population sinking it into the heart of the Galaxy (Horta et al. 2021; Pfeffer et al. 2021).

The chemical compositions of the stars associated with Heracles are in broad agreement with this scenario. The distribution of Heracles in the α-Fe plane does not display the α-Fe knee or shin components of chemically evolved systems (McWilliam 1997), which suggests that its star formation was quenched before core–collapse supernovae contributed significantly to chemical enrichment. In Horta et al. (2021), we checked that this result was not an effect of the chemical composition criteria adopted in the selection of Heracles stars by comparing them with a similarly selected GES sample, which did display a clear α-Fe knee signature. Furthermore, Heracles occupies a locus in different chemical planes that resembles that of low mass galaxies of the Milky Way and/or accreted populations. In particular, Horta et al. (2021) show that the stars associated with Heracles make up a clump in the [Mg/Mn]–[Al/Fe] plane in the region occupied by accreted and/or chemically unevolved populations, characterized by low [Al/Fe] and high [Mg/Mn].

Recent work by Lane, Bovy & Mackereth (2021) suggests that the concentration associated by Horta et al. (2021) with Heracles in E–L_z space could be an artefact of the APOGEE selection function, so the authors caution that the reality of this halo substructure should be further tested. As discussed by Horta et al. (2021), phase mixing makes it especially hard to discriminate accreted systems from their in situ counterparts co-located in the inner few kpc of Milky Way on the sole basis of kinematics. None the less, numerical simulations predict their existence (e.g. Fragkoudi et al. 2020; Kruijssen et al. 2020; Pfeffer et al. 2021; Horta et al., in preparation). Detailed chemistry is thus crucial to tease out the remnants of accreted systems from the maze of in situ populations overlapping in the inner few kpc of the Galaxy.

For that reason, we examine closely the comparison between the abundance patterns of Heracles data and the inner high-α in situ population at the same [Fe/H] (Fig. 17). Our χ² analysis shows that the two populations differ chemically with high statistical significance (χ² = 55.9, |$p_{\chi ^{2}}$| = 0, similar, for instance, to the value obtained when comparing the LMC to the SMC, see Table 2). To check whether this result is sensitive to the choice of [Fe/H]_comp, we reran the analysis adopting [Fe/H]_comp = –0.9 and [Fe/H]_comp = –0.95 obtaining χ² value of 36.28 and 43.40, respectively, which corresponds to a probability value of |$p_{\chi ^{2}}$| = 0 in both cases.

In order to further test this claim, and to ensure that our [Al/Fe] cut employed to select Heracles does not bias our sample and results, we repeat the statistical comparison between Heracles and the inner high-α stars by selecting both populations in the [Mg/Mn]–[Al/Fe] plane adopting various [Al/Fe] criteria (namely, by shifting the diagonal line selection restricting the [Al/Fe] values⁸). The motivation behind this is two-fold: (i) if the diagonal [Al/Fe] selection is the sole reason for the differences in the chemical abundance properties between Heracles and the in situ high-α stellar population, then varying the position of this line should affect the resulting χ² and associated |$p_{\chi ^{2}}$| at a significant level, and should lead to Heracles and the inner high-α disc becoming chemically indistinguishable when loosening the [Al/Fe] restriction; (ii) if Heracles is a distinct population from the high-α stars, this test should provide a deeper understanding of the level of contamination in our Heracles sample by in situ populations which might be caused by this selection in the [Mg/Mn]–[Al/Fe] plane.

In order to perform this test, we select Heracles and inner high-α populations in the [Mg/Mn]–[Al/Fe] plane by shifting the diagonal line by ±0.05 dex and ±0.1 dex, respectively. Upon repeating the quantitative analysis (at [Fe/H] = –1), we find that Heracles and the inner high-α stars are still statistically distinct, and that the resulting χ² and |$p_{\chi ^{2}}$| vary by small amounts. Specifically, the resulting χ² values we obtain when shifting the diagonal line in the [Mg/Mn]–[Al/Fe] plane by (+0.10,  +0.05, −0.05, −0.10) are χ² = (51.32, 58.47, 58.93, 59.43), leading to a |$p_{\chi ^{2}}$| = 0 in all cases. As expected, when increasing (decreasing) the sample of Heracles stars to include more (less) chemically evolved populations, the resulting χ² value decreases (increases) slightly, despite these systems still being statistically different. This is likely due to the fact that we are selecting more (fewer) high-α disc contaminants in our Heracles sample. While the overlap between Heracles and inner Galaxy high-α stars is small at [Fe/H] = –1, given our results it is evident that their chemical compositions are statistically, given our methodology, distinct.

In a recent study Myeong et al. (2022) state that the chemical properties of Heracles are consistent with the recently reported Aurora population (see also Rix et al. 2022). Aurora was originally proposed by Belokurov & Kravtsov (2022) as an in situ component that is kinematically hot, with an approximately isotropic velocity ellipsoid, and very modest net rotation (see also Conroy et al. 2022). In a follow up study, Myeong et al. (2022) adopt a Gaussian Mixture Model (GMM) approach to describe the distribution of solar neighbourhood populations in chemodynamical space, including elemental abundances from APOGEE DR17 and GALAH DR3 (Buder et al. 2021), as well as Gaia-based orbital energy. A key aspect of the Myeong et al. (2022) methodology is that it refrains from adopting any selection criteria in chemical space. As a result, one of the components emanating from their GMM fit is a highly eccentric (i.e. e >0.85) population of APOGEE-Gaia stars deemed to be associated with Aurora, which occupy a locus in chemical abundance space that is qualitatively similar to that inhabited by the inner Galaxy Heracles population.

We examine this claim by applying our statistical method to compare the chemical compositions of Heracles and Aurora, adopting the Myeong et al. (2022) sample for the latter population. To maximize the statistical robustness of the test, the comparison is performed at [Fe/H] = –1.2, where the samples for both substructures are the densest. As shown in Fig. 19, we find a resulting χ² = 47.69 and |$p_{\chi ^{2}}$| = 0, which is similar to the value we have obtained when performing the chemical comparison between Heracles and inner high-α in situ stars. Moreover, we repeat this test restricting the Heracles sample to the stars with e >0.85, for consistency with the Aurora definition by Myeong et al. (2022). Similarly, we obtain a χ² = 40.40 and |$p_{\chi ^{2}}$| = 0. Interestingly, the chemical compositions of Aurora stars and those of the inner Galaxy high-α in situ populations, are statistically similar, yielding χ² = 10.17 and |$p_{\chi ^{2}}$|=0.52 when compared at [Fe/H] = –1. These quantitative results support the notion that Heracles and Aurora are chemically distinct. They also attest for the abundance pattern homogeneity of in situ stellar populations in the inner Galaxy and solar neighbourhood.

All in all, the overarching conclusion that should be taken from all these additional tests is that Heracles is a stellar population that is distinct, given its chemical abundance values available in APOGEE and the methodology employed in this work, to neighbouring in situ populations at the same [Fe/H].

Examining more closely the contrast between the abundance patterns of Heracles and in situ populations, including Aurora and their inner high-α counterparts, we find that they differ in interesting ways. By far the elements displaying the largest differences are oxygen, magnesium, and silicon, whose abundances are lower in Heracles than in in situ populations. At face value, this difference implies less SNII enrichment in Heracles than in the inner high-α disc, possibly reflecting a lower star formation rate.

We stress that it is possible that the chemical properties ascribed to Heracles may be to some extent influenced by our selection method, which is partly based on chemistry. That selection, however, is far from arbitrary. It is rather informed by the fact that the distribution of inner Galaxy stellar populations form a clear clump in the accreted/chemically unevolved region of the [Mg/Mn]–[Al/Fe] plane (see fig. 1 of Horta et al. 2021). Our χ² analysis shows that Heracles presents an almost unique abundance pattern, differing in a (sometimes small, yet) statistically significant way from most stellar populations under study. For instance, we find that the abundance patterns of Heracles and GES are different (i.e. |$p_{\chi ^{2}}$|  = 0). This is also the case when comparing Heracles to the Sequoia (all three samples), Arjuna, Thamnos, and Nyx.

We conclude by stating that our data are consistent with Heracles being the remnant of a major building block of the halo that merged with the Milky Way in its early history. A possible scenario that could accommodate our quantitative analysis with the claims by Myeong et al. (2022) and Rix et al. (2022) is the following: as the local manifestation of a large in situ halo, Aurora extends all the way to the inner Galaxy, where it overlaps with Heracles, a likely more compact remnant of an early accretion event. The two systems originally have distinct chemical compositions, with Aurora being characterized by a higher [α/Fe] ratio, likely due to a higher early star formation rate. Our measurements most likely underestimate the real chemical differences between Heracles and the in situ halo. That is because our Heracles sample is likely contaminated by inner stellar halo counterparts at a level that may be as high as ∼40 per cent, according to the estimate by Horta et al. (2021). Further studies based on an expanded set of elemental abundances for a larger sample, as well as detailed modelling, based both on cosmological numerical simulations and standard chemical evolution prescriptions, are required to definitively establish the origin of Heracles and the nature of early in situ populations.

6.1.2 Gaia-Enceladus/Sausage

Since its discovery (Belokurov et al. 2018; Helmi et al. 2018), the Gaia-Enceladus/Sausage (GES) substructure has been extensively studied, both from an orbital and chemical compositions perspective (e.g. Hayes et al. 2018; Mackereth et al. 2019; Vincenzo et al. 2019; Koppelman et al. 2019a; Feuillet et al. 2020; Monty et al. 2020; Simpson et al. 2020; Aguado et al. 2021b; An & Beers 2021; Hasselquist et al. 2021; Horta et al. 2021; Buder et al. 2022; Carrillo et al. 2022). In this work, we have identified a large sample of GES stars, and have shown that stars belonging to this population are characterized by intermediate-to-high orbital energies and high eccentricity, displaying no significant systemic disc-like rotation.

The chemical compositions of the GES substructure are characterized by lower [α/Fe] at [Fe/H] ≳ −1.6, than high-α disc for most α elements (namely Mg, O, Si, Ca, S), in agreement with previous work (e.g. Hayes et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Mackereth et al. 2019; Horta et al. 2021; Buder et al. 2022), and displays an α-knee at [Fe/H] ∼ −1.1 (see Fig. 20). The stellar populations of GES are also characterized by lower Al, C, and Ni than in situ populations, resembling the abundance patterns of stars from satellites of the Milky Way (Hasselquist et al. 2021; Horta et al. 2021). They also occupy the accreted/unevolved region of the [Mg/Mn]–[Al/Fe] plane. In summary, the chemical compositions of GES confirm the results from previous studies, and reinforce the idea that this halo substructure is the remnant of an accreted satellite whose debris dominate the local/inner regions of the stellar halo.

6.1.3 Sequoia

The Sequoia substructure was initially discovered due to the highly unbound and retrograde orbits of its constituent stars and associated globular clusters (e.g. Barbá et al. 2019; Matsuno et al. 2019; Myeong et al. 2019), which made it easily distinguishable from in situ populations. Since its discovery, several groups have sought to identify Sequoia within data from different surveys. Based on an analysis of the Villalobos & Helmi (2008) N-body simulations, Koppelman et al. (2020) suggested that, rather than a separate system, Sequoia may constitute a fringe GES population comprised of stars in low eccentricity retrograde orbits left over after the merger of that system. A different scenario was proposed by Naidu et al. (2021) who suggest that Sequoia was instead a satellite of GES.

Whether Sequoia is the remnant of a distinct accretion event or a component of GES, or even its satellite is obviously a very important question, which in principle one should like to address by contrasting the chemical compositions of the two systems. In order to achieve that goal theoretical predictions for the differences in chemical composition under the various scenarios are required. Myeong et al. (2019) found that the peak [Fe/H] of their Sequoia sample is lower by ∼0.3 dex than that of their GES counterparts. They argue that the mass–metallicity relation at the relevant redshift implies a 1:10 mass ratio, which in turn is consistent with a Sequoia mass estimate based on the total mass in globular clusters presumably associated with the system. Koppelman et al. (2020), in turn, argue that the metallicity gradient expected for a galaxy with the pre-merger size and mass of GES is consistent with the mean metallicity estimates for GES and Sequoia. However, their metallicity gradient expectation is based on measurements made by Ho et al. (2015) on the gas component of nearby star forming galaxies with M_⋆ < 10^9.6M_⊙. It is unclear how accurately those are representative of typical metallicity gradients in dwarf galaxies around 10 Gyr ago.

As discussed in Section 2.1.4, our study of the Sequoia system is based on three independent definitions, derived from the Myeong et al. (2019, ‘M19’), Koppelman et al. (2019c, ‘K19’), and Naidu et al. (2020, ‘N20’) studies. Because the N20 sample is by definition limited to fairly metal-poor stars, it is not considered in our quantitative analysis. The mean metallicities we infer when adopting either the M19 or the K19 Sequoia samples are 〈[Fe/H]〉 = −1.41 and −1.53, respectively. In contrast, we obtain 〈[Fe/H]〉 = −1.18 for our GES sample, suggesting a difference of ∼0.2–0.3 dex, which is roughly consistent with the estimate by Myeong et al. (2019). However, as explained in Section 1, our sample is not corrected for selection effects, so we do not ascribe much weight to estimates of any moments of the MDFs of either system. We thus focus on the differences in abundance pattern, for which unfortunately there are no theoretical predictions that can be confronted by our data. Consequently, the scope of our discussion is limited to comparing the abundance patterns of the two systems, in hopes that the ramifications of that comparison for the nature of Sequoia will be explored in the future on the basis of a more solid theoretical framework.

Qualitatively, we find that the three Sequoia samples occupy a similar locus in all chemical composition planes, resembling that of low mass satellite galaxies of the MW and/or accreted populations, displaying low Mg, Al, C, and Ni, that overlap with the GES substructure. As in the case of the systems discussed above, Sequoia occupies the accreted/unevolved region of the [Mg/Mn]–[Al/Fe] plane, which is encouraging given that two of those samples were selected purely on the basis of orbital parameters. Perhaps not unexpectedly, we find that the N20 sample appears to be simply the metal-poor tail of the M19 and K19 samples.

When running a quantitative χ² comparison between the M19 and K19 Sequoia samples, we find that their chemical compositions are technically indistinguishable from each other, yielding a probability value of |$p_{\chi ^{2}}$| = 0.23. Therefore we conclude that the chemical composition we obtain for the Sequoia system at the reference metallicity of [Fe/H]_comp = −1.3 does not depend on the selection criterion adopted. Despite this similarity, a quantitative comparison between Sequoia and GES leads to conflicting results. When contrasting GES with the M19 sample, the probability obtained is |$p_{\chi ^{2}}$| = 0.64, suggesting that Sequoia and GES are chemically indistinguishable. Conversely, when comparing the K19 Sequoia sample with GES, the resulting probability value is |$p_{\chi ^{2}}$| = 0.0 (with a χ² = 47.9), which implies a statistically significant difference.

Our samples for GES and Sequoia cover in a statistically meaningful way a wide range of metallicities, so that they lend themselves nicely to a more detailed comparison of the distributions of those two systems in the α-Fe plane. A crucial diagnostic is the metallicity of the ‘α-knee’, [Fe/H]_knee (Section 4.1), which is strongly sensitive to the details of the star formation history of the system. In recent studies, Matsuno et al. (2019), Monty et al. (2020), and Aguado et al. (2021b) have obtained different values for [Fe/H]_knee, ranging from ∼−2.5 to ∼−1.5, depending on the sample and/or α element considered. We take advantage of our homogeneous sample in order to revisit this measurement. We perform a piece-wise linear fit to our GES and Sequoia samples in order to accurately determine the position of the [Fe/H]_knee in those two systems, in a manner similar to the approach followed by Mackereth et al. (2019).

The resulting fits (solid lines) and 1-σ dispersions (shaded regions) are displayed along with samples for GES (blue), Sequoia M19 (green) and K19 (red) in Fig. 20. The piece-wise function was determined using the PiecewiseLinFit function included as part of the pwlf package (Jekel & Venter 2019). The [Fe/H]_knee values for GES and the M19 Sequoia sample are within ∼0.02 dex from each other. Conversely, the fit to the K19 sample yields a much lower [Fe/H]_knee ≲ −2. Although this result is in relatively good agreement with some of the previous work (e.g. Matsuno et al. 2019; Monty et al. 2020), we deem it quite uncertain, because our sample has very few stars with [Fe/H]≲ −2, so that the slope of the ‘plateau’ part of the distribution is poorly constrained. On the other hand, the slope of the Sequoia sequence on the Mg–Fe plane according to the M19 sample is significantly different from that of the K19 sample, suggesting a genuine difference in the star formation history inferred from each sample.

In summary, the situation regarding the nature of the Sequoia system remains somewhat inconclusive. The fundamental question one would like to answer is whether Sequoia is the remnant of a distinct system, or a component of GES. While chemistry can provide clues, a definitive answer depends on knowledge of typical or expected pre-merger internal chemical composition variations of dwarf satellites of the Milky Way. To our knowledge, such information is not available either in the form of observational measurements or theoretical predictions. When turning to the lower level question of whether Sequoia and GES have similar chemical compositions, we find additional uncertainty, as the result depends on how the Sequoia system is defined. Out of three definitions adopted in this paper, one delivers a Sequoia sample that differs from GES both in terms of abundance pattern and [Fe/H]_knee (Koppelman et al. 2020).⁹ Conversely, the sample defined as in Myeong et al. (2019) results in a Sequoia system that is much more similar to GES. To reach resolution of this impasse one will need a bigger and better defined observed sample of the Sequoia system, and preferably one that is not affected importantly by selection effects.

6.1.4 Helmi stream

The Helmi stream is a halo substructure that appears to jut out of the Galactic disc. It is characterized by stars on highly perpendicular (i.e. high L_⊥ and v_z) and prograde orbits (Helmi et al. 1999; Koppelman et al. 2019b), which appear to form a pillar at high orbital energies in the prograde wing of the E–L_z plane (see Fig. 4). Despite this substructure being discovered decades ago, its chemical compositions have not been studied in great detail until recently, largely due to the difficulty of obtaining a high confidence sample with reliable chemical composition information. Such studies suggest that the chemical abundance ratios of the Helmi stream resemble that of the rest of halo stars (Roederer et al. 2010; Gull et al. 2021; Limberg et al. 2021; Nissen, Silva-Cabrera & Schuster 2021). However, lower abundance ratios in Na and α elements (Matsuno et al. 2022a), as well as heavier neutron capture elements (Aguado et al. 2021a), have been recorded, suggestive that the Helmi streams have a smaller contribution from massive star nucleosynthesis.

Our results on the chemical compositions of the Helmi stream imply that, unsurprisingly, this substructure presents chemistry that is typical of accreted populations and/or dwarf satellites (i.e. low Mg, Al, C, and Ni). We also find that, despite it being selected purely on a kinematic and position basis, it occupies the accreted/unevolved region of the [Mg/Mn]–[Al/Fe] which, combined with its orbital properties, confirms its accreted nature. Furthermore, we find that, when comparing this halo substructure with the others studied in this work, the Helmi stream differs statistically from all other halo substructures. In Fig. 20 we display the data for the Helmi stream alongside a piecewise polynomial fit performed in the same way as described for GES and the Sequoia samples. Interestingly, the best fit for the Helmi stream indicate the occurrence of an ‘inverted knee’, whereby the slope of the relation between Mg–Fe becomes less negative. As discussed above for the case of the Sgr dSph, this is the signature of a burst of star formation. This is a very interesting result that is in line with recent work by Ruiz-Lara et al. (2022), who looked at the colour–magnitude diagram of Helmi Stream star candidates in the Gaia DR3 data and found hints that this halo substructure experienced a star formation burst approximately ∼8 Gyr ago. We encourage this finding merits further investigation on the basis of a larger sample.

We note that in a recent paper Aguado et al. (2021a) examined the chemical compositions of the metal-poor part of the Helmi stream (which they call S2). Their sample is comprised of stars within −3 < [Fe/H] < −1.5, and thus only overlap with our Helmi stream sample in the metal-poor regime. Aguado et al. (2021a) find an α-knee for this substructure at [Fe/H] ∼ −2, which is different in both shape and value to the one shown in this work. We attribute this difference to the lack of Helmi stream stars in our sample below [Fe/H] < −2. Encouragingly, however, we find that at [Fe/H]  = −2 our Helmi stream sample has [Mg/Fe] ∼ 0.3 dex, which corroborates the findings from Aguado et al. (2021a).

6.1.5 Arjuna

The existence of this substructure was proposed by Naidu et al. (2020) as part of the H3 survey (Conroy et al. 2019). Naidu et al. (2020) show that the MDF of the retrograde component of the halo displays three peaks, which they ascribe to Arjuna (the most metal-rich), Sequoia, and I’itoi (the most metal-poor, see Fig. 3). Naidu et al. (2021) argue that Arjuna corresponds to the outer parts of the GES progenitor, which, according to their fiducial numerical simulation was stripped early in the accretion process, thus preserving the highly retrograde nature of the GES approaching orbit. In additional support to that proposition, Naidu et al. (2021) point out that the peak [Fe/H] and mean [α/Fe] of Arjuna are in excellent agreement with those of GES, which in turn should be consistent with a much lower metallicity gradient in the GES progenitor than suggested by Koppelman et al. (2020).

The detailed quantitative comparison of the chemical compositions of Arjuna and GES (Fig. 15) shows that the similarity of these two systems indeed encompasses a broader range of elemental abundances, leading to a |$p_{\chi ^{2}}$|  = 0.92. In addition, the distributions of Arjuna and GES stars in the α-Fe plane are also very similar, with the two values for [Fe/H]_knee agreeing within 0.06 dex (see Fig. 20).

Therefore, our results are at face value in agreement with the suggestion by Naidu et al. (2021) that the stars associated with the Arjuna substructure were originally part of GES. That association predicts a very low metallicity gradient for GES at the time of the merger with the Milky Way. Further theoretical and observational work is required to ascertain the reality of that prediction (see discussion in, e.g. Horta et al. 2021; Naidu et al. 2021).

6.1.6 I’itoi

Similarly to Arjuna, the I’itoi substructure is a high-energy retrograde substructure identified by Naidu et al. (2020). However, it is comprised by more metal-poor stars (see Fig. 3) than its high-energy retrograde counterparts, Arjuna and Sequoia. Naidu et al. (2021) propose that I’itoi was, in fact, a satellite of GES, based on its low metallicity and high energy retrograde orbit. Our detailed comparison of the chemistry of GES and I’itoi suggests that their abundance patterns are consistent with a |$p_{\chi ^{2}}$|  = 0.3 (Fig. 16). We point out, however, that this result is highly uncertain, given the relatively small size of our I’itoi sample and its low metallicity ([Fe/H]_comp = −2.1), which places its stars in a regime where ASPCAP abundances are relatively uncertain. The matter needs revisiting on the basis of more detailed chemical composition studies applied to a larger sample.

6.1.7 Thamnos

Initially conjectured by its discoverers to be the amalgamation of two smaller systems (Koppelman et al. 2019c), Thamnos is a substructure that occupies a locus in the retrograde wing of the velocity and IoM planes. It is comprised of stars with intermediate orbital energy (i.e. E∼–1.8 × 10⁵ km² s⁻²) and fairly eccentric orbits (e ∼ 0.5), that occupies a position at the foot of GES in the Toomre diagram.

As in the case of most orbital substructures in this study, we find that the locus occupied by Thamnos in the Ni–Fe, C–Fe, and [Mg/Mn]–[Al/Fe] chemical planes resembles that of low-mass satellite galaxies and accreted populations of the Milky Way. However, Thamnos distinguishes itself from other substructures by showing a relatively high [α/Fe] ratio, although not as high as Heracles. In fact, Thamnos does not match the abundance pattern of any other substructure in this study.

6.1.8 Aleph

Aleph was identified in a study of the stellar halo based on the H3 survey (Naidu et al. 2020). In this work, we identify Aleph members by selecting from our parent sample stars that satisfy the selection criteria outlined in Naidu et al. (2020). We have also searched for stars that are included in both the sample by Naidu et al. (2020) and in APOGEE DR17 (highlighted in the chemical abundance figures with purple edges). As described in their work, this substructure is comprised of stars on very strongly prograde orbits with low eccentricity, which appear to follow the same distribution as the Galactic disc component, although at higher J_z values.

We find that the locus occupied by Aleph in various chemical planes is placed somewhere in between low- and high-α disc stars, while sitting squarely within the in situ region of the [Mg/Mn]–[Al/Fe] plane. This is also the case for the two stars common with the sample by Naidu et al. (2020). Our quantitative comparison of the Aleph chemistry with that of the low-α disc suggests a statistically significant difference, albeit on the borderline (χ²  = 19.8 and |$p_{\chi ^{2}}$|  = 0.05). This is in fact not surprising, seeing as the distribution of Aleph stars on the Mg–Fe plane straddles both the high and low-α discs (Fig. 6). These results suggest that Aleph may be a stellar population comprised of a mix of warped/flared low-α disc, and high-α disc stars, which also explains its location within the locus of in situ stellar populations in the [Mg/Mn]–[Al/Fe] plane.

In addition, the chemical abundances of Aleph are in good qualitative agreement with the chemical abundances of the Anticenter Stream (ACS) (Grillmair 2006), as reported by Laporte et al. (2020). Laporte et al. (2020) argue that the ACS is the remains of a tidal tail of the Galactic disc excited during the first pericentric passage of the Sagittarius dwarf galaxy. If this scenario is correct, the fact that Aleph presents chemical compositions that are barely distinguishable from the low-α disc suggests that the low-α disc of the Galaxy was in place sometime before this satellite interaction. We encourage future work to test this hypothesis.

6.1.9 LMS-1

LMS-1 is a metal-poor substructure comprised of stars on mildly prograde orbits at intermediate/high orbital energies. Yuan et al. (2020) identified this substructure by applying a clustering algorithm to SDSS and LAMOST data in the E–L_z plane. Although it presents great overlap with the GES in IoM space, it is suggested to be an independent substructure based on the detection of a metallicity peak in the MDF of the stars included within the E–L_z box defining this system (Naidu et al. 2020). Yuan et al. (2020) also note that there are potentially several globular clusters with similar metallicity and orbital properties. In this work, we identified LMS-1 candidates adopting the same selection criteria as Naidu et al. (2020)’s, obtaining a relatively small sample of only 31 stars.

We examine the distribution of LMS-1 stars in various chemical planes, concluding that its chemistry is consistent with those of other accreted systems, with all its stars falling in the ‘accreted/unevolved’ region of the [Mg/Mn]–[Al/Fe] plane. Furthermore, the ∼0.5 probability value obtained when comparing this halo substructure with GES suggests that at face value LMS-1 could be part of either of the omnipresent GES substructure. However, these comparisons are made at [Fe/H]_comp = −2.1, where our samples are small and ASPCAP abundances are relatively more uncertain and in situ and accreted structures tend to converge towards the same locus in the regions of chemical space sampled by APOGEE (e.g. Horta et al. 2021). Moreover, given its location in IoM space, it is possible that our LMS-1 sample is importantly contaminated by GES stars.

As mentioned in Section 6.1.4, Jean-Baptiste et al. (2017) show that a single accretion event can lead to multiple overdensities in orbital space (see also Koppelman et al. 2020). Due to the chemical similarity between LMS-1 and GES, and the close proximity between these two halo substructures in orbital planes (see Fig. 4), an association between these two systems seems tempting. However, we defer a firmer conclusion to a future when more elemental abundances are obtained for a larger sample of both LMS-1 and GES candidates.

6.1.10 Nyx

The Nyx substructure is conjectured to be a stellar stream in the solar vicinity of the Milky Way (Necib et al. 2020). Given the chemical compositions obtained for this substructure in this work and its strong overlap with the high-α disc in all the chemical planes studied, we suggest that the Nyx is likely comprised by in situ high-α disc populations. Our quantitative estimate of the similarity between Nyx and the stars from the high-α disc yields a very low χ² with associated likelihood probability of 0.75. We note that the stars in common between our sample and that of Necib et al. (2020) seem to boldly confirm this result. Along these lines, we note that our result is in agreement with a recent study focused on studying the chemical compositions of the Nyx substructure (Zucker et al. 2021), who also conjecture Nyx to be comprised of Galactic (high-α) disc populations.

7 CONCLUSIONS

The unequivocal association with accretion events of halo substructures identified on the basis of orbital information alone is extremely difficult, as demonstrated by recent numerical simulations (e.g. Jean-Baptiste et al. 2017; Koppelman et al. 2020; Naidu et al. 2021). Substructure in integrals of motion space can also be influenced or even artificially created by survey selection effects (Lane et al. 2021). Because the chemical compositions of halo substructures contain a fossilized record of the evolutionary histories of their parent galaxies, abundance pattern information can help linking substructure in orbital space to their progenitor systems.

In this work, we have utilized a cross-matched catalogue of the latest APOGEE (DR17) and Gaia (EDR3) data releases in order to study the chemodynamic properties of substructures previously identified in the stellar halo of the Milky Way. We have successfully distinguished stars in the APOGEE DR17 catalogue that are likely associated with the following substructures: Gaia-Enceladus/Sausage, Sagittarius dSph, Heracles, Helmi stream, Sequoia, Thamnos, Aleph, LMS-1, Arjuna, I’itoi, Nyx, Icarus, and Pontus. Using the wealth of chemical composition information provided by APOGEE, we have examined the distributions of the stellar populations associated with these substructures in a range of abundance planes, probing different nucleosynthetic channels. We performed a quantitative comparison of the abundance patterns of all the substructures studied in order to evaluate their mutual associations, or lack thereof. Below we summarize our main conclusions:

• We show that the chemical compositions of the majority of the halo substructures so far identified in the Galactic stellar halo (namely, Gaia-Enceladus/Sausage, Heracles, Sagittarius dSph, Helmi stream, Sequoia, Thamnos, LMS-1, Arjuna, I’itoi) present chemical compositions which resemble those of low-mass satellites of the MW and/or accreted populations. There are however a couple of exceptions, namely Nyx and Aleph, that do not follow this pattern and instead present chemical properties similar to those of populations formed in situ. Furthermore, in Appendix  A we discuss the nature of Icarus, concluding that this substructure is likely comprised of stars formed in situ, for which the ASPCAP abundances are unreliable.

• The chemical properties of the inner-Galaxy structure, Heracles, differ from those of its in situ metal-poor counterparts, both in the inner halo and in the solar neighbourhood (Aurora). In particular, Heracles displays lower abundances of α-elements oxygen, magnesium, and silicon than the in situ populations. A possible interpretation of this result is that the star formation rate of Heracles was lower than that of the early in situ halo. By the same token, Heracles is found to have higher [α/Fe] ratios than GES which, as suggested by Horta et al. (2021), is an indication of an early truncation of star formation and the resulting absence of an α-knee in the former system. The abundance pattern of Heracles is indeed found to differ from all of the other substructures studied in this work. We propose a scenario according to which Heracles and Aurora overlap spatially in the innermost Galaxy. Further studies based on additional elemental abundances for larger samples, as well as detailed numerical modelling, are required to disentangle the in situ and accreted metal-poor populations cohabiting the innermost Galactic halo (i.e. R_GC < 4 kpc).

• We show that a large fraction of the substructures studied (namely, (M19/N20) Sequoia, Arjuna, LMS-1, I’itoi) present chemistry indistinguishable from that of the omnipresent Gaia-Enceladus/Sausage. These findings bring into question the independence of these substructures, which are at least partially overlapping with GES in kinematic planes (see Fig. 4). In view of these similarities, claims in the literature about the nature of Sequoia as being originally a higher angular momentum component located in the outskirts of GES (Koppelman et al. 2020), a satellite of GES (Naidu et al. 2021), or an altogether unrelated system (Myeong et al. 2019) may need to be reexamined. The possibility that Sequoia (selected as done in Myeong et al. (2019) and Naidu et al. (2020)) was a detached, but much less massive galaxy than GES is likely challenged by their chemical similarity. On the other hand, confirmation that it, or Arjuna, might be the remains of populations originally located in the outskirts of GES depends crucially on the magnitude of chemical composition gradients one should expect for dwarf galaxies at z ≈ 2, and on whether that is a sufficiently discriminating criterion.

• We found that the halo hosts substructures which differ from GES in a statistically significant way. Three among those substructures display chemistry that is genuinely suggestive of an accreted nature, namely Heracles, Thamnos, and the Helmi stream (although for the latter this conclusion is not firm due to uncertainties in the chemistry and small sample size).

• Conversely, the chemistry of the two remaining substructures, Nyx and Aleph, is indistinguishable from that of in situ populations. We conjecture that Nyx forms part of the high-α disc. For the case of Aleph, we suggest that it could likely comprised of stars both from the low-α (flared/warped) disc as well as stars from the high-α disc, or be part of the Anti-Centre stream.

• The situation regarding the nature of the Sequoia system remains somewhat inconclusive given our results. Out of three definitions adopted in this paper, one delivers a Sequoia sample that differs from GES both in terms of abundance pattern and [Fe/H]_knee (Koppelman et al. 2019c). Conversely, the sample defined as in Myeong et al. (2019) results in a Sequoia system that is much more similar to GES. For the case of the Sequoia as defined in Naidu et al. (2020), the hard [Fe/H] cut prevents us from performing any quantitative comparison.

• Our results suggest that the local/inner (R_GC ≲ 20 kpc) halo is comprised of the debris from at least three massive accretion events (Heracles, GES, and Sagittarius) and two/three lower mass debris (Thamnos, Helmi stream, Sequoia). Upcoming large spectroscopic surveys probing deeper into the outer regions of the stellar halo (beyond R_GC ∼20 kpc) will likely uncover the debris from predicted lower mass/more recent accretions (e.g. van den Bosch et al. 2016; Horta et al., in preparation), and will enable the full characterization of those already known. Conversely, the precise contribution of massive accretion to the stellar populations content of the inner few kpc of the halo is still to be fully gauged. Heracles is likely the result of a real building block of the halo, likely the most massive (in relative terms) to have ever happened in the history of the Milky Way, but we are still scratching its surface.

This paper presents a chemical composition study of substructures identified (primarily) on the basis of phase-space and orbital information in the stellar halo of the Milky Way. Current and upcoming surveys will continue to map the stellar halo and will provide further clues to the nature and reality of phase-space substructures discovered in recent years. For the inner regions of the Galaxy and the local halo, the Milky Way Mapper (Kollmeier et al. 2017) and the Galactic component of the MOONS GTO survey (Gonzalez et al. 2020) will provide revolutionizing chemodynamical information for massive samples. For the outer regions of the stellar halo, the WEAVE (Dalton et al. 2012), 4MOST (de Jong et al. 2019), and DESI (DESI Collaboration 2016) surveys will provide spectroscopic data for millions of stars in both high and low resolution. In addition, in this work we have only studied the chemistry of the most massive accretion events. However, there is a plethora of lower mass halo substructures in the form of coherent stellar streams that has not been studied in this work and also require to be fully examined (see Li et al. 2021 for a recent example). The advent of surveys like S⁵ (Li et al. 2019) will aid in such endeavours. All this information, when coupled with the exquisite astrometry and upcoming spectroscopic information from the Gaia satellite, will provide the necessary information for further discoveries and examinations of substructure in the stellar halo of the Milky Way.

ACKNOWLEDGEMENTS

The authors thank GyuChul Myeong and Vasily Belokurov for making the IDs of the Aurora stars available, Vasily Belokurov for helpful comments on earlier versions of the manuscript, Rohan Naidu and Charlie Conroy for making available the H3 catalogue in digital format, and Paola Re Fiorentin and Alessandro Spagna for sharing their Icarus star’s APOGEE IDs. DH thanks Tadafumi Matsuno, Melissa Ness, Holger Baumgardt, and Cullan Howlett for helpful scientific discussions, and Sue, Alex, and Debra for their constant support. DG gratefully acknowledges financial support from the Dirección de Investigación y Desarrollo de la Universidad de La Serena through the Programa de Incentivo a la Investigación de Académicos (PIA-DIDULS). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of china, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

Software:astropy (Astropy Collaboration 2013; Price-Whelan et al. 2018), scipy (Virtanen et al. 2020), numpy (Oliphant 2006), matplotlib (Hunter 2007), galpy (Bovy 2015; Mackereth & Bovy 2018), topcat (Taylor 2005).

Facilities: Sloan Foundation 2.5 m Telescope of Apache Point Observatory (APOGEE-North), Irénée du Pont 2.5 m Telescope of Las Campanas Observatory (APOGEE-South), Gaia satellite/European Space Agency (Gaia).

DATA AVAILABILITY

All APOGEE DR17 data used in this study is publicly available and can be found at: https://www.sdss.org/dr17/

Footnotes

REFERENCES

APPENDIX A: THE REALITY OF ICARUS

The Icarus substructure is a metal-poor, [Mg/Fe]-poor structure that was identified by Re Fiorentin et al. (2021) using the APOGEE DR16 and GALAH DR3 data. We have attempted to identify Icarus star members by employing the criteria listed by Re Fiorentin et al. (2021), and identified only one potential member within the parent sample described in Section 2.1.9. We thus decided to retrieve the DR17 properties of the 41 Icarus stars presented by Re Fiorentin et al. (2021) on the basis of APOGEE DR16 data.

We found that only 3 out of the 41 Icarus candidates have T_eff and log g characteristic of red giant stars. The remaining 38 stars have stellar parameters consistent with a main sequence or subgiant nature (Fig. A1), which is the reason why the selection criteria described in Section 2 excluded them from our parent sample. This finding is confirmed by the positions of the Icarus member sample in the Gaia colour–magnitude diagram shown in Fig. A2. Icarus was conjectured to be an accreted substructure in the Galactic disc based largely on chemical abundance information. However, elemental abundances for subgiant and main-sequence, and particularly for M dwarfs, were not very reliable in APOGEE DR16 (although see Birky et al. 2020; Souto et al. 2020). Therefore, we decided to perform some checks to see if the resulting abundances for these (primarily) M dwarf stars are reliable in APOGEE DR17.

We checked the STAR₋FLAG and ASPCAP₋FLAG flags of the 41 Icarus stars determined by Re Fiorentin et al. (2021). The results from this examination showed that: (i) 25/41 stars had a ‘BRIGHT₋NEIGHBOUR’, ‘PERSIST₋HIGH’, ‘PERSIST₋LOW’, ‘LOW₋SNR’, SUSPECT₋BROAD₋LINES, SUSPECT₋ROTATION’, MULTIPLE₋SUSPECT’, flag set in their STAR₋FLAG, suggesting the occurrence of various issues with the observed spectra of these stars, which can lead up to uncertainties in stellar parameters and elemental abundances; (ii) 40/41 stars had one or more of the following ASPCAP_FLAGs raised: ‘STAR₋WARN’, ‘COLOUR₋TEMPERATURE₋WARN’, ‘SN₋WARN’, ‘VSINI₋WARN’, ‘VMICRO₋WARN’, ‘N₋M₋WARN’, ‘TEFF₋WARN’, ‘LOGG₋WARN’, suggesting that the results from the ASPCAP analysis for these stars are likely uncertain.

In light of these results, we chose to exclude Icarus from this study, and suggest that this substructure is not the debris from a cannibalized dwarf galaxy, but is rather likely comprised of (primarily) M-dwarf stars with unreliable APOGEE abundances in the disc of the Milky Way.

APPENDIX B: THE CHEMICAL COMPOSITIONS OF PONTUS ACCORDING TO APOGEE

As we have only been able to identify two star members belonging to the Pontus halo substructure, we will refrain from making any strong statement about its chemistry. However, based solely on the chemical compositions of the two member stars we have identified in this work, we find that the Pontus substructure presents elevated abundances in α-elements (namely, in O, Mg, Si), low C, Al, Ti, and Ce, and approximately solar Mn and Ni. The two stars we find associated with the Pontus substructure display similar chemical compositions to that of other halo substructures and/or satellite galaxies of the Milky Way. Although it is impossible to attribute any similarities or differences on the basis of two stars alone, given the close proximity between these two Pontus stars and Thamnos in the E–L_z plane, it is tempting to associate the former to the later. However, a larger sample of Pontus stars with reliable chemistry is required to ascertain this speculation.

APPENDIX C: CHEMICAL AND KINEMATIC PROPERTIES OF HALO SUBSTRUCTURES

 

  

 

APPENDIX D: α-ELEMENTS

 

 

 

 

APPENDIX E: (C + N)/Fe

 

 

 

APPENDIX F: IRON-PEAK ELEMENTS

 

 

 

APPENDIX G: ODD-Z ELEMENTS

 

 

 

APPENDIX H: [Fe/H] USED TO COMPARE SUBSTRUCTURES