© T. Ruiz-Lara et al. 2022

1. Introduction

Satellite accretion and galaxy interactions appear to have been key actors in the formation and evolution of galaxies. Although mergers must have happened throughout the history of the Universe (Ibata et al. 1994; Knapen et al. 2015; Paudel et al. 2018), at high redshift they were likely much more ubiquitous because of the higher average density of the Universe (Abraham et al. 1996; Conselice 2007). This is naturally incorporated into the current cosmological paradigm, which predicts that galaxies grow in mass hierarchically through the accretion of smaller systems (e.g. Qu et al. 2017). For many systems, and particularly for large disc galaxies, their accretion history is expected to have been violent in the early stages, and more quiescent at later times (see e.g. Blumenthal et al. 1984; Conselice 2014; Grand et al. 2017).

Reconstruction of the detailed accretion histories of external galaxies is unfortunately beyond current observational capabilities (although see e.g. D’Souza & Bell 2018; Mackey et al. 2019; Zhu et al. 2020; Davison et al. 2021). However, our privileged position within our own Galaxy enables the characterisation of its accretion history. This is because individual stars retain information of their formation sites in their dynamics and chemical abundances (see Freeman & Bland-Hawthorn 2002; Helmi 2020 and references therein). Debris from ancient accretion events is expected to have been deposited mostly in the Galactic halo (Searle & Zinn 1978; Bullock & Johnston 2005; Matsuno et al. 2019), with the most massive mergers dominating the inner regions (Helmi et al. 2002; Cooper et al. 2010).

With the goal of disentangling accretion events in the Galactic halo, analyses have focused on the space defined by integrals of motion (IoM) such as energy or angular momentum (Helmi & de Zeeuw 2000) or actions (McMillan & Binney 2008), based on the fact that stars from accreted systems remain clumped in such a space even after they have phase-mixed (Helmi & de Zeeuw 2000). The advent of full phase-space information for large samples of stars provided by the Gaia mission (Gaia Collaboration 2016, 2018a, 2021), sometimes in combination with ground-based spectroscopic surveys (e.g. LAMOST, Wang et al. 2020; Liu et al. 2019; RAVE, Steinmetz et al. 2020a,b; GALAH, Buder et al. 2021; or APOGEE Ahumada et al. 2020), has recently made such studies feasible in the halo near the Sun.

The emerging view from these studies (see e.g. Helmi et al. 2018; Massari et al. 2019; Koppelman et al. 2019a,b; Myeong et al. 2019; Belokurov et al. 2020; Yuan et al. 2020a) is that the stellar halo indeed contains large amounts of debris associated to different accretion events experienced by the Milky Way. Probably one of the most important events in the early history of the Galaxy is the merger with Gaia-Enceladus (GE, with a stellar mass of ∼6 × 10⁸ M_⊙ at the time of merging, ∼10 Gyr ago, Helmi et al. 2018, see also Belokurov et al. 2018). Stars from GE together with an in situ, more metal-rich population that was likely heated up from an ancient disc by this event, dominate the inner halo near the Sun (Gaia Collaboration 2018b; Koppelman et al. 2018; Gallart et al. 2019; Belokurov et al. 2020). Apart from GE, several other Galactic building blocks have been uncovered in the inner halo. The Helmi streams discovered by Helmi et al. (1999) are debris from a system with approximately 10⁸ M_⊙ in stars that was likely accreted to the Milky Way around 5 to 8 Gyr ago (Kepley et al. 2007; Koppelman et al. 2019b). More recently, Myeong et al. (2019) reported the existence of another accreted system, Sequoia. With a stellar mass of 5 × 10⁷ M_⊙, Sequoia may have provided the Milky Way with most of the nearby retrograde halo stars, including several globular clusters (Myeong et al. 2018). Thamnos, another possible building block of the Milky Way halo and composed of stars with high binding energy, relatively low inclination, and retrograde orbits, was first reported in Koppelman et al. (2019a). Thamnos was originally divided into two components (Thamnos 1 and Thamnos 2), which could have originated from the same progenitor, a stellar system with M_⋆ ≤ 5 × 10⁶ M_⊙. In addition to these events, a growing number of studies have unveiled several small and subtle stellar substructures (Wukong, Arjuna, Rg5, etc.) possibly linked to other minor merger events in the history of our Galaxy (Myeong et al. 2018; Naidu et al. 2020; Yuan et al. 2020a).

Besides studies making use of the distribution of stars in IoM space, the characterisation of the age–metallicity relation of Galactic globular clusters as well as their dynamical and chemical properties led (Kruijssen et al. 2019a) to postulate the presence of a possible major building block in the inner Galaxy¹ (see also Massari et al. 2019; Kruijssen et al. 2020; Horta et al. 2020, 2021). Comparisons with the E-MOSAICS simulations (Pfeffer et al. 2018; Kruijssen et al. 2019b) suggest that this ‘Kraken’-like event might have had a stellar mass of around 2 × 10⁸ M_⊙ at the time of the merger (similar to that of GE), around 11 Gyr ago (Kruijssen et al. 2020).

Therefore, through the application of ever more sophisticated techniques to data sets of increasing quality, we have gained a better and more complete picture of the accretion history of the Milky Way. Nevertheless, some important questions remain, such as how many independent merger events our Galaxy has experienced; or what their characteristics were. By answering these seemingly simple questions, not only will we have a better understanding of the role of accretion in galaxy formation (Cooper et al. 2013; Qu et al. 2017; Pillepich et al. 2018), but the characterisation of the different building blocks of our Galaxy will provide key information on galaxy formation and evolution in different environments and at different times. For example, characterisation of the smallest (accreted) systems would allow us to further constrain the subhalo mass fraction, and therefore cosmological models (e.g. Madau et al. 2008; Wetzel et al. 2009; Geen et al. 2013; Bozek et al. 2013; Fielder et al. 2019).

To reach these goals, statistically robust identification of the various substructures that may be associated to merger debris is necessary, as well as reliable determination of the likelihood of stars to belong to any of the substructures. In this work, we use the outcome of a data-driven hierarchical clustering algorithm (see Lövdal et al. 2022, Paper I from now on) to this end. In Paper I, we identified 67 highly significant clusters in IoM space for nearby stars with halo-like kinematics. Through inspection of the properties of the hierarchical tree of the clusters (via a dendrogram), and comparison to previous identifications of substructures in the nearby halo, we tentatively suggested that it can be divided in six main groups and a number of independent clusters. However, it is only through proper characterisation of the internal properties of all the significant clusters detected by the algorithm, together with a thorough comparison of their properties, that we can establish their independence or relationship in case of a common origin. This is the goal of the present paper.

The organisation of this paper is as follows. In Sect. 2, we present the sample of stars used for our analyses, and briefly summarise the main characteristics of our clustering approach. Section 3 presents the methodology used to characterise and compare the stellar content of clusters, which is based on isochrone fitting to describe the distribution of stars in the colour–absolute magnitude diagram (CaMD), and on their metallicity distribution functions. These comparisons allow us to define possible independent substructures, which we describe in Sect. 4. We then proceed to characterise their average age and chemical properties and provide an updated view of the Galactic halo in Sect. 5. An interesting result from our analysis is that we find that the building blocks contributing stars with very retrograde orbits are not only more metal poor but are also on average older than those in prograde orbits. This and other conclusions are provided in Sect. 6.

Please note that in what follows, we mostly use the word cluster to refer to the different groups of stars identified by the clustering algorithm of Paper I, that is, to clustered data in dynamical space, rather than to gravitationally bound systems, unless stated otherwise. We use the words group or substructure to refer to the different associations of clusters found that could potentially correspond to different events in the history of our Galaxy (i.e. different building blocks).

2. Data and clustering procedure

This work is based on the results presented in Paper I, where we introduced a data-driven clustering algorithm to systematically identify kinematic structures in the stellar halo near the Sun. For completeness, we provide a short summary below and refer the reader to Paper I for further technical details.

We use a sample of stars from the local halo (within 2.5 kpc from the Sun), defining halo stars according to |V − V_LSR|> 210 km s⁻¹, with V_VLSR = 232 km s⁻¹ (McMillan 2017), where V is the total velocity vector of a star. The sample is constructed from the Gaia Early Data Release 3 (EDR3, Gaia Collaboration 2021), and complemented where possible with line-of-sight velocities (vlos) from GALAH DR3 (Buder et al. 2021), APOGEE DR16 (Ahumada et al. 2020), RAVE DR6 (Steinmetz et al. 2020a,b), and the LAMOST DR6 medium (Liu et al. 2019) and low resolution (Wang et al. 2020) surveys. We require precise parallaxes (parallax_over_error > 5), low renormalised unit weight error (ruwe < 1.4), and low vlos uncertainty (below 20 km s⁻¹). This results in a final sample with 51 671 stars.

To this data set we applied (in Paper I) the single linkage algorithm in three-dimensional (IoM) space, including total energy (E) and the angular momentum in the z-direction (L_z) and in its perpendicular component (L_⊥). We refer the reader to Sect. 3.2 of Paper I for the exact definitions and computation of these features. The single linkage algorithm defines a series of connected components (or potential clusters) in the data set, as each step of the algorithm forms a new group by connecting the two closest data points not yet in the same group. We then evaluate the significance of each potential cluster produced in this process by comparing the density in an ellipsoidal region around the potential cluster with the density in the same region in an artificially generated halo. This procedure returns a set of (potentially hierarchically overlapping) statistically significant clusters, and final (non-overlapping) clusters are extracted by selecting the subsets with the highest statistical significance. We then refine those clusters by computing Mahalanobis distances (D_ij) for each star x_j with respect to the centre of every cluster C_i (assuming that this can be fitted with an ellipsoid, i.e. a multivariate Gaussian in 3D characterised by a covariance matrix Σ_i). The Mahalanobis distance² between star x_j and cluster C_i is defined as

$$\begin{array}{c}\hfill D_{\mathit{ij}}=\sqrt{{({\mathit{\mu }}_i-{\mathit{\mu }}_j)}^T{\mathrm{\Sigma }}_i^{-1}({\mathit{\mu }}_i-{\mathit{\mu }}_j)},\end{array}$$(1)

where μ_j denotes the location of star x_j and μ_i that of cluster C_i (mean of the cluster distribution to be precise). In what follows, we consider that a star x_j belongs to cluster C_i (core member) if their Mahalanobis distance is below a certain limit (D_ij <  2.13, which corresponds to the 80th percentile of the distribution of Mahalanobis distances for the stars originally in any cluster; see Paper I for more information).

The above procedure identifies 67 clusters in our data set with a significance of at least 3σ above that expected for a randomised halo distribution³. Table A.1 provides an overview of the main properties of these refined clusters. Around 13.8% of the stars in this sample (7127 stars) are within a Mahalanobis distance of 2.13 from any of these clusters⁴. To establish possible relations between clusters, we compute the Mahalanobis distance between them using the above definition (Eq. (1)), but where the covariance matrix is Σ_ij = Σ_Ci + Σ_Cj, with Σ_Ci and Σ_Cj describing the covariance matrices of the two cluster distributions and μ_i and μ_j the mean location of clusters C_i and C_j, respectively.

Figure 1 (adapted from Paper I) shows the relationship between the 67 extracted clusters in the form of a dendrogram based on the Mahalanobis distances between the clusters in IoM space. This figure clearly shows that the clusters identified are not fully independent, as they are grouped together in larger structures that even in some cases present a finer hierarchical substructure. Some of these (marked in colour) large structures correspond to previously studied Galactic halo building blocks as reported in Paper I. In that paper, we tentatively drew the limit for clusters to be linked at a Mahalanobis distance of ∼4.0. In Fig. 2 we show the distribution in IoM of such structures including those clusters that are kept as separate entities (according to our tentative cut). Using this analysis as a starting point, in this paper, we further investigate the independence of the substructures by comparing their stellar population characteristics.

Fig. 1. Relationship between the significant clusters according to the single linkage algorithm, depicted according to their Mahalanobis distance in the IoM space (Paper I). The different colours represent a tentative classification into different substructures, namely six large groups (named A to F) and a set of individual clusters (in grey; see text for details). The order of the clusters along the x-axis is such that it facilitates the interpretation and readability of the dendrogram, and does not reflect internal properties. In our analyses (e.g. tables, etc.), we broadly follow this same order.

Fig. 2. Distribution of stars for the tentative substructures in IoM space: E vs. Lz (top two rows) and L⊥ vs. Lz (bottom two rows). Each panel plots one of the six large groups and the independent clusters from Paper I. Colours are as in Fig. 1 except for the independent clusters panel. The full original sample upon which this study is built is in grey in all panels.

3. Methodology

In this section, we explain the main tools used to fully assess whether two clusters or groups of clusters – identified as described above – can be linked together in a substructure (i.e. they have a common physical origin), or whether they are fully independent. We make use of Gaia EDR3 astrometric and photometric data as well as information provided by the aforementioned ground-based spectroscopic surveys.

3.1. Isochrone fitting

The precision of the photometric data and distances to individual stars provided by Gaia allows the construction of accurate CaMDs down to the main sequence turn-off and below. Although issues related to completeness and uncharacterised selection effects prevent us from easily inferring full star formation histories from CaMD fitting techniques (Cignoni & Tosi 2010; Gallart et al. 2005, 2019), the quality of these observed CaMDs enables the estimation of an average age and metallicity via comparison with stellar isochrones. To this end, we use the updated BaSTI isochrones (Hidalgo et al. 2018) in the Gaia EDR3 photometric system. In particular, we employ the α-enhanced version of the BaSTI models ([α/Fe] = 0.4), with the canonical He abundance of Y = 0.247⁵.

The observed colour–magnitude diagrams from Gaia are transposed to the absolute magnitude plane by making use of distances computed by inverting the observed parallax after correcting for an offset of 0.017 mas (Lindegren et al. 2021). Although it is possible to improve the characterisation of this parallax offset in Gaia EDR3 (see e.g. Lindegren et al. 2021; Ren et al. 2021; Groenewegen 2021), this turns out to have a negligible effect on our analyses. We also correct for the effect of reddening on the magnitudes of individual stars using the dust map from Lallement et al. (2018) and the recipes to transform to Gaia magnitudes provided in Gaia Collaboration (2018b) ⁶.

Furthermore, we define the observed CaMD by considering only stars fulfilling the extra condition that the BP/RP excess factor (phot_bp_rp_excess_factor) verifies:

$$0.001+0.039\times \mathtt{bp\_rp}<log(\mathtt{phot\_bp\_rp\_excess\_factor})$$

and

$$log(\mathtt{phot\_bp\_rp\_excess\_factor})<0.12+0.039\times \mathtt{bp\_rp}$$

where bp_rp is the apparent (G_BP − G_RP) colour. We also avoid stars affected by high extinction, A_G >  0.7, as their position in the CaMD will be less robust, as well as horizontal branch stars and blue stragglers (that occupy the region to the left of the dashed grey line indicated in the right-hand panel of Fig. 3).

Fig. 3. Isochrone fitting procedure applied to clusters #60 and #61 together, which are associated with the Helmi streams. Left: age–metallicity plane, where the background colours show the distribution of the d2 parameter used for assessing the quality of the fit obtained. The average age and metallicity, together with the corresponding errors, are represented by the black point with errorbars. To compute these average values, we use isochrones with a d2 parameter within its lowest 10th percentile (small black dots). Horizontal dotted lines show the [Fe/H] values adopted for the isochrones. The histogram of d2 values is displayed in the bottom left for all the different isochrones considered. Right: CaMD for stars in the union of clusters #60 and #61. The best-fit isochrone is depicted in red, and all compatible isochrones (black dots in the left-hand panel) are overplotted in black with transparency. Stars bluer than the grey dashed lines (horizontal branch and blue stragglers, shown with transparency) are avoided in the fit. Every star is colour-coded according to its weight in the fit (see text for details).

To determine the average age and metallicity (hereafter [Fe/H]) of any group of stars, be it individual clusters or several clusters linked together as a group or substructure, we developed an isochrone-fitting procedure which is carried out in a series of steps. The set of isochrones considered in the fitting procedure cover the full range of ages (200 Myr to 14 Gyr, in steps of 0.2 Gyr) but we restrict the range in [Fe/H] to values according to our knowledge from spectroscopic surveys (LAMOST LRS), considering only isochrones corresponding to metallicities from ⟨[Fe/H]⟩−σ([Fe/H])_{LAMOST, LRS} to ⟨[Fe/H]⟩+σ([Fe/H])_{LAMOST, LRS} (with ⟨[Fe/H]⟩ being the mean and σ([Fe/H])_{LAMOST, LRS} the standard deviation of the sample). This range encompasses on average 70% of the stars with [Fe/H] determinations in each of the substructures under analysis⁷; see Sect. 4. This constraint is only applied if at least ten stars in the cluster, group, or substructure have [Fe/H] determinations from LAMOST LRS; otherwise the full set of isochrones is used (ranging from [Fe/H] = −2.5 to [Fe/H] = 0.06⁸).

We then measure how well a given isochrone reproduces the distribution of stars in the CaMD with the following parameter:

$$\begin{array}{c}\hfill d_2=\sqrt{\frac{\sum _{i=1}^n|X_i-X_{i,\mathrm{iso}}{|}^2\times W_i}{\sum _{i=1}^n{W}_i}},\end{array}$$(2)

where |X_i − X_i, iso|² represents the distance from the CaMD position (X) occupied by star i to the closest point of the isochrone (X_i, iso) and W_i is the weight given to this particular star, defined as the inverse of the sum in quadrature of the photometric errors (in colour and absolute magnitude).

To provide an average value and an uncertainty for the age and [Fe/H] of a given structure, we consider all possible isochrones compatible with the observed distribution of stars. To this end, we select those within the lowest 10th percentile of the d₂-distribution. The best estimate of the age and [Fe/H] (and its uncertainty) is then computed using a weighted average (and respectively, dispersion) where the weights are given by the inverse of the d₂ values.

We now discuss one specific case as an example of how the outlined procedure works. In Paper I, we identified a link between clusters #60 and #61 (substructure E in Fig. 1), namely that they are part of the Helmi streams. Figure 3 shows the distribution of the d₂ parameter in the age–metallicity plane sampled by the BaSTI isochrones (left), as well as the CaMD of the 143 core member stars from clusters #60 and #61 fulfilling our selection criteria. From this isochrone fit, we find that the Helmi stream stars are on average 12.1 ± 1.1 Gyr old and have an average metallicity of [Fe/H] = −1.45 ± 0.11 (while the average spectroscopic [Fe/H] of the 100 stars with LAMOST LRS spectra is −1.6 with a dispersion of 0.3). These values are in very good agreement with previous studies (e.g. Koppelman et al. 2019b).

We note that we expect most of the clusters identified in Paper I to be complex systems, composed of multiple populations. In this sense, the isochrone fitting approach must be understood as an attempt to compute an average age of the system, which will typically not provide a perfect fit to the CaMD.

3.2. How to assess the independence of halo substructures

The clustering procedure presented in detail in Paper I identified 67 significant clusters in a 3D IoM space. As mentioned above, on the basis of the dendrogram shown in Fig. 1 and a comparison to previous work, we established in that paper that these 67 clusters could tentatively be grouped into six larger substructures (named A to F) and several individual clusters (#1, #2, #3, #4, and #68).

Figure 2 shows the distribution in the E − L_z plane (top two rows) and in L_⊥ − L_z (bottom two rows) of the stars belonging to these tentative structures. We can roughly identify substructure A as part of the hot thick disc; associate B with Thamnos 1+2 (Koppelman et al. 2019a); C with GE (Helmi et al. 2018, see also Belokurov et al. 2018); D with the region reported to be occupied by Sequoia (Myeong et al. 2019); E as part of the Helmi streams (Helmi et al. 1999); and find that F is located in the region of the E − L_z plane that is expected to be populated by debris from the Sagittarius dwarf galaxy (Sgr, Naidu et al. 2020) and that is approximately occupied by globular clusters linked with Sgr in Massari et al. (2019).

Apart from this rough identification, we note a large amount of complexity in the dendrogram of Fig. 1. For example, substructure A presents two groups of clusters (A1 and A2) linked at small Mahalanobis distances. Substructure C (GE) can be clearly split into several subgroups (C1 to C6). In both cases, several other clusters are linked at larger Mahalanobis distances. Other examples are substructures B (Thamnos 1+2) and D (Sequoia region), which show two core members with a third cluster at a larger distance.

To fully assess if the substructures with similar dynamical properties have a common origin or whether they can be split further, we developed a thorough, exhaustive, and quantitative approach that makes use of their internal stellar population properties. We compare the metallicity distribution function (MDF) of their stars and their colour distribution in the CaMD via Kolmogorov–Smirnov (KS) statistical tests. For the MDF analysis, we concentrate on stars with determinations of [Fe/H] from LAMOST LRS (the spectroscopic survey with the best coverage of our sample, Wang et al. 2020), thereby avoiding heterogeneity from surveys that may use different metallicity scales. For the colour distribution analysis, first we find the isochrone that best represents the group or structure that we are studying (understood as the union of the potential candidate clusters). Then, we compute the colour distance from this isochrone to the individual stars making up each cluster. Specifically, we proceed as follows:

– Step 1: We identify the smallest units or subgroups within each structure, i.e. composed by clusters linking at the smallest Mahalanobis distances (e.g. A1, A2, C1, ..., C6, B, D, E, F, see Fig. 1).

– Step 2: For each (sub)group, we compare the MDFs and colour distributions of each of its constituent clusters with that of the rest of the (sub)group. For example, for the A1 group, composed of clusters #34, #28, and #32, we make three comparisons: #34 with (#28 ∪ #32); #28 with (#32 ∪ #34); and #32 with (#28 ∪ #34). We decide whether two distributions (ϕ₁, ϕ₂) are compatible if the computed KS p-value P_ϕ1, ϕ2 >  0.05 when comparing their MDFs. When comparing the colour distributions, we lower the threshold to 0.03 to claim that two distributions are different (because of the lower discerning power of colours).

– We note that for clusters with less than 20 stars, Poisson noise might affect the results of the KS tests. In such cases, we first compute P_ϕ1, ϕ2 as above. To make the comparison meaningful between the two structures C₁ and C₂ where N₂ <  20, we randomly draw N₂ stars from C₁ (${\varphi }_1^{\mathrm{\prime }}$). We then run a second KS test comparing ${\varphi }_1^{\mathrm{\prime }}$ with its parent distribution, ϕ₁, that is, we determine $P_{{\varphi }_1,{\varphi }_1^{\mathrm{\prime }}}$ to assess the effect of small number statistics. This procedure is repeated 100 times. If P_ϕ1, ϕ2 <  0.05, and the recovered $P_{{\varphi }_1,{\varphi }_1^{\mathrm{\prime }}}>0.05$ for at least 80% of the cases, we cannot consider random noise to be responsible for the low measured p-value, and therefore we can conclude that ϕ₁ and ϕ₂ are indeed drawn from different distributions.

– Step 3: Clusters kept as independent because of their larger linking Mahalanobis distances, or those discarded from the groups studied in Step 2, are evaluated to assess their relation with the coherent groups identified, whilst keeping an eye on the distribution of their stars in IoM. No clusters will be added to any group if they are not clearly related in IoM. This step allows us to refine the groups defined in Step 1.

– Step 4: All groups and independent clusters (those that remain as such after Step 3) are compared via KS tests again.

As an additional consistency check and whenever possible, we make use of chemical abundance ratios from APOGEE DR16 if there are enough stars in each cluster or substructure with abundance measurements to define a clear trend after removing stars that have STAR_BAD flag in ASPCAPFLAG, FE_H_FLAG ≠ 0, or MG_FE_FLAG ≠ 0. In this study, we focus on the behaviour of [Mg/Fe] as a function of [Fe/H].

As in the previous section, we use the Helmi streams again to exemplify the procedure outlined above. The algorithm presented in Paper I associates⁹ 95 and 48 stars with clusters #60 and #61, respectively. Of these, 66 and 34 stars respectively have [Fe/H] determinations from LAMOST LRS. The right panel of Fig. 4 shows that they present similar metallicity distributions, displaying a p-value of 0.99 according to the KS test when compared to one another. The left panel of Fig. 4 shows their distribution in the CaMD colour-coded by the distance in colour to the best-fitting isochrone for the whole structure (when the clusters are merged). The distribution of distances from the best-fitting isochrone is shown in the central panel, and is consistent between the clusters, presenting a p-value of 0.6. Thus, as already mentioned in Paper I, we can conclude that both clusters are part of the same structure. The origin of the splitting of the Helmi streams into two significant clusters is probably the presence of a resonance induced by the Galactic potential (see Dodd et al. 2022).

Fig. 4. Analysis of the stellar populations of the Helmi streams. Left: distribution in the CaMD of the stars classified as part of the Helmi streams (clusters #60 and #61) colour-coded according to the colour distance to the best-fitting isochrone (black dashed line, see Sect. 3.1). Center: histograms of the difference in colour for stars in clusters #60 and #61, and showing compatible distributions with a KS p-value of ∼0.6. Right: metallicity distribution function from LAMOST LRS for both clusters (p-value of 0.99). Error bars were computed assuming Poisson statistics.

4. Analysis: identification of independent halo substructures

To obtain a more precise and updated view of the history of the Galactic inner halo, in this section we determine which substructures can be considered independent according to all the information at our disposal and using the methodology described above. In Sect. 5 we use this information to characterise their chemical and stellar population properties and interpret our findings.

4.1. Substructure A

Figure 1 shows that substructure A can first be split into two groups A1 and A2 with clusters (#34, #28, #32) and (#31, #15, #16), respectively, linking at small Mahalanobis distances. A further eight clusters are linked at larger Mahalanobis distances. We now proceed to apply the methodology described in Sect. 3.2.

The results of Steps 1 and 2 (Sect. 3.2) applied to groups A1 and A2 are given in Table 1. According to this table, cluster #28 from A1 should be considered independent, because the KS test p-value for the comparison of its MDF to that of the union of clusters #34 and #32 (in A1) is ∼0.003. We checked that this value is not affected by Poisson noise given that #28 has only 19 stars with [Fe/H] determinations from LAMOST LRS. We find such a low value in only 3% of the random realisations. The comparison of the colour distributions also yields a relatively low p-value (0.19). Based on this, and the information in Table 1, we conclude that A1 = (#34 ∪ #32) and A2 = (#31 ∪ #15 ∪ #16)¹⁰.

We now proceed to Step 3 of the method and consider the independent clusters #33, #12, #40, #42, #46, #38, #18, and #30, as well as cluster #28, and whether they should be associated to A1 or A2. The results are given in Table 2, where we see that, in addition to cluster #28, now clusters #38, and #40 also present significant differences in their MDF and should be considered independently. Cluster #12 is a special case, as it displays clear differences in the colour distributions; we return to this cluster below.

The remaining clusters are compatible with both A1 and A2. Indeed, if we compare these two groups (Step 4 of our method), we see that they should form a single substructure, which we refer to as substructure A.

Using all this information, Fig. 5 summarises the final characterisation of this substructure. This figure shows the distribution in IoM, as well as the MDFs of substructure A and of clusters #12, #38, #28, and #40 separately. We can see that these clusters present distributions in IoM and/or MDF that clearly differ from that of A. For example, the peak of the MDF of cluster #38 is shifted to lower metallicities (with this cluster having higher L_⊥ than typical of A), while the MDFs of clusters #28 and #40 are somewhat narrow (peaking at higher and lower metallicities with respect to A, respectively). On the other hand, cluster #12 displays a clear peak at low metallicity. Such stars are responsible for the low KS p-value reported above, and are apparent in Fig. 6 given that they display bluer colours in the CaMD, particularly on the main sequence (represented in the figure as crosses). It seems plausible that cluster #12 suffers from a large amount of contamination from substructure A and that one should consider only the subset of stars associated to the metal-poor peak as truly independent. We return to this possibility later on.

Fig. 5. Characterisation of substructure A and the independent clusters #12, #38, #28, and #40. Left: distribution in the L⊥ − Lz plane. The original sample of halo stars is shown as a grey background. Middle: normalised MDFs for the three largest (in number of stars) structures/clusters. As an inset we add the MDF of cluster #12 alone to highlight its apparently double population. See text for details. Right: normalised MDFs for clusters #28 and #40, showing a narrow MDF. We add the MDF of substructure A (dashed, blue line) for comparison. Error bars in all plots are computed assuming Poisson statistics.

Fig. 6. Comparison of the distribution of stars in CaMD for cluster #12 (blue) and substructure A (red). The black line represent the isochrone that best reproduces the distribution of stars in the CaMD for substructure A (see Sect. 3.1). The inset panel shows the distributions of colours with respect to the best fit isochrone for both structures. We mark low-metallicity stars ([Fe/H] < −1.4) with blue crosses.

In summary, substructure A is the union of (#15, #16, #46, #31, #32, #30, #42, #18, #33, #34), with clusters #12, #38, #28, and #40 to be considered separately. This classification is supported by the final comparison given in Table 3. We note that substructure A, which lies in the hot thick disc, does not overlap substantially with the previously reported Nyx overdensity – that has prograde thick disc-like kinematics (Necib et al. 2020a,b) – because of the dynamical selection of our sample of halo stars.

The low number of stars populating clusters #28 and #40 (29 and 10 stars, respectively) prevents a more thorough characterisation, and therefore we not study them in more detail in the following sections. We deem their origin to be unclear.

4.2. Substructure B: Thamnos 1 and Thamnos 2

Substructure B (overlapping with Thamnos1 and 2 in IoM) shows a configuration in which two core clusters (#8, #11) are linked at a larger Mahalanobis distance to a third one (#37). Table 4 shows the p-values of KS tests comparing the MDFs of these three clusters, and confirms that there indeed seem to be two distinct subgroups. Clusters 8 and 11 display the more similar MDFs, which are clearly different from that of cluster #37, although the distributions of stars in the CaMD are compatible. Consequently, we divide substructure B into B1 (cluster #37) and B2 (#8 ∪ #11).

These substructures are in qualitative agreement with Thamnos 1 and Thamnos 2 respectively, as defined in Koppelman et al. (2019a), although their mean E and L_z are different from those originally proposed. This likely is due to improvements in the astrometry in Gaia EDR3 compared to Gaia DR2 which moves the original selection of Thamnos stars to lower energies and less retrograde orbits. As in Koppelman et al. (2019a), B2 (Thamnos 2) shows an overall higher metallicity than B1 (Thamnos 1).

Figure 7 shows that there may be some amount of contamination in B1 from the more metal-rich stars in B2. It is likely this contamination that blurs the differences when comparing the colour distributions of these structures.

Fig. 7. Characterisation of the substructures identified within the tentative structure B. Left: distribution in the E − Lz plane. The original sample of halo stars is shown as a grey background. Right: normalised MDFs. We use blue for B1 (Thamnos 1) and cyan for B2 (Thamnos 2; see text for details). Error bars are computed assuming Poisson statistics.

4.3. Substructure C: Gaia-Enceladus

Gaia-Enceladus (GE) is considered the most prominent merger that the Milky Way experienced in its early history (Helmi et al. 2018). As such, it is no surprise that our clustering algorithm appears to have linked together nearly 40 individual clusters as a single substructure (tentative substructure C) in the region of IoM space thought to be occupied by GE. As in the case of substructure A, our clustering algorithm reveals a rich fine structure within C as shown in Fig. 1, with six groups composed of clusters linked at small Mahalanobis distances¹¹. Twelve clusters are tentatively marked as independent¹² (#6, #52, #14, #20, #19, #17, #43, #39, #36, #23, #5, and #21).

Tables 5 and 6 study the similarities in MDF and CaMD, respectively, between the six groups in substructure C. Based on these tables, we see that the six subgroups are indeed composed of the clusters with a probable common origin. The only notable exceptions being cluster #55 in group C3 and group C6 on the basis of their MDFs. Therefore, in Step 3, where we evaluate the independent or more isolated clusters, we consider the clusters originally in C6 (#9, #24, #35) and #55 separately.

This is exactly what is shown in Table 7, where we compare all these clusters with groups C1 to C5. Cluster #6 is particularly striking as it has an MDF that is clearly different from that of groups C1 to C5 (with the only exception being C2). In the remaining cases, visual inspection of their position in IoM and similarities between their MDFs and CaMDs allow us to assign them to groups C1 to C5, establishing the following subgroup configuration for substructure C:

– C1: Clusters #56 and #59,

– C2: Clusters #10, #7, #13, #20, #19, and #43,

– C3: #57, #58, #54, #50, #53, #47, #44, #51, #39, #35, #9, #24, #5, #23, #36, #21, and #17,

– C4: #48, #45, #25, #29, #41, #49, and #55,

– C5: #27, #22, and #26.

Apart from cluster #6 (as commented above), the more isolated cluster #52 is too far in the IoM space from any of these five subgroups in C (despite the similarity in MDFs and CaMDs). This is shown in Fig. 8, which plots their distribution in the E − L_z and L_⊥ − L_z planes. Cluster #14 on the other hand is close to subgroup C2 with both depicting the more retrograde motions of substructure C, but their MDFs are not consistent. Their position in E − L_z space roughly coincides with that of the Rg5 structure identified in Myeong et al. (2019; see also Paper I, and Yuan et al. 2020a; notably Naidu et al. 2020 associate this structure to Thamnos). However, the L_⊥ found by these latter authors appears quite different (< 500 kpc km s⁻¹ compared to > 1000 kpc km s⁻¹ of Rg5 stars in Yuan et al. 2020a). Therefore, either C2 or #14 is not related at all to Rg5, or they all stem from a fairly large progenitor (in which case, we would have expected it to also extend further in E − L_z). As we see below, all evidence suggests that both C2 and cluster #14 are part of GE.

Fig. 8. E − Lz (top) and L⊥ − Lz (bottom) distribution of the subgroups making up GE (substructure C). As in Fig. 2, the panels present each of the groups identified and the three independent clusters. The full original sample upon which this study is built is represented as a grey background in all panels. The full extension of substructure C is highlighted in silver.

Now we compare (pairwise) the MDFs of all the subgroups and independent clusters, again using KS statistical tests (our Step 4 of Sect. 3.2). As can be seen from Table 8, all groups and independent clusters can be considered compatible with each other. The only cluster that remains independent is #6 (with 122 stars, 88 with [Fe/H] determinations from LAMOST LRS), whose MDF is clearly different from those of the others¹³. This analysis suggests that substructure C, which unites all of the groups, can be firmly associated with GE.

Abundance determinations from APOGEE are available for a large number of the stars belonging to C/GE according to our algorithm. This allows us to study their [Mg/Fe] abundance ratios as a function of [Fe/H] as extra proof of their common origin. As shown in Fig. 9, all of the subgroups seem to follow the same chemical abundance track, which is compatible with that known for GE (Helmi et al. 2018; Mackereth et al. 2019). Therefore, the chemical abundance ratios also support our analysis and interpretation. Interestingly, Fig. 9 shows that the ⟨[Mg/Fe]⟩ at low [Fe/H] is lower than that of other stars in our sample at similar metallicity (see also Sect. 5.3).

Fig. 9. [Mg/Fe] abundance ratios as a function of [Fe/H] for substructure C (GE) from APOGEE DR16. We use different symbols for the different clusters. Colours for the groups (C1 to C5) follow Fig. 8. The small dots show abundances of all the halo stars within 2.5 kpc in APOGEE DR16.

Now that our clustering algorithm has clearly identified GE, in Fig. 10 we compare the MDF of all the stars that we have just established belong to GE with the overall MDF of our halo sample. This halo sample is composed of both accreted and in situ stars as discussed in previous works (e.g. Carollo et al. 2007, and especially the blue and red inner halo sequences reported in Gaia Collaboration 2018b; Haywood et al. 2018; Gallart et al. 2019; Di Matteo et al. 2019). The presence of two main components is clearly reflected in the bimodal shape of the halo MDF (shaded area). It is reassuring that the MDF of the stars identified as GE corresponds to only the more metal-poor component, which provides further support to the good performance of our approach in identifying accreted structures.

Fig. 10. Combined metallicity distribution function of all clusters in substructure C/Gaia-Enceladus. We overlay the [Fe/H] distribution for all 34 532 stars with [Fe/H] from LAMOST LRS in our initial halo sample (containing 51 671 stars). We note that our selection of GE stars based on our clustering algorithm naturally does not suffer from contamination from the hot thick disc (red sequence). Error bars were computed assuming Poisson statistics.

4.4. Substructure D: not only Sequoia

Substructure D contains some of the most retrograde stars in our halo sample and, based on the location in IoM, may be associated (at least in part) with Sequoia (Myeong et al. 2019). A quick visual inspection of Fig. 2 reveals three clusters in this region: two large ones with different E, L_z, and L_⊥ (#63 and #64, the latter being the most retrograde) and a third one (#62) with low L_⊥ but similar L_z and energy to those of the larger cluster #63 but easily distinguishable in the L_⊥ − L_z plane in Fig. 2.

Even though Table 9 suggests that the three clusters can be linked (p-values above 0.05, with cluster #62 being the most different according to their MDFs), their position in IoM and velocity space suggests otherwise. The right-hand and middle panels of Fig. 11 show that cluster #64 (in red) is detached from clusters #62 (in blue) and #63 (in green). In particular, the kinematics of its stars in the v_R − v_ϕ plane (middle panel) are clearly different from the rest. Furthermore, the three clusters identified in the Sequoia region are disconnected in IoM as shown in the rightmost panel of the figure. Therefore, we decided to divide substructure D into the three clusters identified: (a) Cluster #64, the most retrograde; (b) cluster #63, which we associate with what has been reported until now as Sequoia (Myeong et al. 2019); and (c) cluster #62, which displays a higher average [Fe/H] when compared to clusters #63 and #64 (see Fig. 12). Given that cluster #62 seems to be complementing cluster #63 in the v_z − v_R plane, the option that cluster #62 is really a substream of Sequoia cannot be discarded.

Fig. 11. Detailed distribution in the E − Lz (left), the vR − vz (middle), and the L⊥ − Lz (right) planes of the tentative substructure D, which is colour-coded according to the different clusters associated to it. The original sample of halo stars is shown in grey in all panels.

Fig. 12. Metallicity distribution functions of the 12 extended substructures identified in this work. These MDFs are compared to that of the whole halo population under analysis (blue shaded area) and to that of Gaia-Enceladus (grey dashed histogram). Substructure ID and number of associated stars with a [Fe/H] determination from LAMOST LRS are listed. Error bars were computed assuming Poisson statistics.

4.5. Substructure E: Helmi streams

Substructure E, which is related to the Helmi streams, is discussed in Sect. 3 and is used to show our methodology. We refer the reader to that section (as well as to Figs. 3 and 4) for more information.

4.6. Substructure F

Clusters #65 and #66 could be grouped together according to their distribution in IoM at a Mahalanobis distance of ∼3.9. Unfortunately, both clusters have very few stars, with cluster #65 comprising eight stars with [Fe/H] measurements, and cluster #66 nine such core members. Comparisons of both MDFs and CaMDs suggest that they are drawn from different distributions; see Table 10. Although this conclusion is tentative because of the small number of stars, it is supported by their different average metallicities of ⟨[Fe/H]⟩ = −1.7 and −1.3, respectively. Given these results, we do not consider substructure F to be a physically independent structure.

It is interesting to note that the location in the E − L_z plane of clusters #65 and #66 overlaps with the region that would be associated to debris from the Sagittarius dwarf galaxy (Naidu et al. 2020) and indeed with some of the globular clusters linked with Sgr in Massari et al. (2019). However, that is not the case in other IoM projections, with the clusters identified in this work displaying a higher L_⊥ than the globular clusters associated with Sgr in Massari et al. (2019). With the information at hand, it is not possible to properly further assess the nature of the clusters or their tentative association to Sgr, except to note that the average metallicity of cluster #66 is only slightly more metal-poor than that reported for the leading arm of Sgr in Hayes et al. (2020).

4.7. Cluster #3

Cluster #3 is not only detached from any other structure according to the dendrogram of Fig. 1 and our tentative Mahalanobis distance threshold, but is also composed of a large number of stars, which drives us to consider it as a potential substructure on its own. Cluster #3 is located in the low-energy part of the densely packed region where halo and thick disc cohabit. It is the largest cluster detected by the algorithm, with 3032 stars in total, and presents a statistical significance of 6.3σ (see Paper I). Of its 3032 stars, 2075 have metallicity determinations from LAMOST LRS.

Cluster #3 displays a complex MDF (see Fig. 12). While most of the stars are compatible with being thick-disc stars ([Fe/H] > −1.0, peaking at ∼ − 0.7), there exists a non-negligible population of metal-poor stars (with [Fe/H] from −1.0 to −2.2). Thus, cluster #3 contains a mix of populations, with the more metal-poor one possibly of accreted origin, although we cannot rule out the possibility that this population is associated to an old in situ halo or disc. Comparing the energy and angular momentum distributions of the globular clusters associated to Kraken/Heracles with cluster #3, we find that cluster #3 is a more weakly bound (displaying higher energies) system than Kraken/Heracles, and thus, they are not related to one another (see Kruijssen et al. 2019a; Massari et al. 2019; Horta et al. 2021).

4.8. Individual small clusters

From the analysis presented so far in the previous sections, we identify 12 large substructures, and isolate 4 small clusters that are not easily associated to any of those large structures. These are clusters #28 and #40 (originally associated with substructure A), and clusters #65 and #66 (originally postulated to form substructure F). Apart from these, there are 4 additional clusters that are not linked in IoM by our approach unless we allow a large Mahalanobis distance for the union (see Fig. 1). These clusters are #1, #2, #4, and #68, and are shown in the bottom right panels of Fig. 2. We briefly comment on each of them below:

– Cluster #1 is formed by 33 stars that tightly clump together, not only in IoM (very low energy), but especially in velocity space and on the sky. As discussed in Paper I, these stars are likely members of the globular cluster M 4.

– Cluster #2 is composed of 12 stars (8 with [Fe/H] values from LAMOST LRS) located in between the hot thick disc and the low-energy half of the halo population, with L_z ∼ 500 kpc km s⁻¹. Its stars show an overall metallicity of [Fe/H] ∼ −0.86, with a wide dispersion of 0.4 dex. Approximately 63% of its stars have a narrower metallicity range, namely [Fe/H] = [−0.75, −0.4]. The velocity distribution is particularly tight in the v_ϕ velocity component.

– Cluster #4 is also located in the region dominated by thick-disc stars (close to cluster #2), and contains only 14 stars, of which 9 have [Fe/H] determinations from LAMOST LRS. Their average [Fe/H] ∼ −1.2 and dispersion is 0.4 dex.

– Cluster #68 contains 29 stars (represented in pink in Fig. 2, panels for independent clusters) and seems to be related to the outer disc (given its z-angular momentum and low L_⊥).

We do not further consider the details of these small individual clusters (#1, #2, #4, #28, #40, #65, #66, and #68) because of their limited sample size and given that we are interested in unveiling the main constituents of the Milky Way halo. We will revise these small individual clusters together with other tight clumps in IoM and velocity space in Balbinot et al. (in prep.).

5. An updated view of the Galactic halo

The thorough analysis carried out in the previous section allows us to add support to some of the tentative conclusions of Paper I, where we identified 67 clusters with high significance in the nearby halo with a single-linkage-based clustering algorithm. By comparing the properties of these clusters, we are able to isolate 12 extended (in IoM) separate substructures, some of which are constituted by a single cluster. The substructures include the Helmi streams (substructure E, Helmi et al. 1999); three dynamically distinct clusters in the Sequoia region (clusters #62, #63, and #64, see also Myeong et al. 2019); Thamnos 1 and Thamnos 2 (substructures B1 and B2 respectively, see also Koppelman et al. 2019a); GE (substructure C, Helmi et al. 2018); one independent structure with similar angular momentum to GE but higher binding energy (cluster #6); a high-binding-energy substructure in a region of overlap between the hot thick disc and halo (cluster #3); and three subgroups (substructure A, and clusters #12 and #38) in the region occupied by the hot thick disc. In Appendix B, we present summary plots of their distributions in IoM and in velocity space. In those figures, the colour coding represents the constituent clusters from Paper I.

Apart from these 12 substructures, we discern eight individual small clusters that, based on their IoM distribution or internal properties, we are not able to assign to any of the identified structures. Their small sizes and in some cases their tightness in IoM or velocity space suggests a possible link to disrupted star clusters. These will be investigated in detail in Balbinot et al. (in prep.). Therefore, we do not consider these small clusters in the remainder of this paper.

We now proceed to the full characterisation of the 12 substructures. To this end, in Sect. 5.1 we first compare their global MDFs as well as their CaMDs using our isochrone-fitting approach, which includes determination of their average age and [Fe/H]. In Sect. 5.2 we characterise their dynamical properties via a principal component analysis (PCA), which permits the identification of additional member stars. Lastly, in Sect. 5.3, we make a first attempt at characterising their chemical properties using APOGEE data. This combined approach provides an updated view of independent substructures in the nearby Galactic halo, as discussed in Sect. 5.4.

5.1. MDF and CaMDs of the substructures

Table 11 displays the p-values of the KS tests comparing the MDFs of every pair of possible independent substructures. This table reveals clear quantitative differences in the metallicity distributions of the various substructures, particularly of substructures that are located in the region of the hot thick disc (substructure A and clusters #3, #12, and #38) compared to the rest.

We also see from Table 11 that some substructures that occupy very different regions of IoM present similar MDFs, such as Sequoia, the Helmi streams, and Thamnos 1. Some of this is expected, because the progenitors of Sequoia and the Helmi streams for example have likely been systems of comparable size (see e.g. Koppelman et al. 2019a), possibly obeying a mas-s-metallicity relation (e.g. Kirby et al. 2013). The MDFs of substructures that are closer in IoM, as in the case of Sequoia and cluster #64, could be similar because they suffer from contamination from each other (see Sect. 5.2 for further discussion). Also, interesting similarities are present between cluster #62 (which is tentatively independent of Sequoia) and Thamnos 1 and 2. Although we have not fully ruled out an association between cluster #62 and Sequoia, the MDF of the former appears to be somewhat more similar to that of smaller structures like Thamnos 1 and 2.

Figure 12 shows the MDFs of the various substructures. There clearly exists a correlation between the location of structures in the E − L_z plane and the properties of their MDF. Structures overlapping with the hot thick disc region (substructure A and clusters #3, #12, and #38) show an important contribution of higher metallicity stars (with peak [Fe/H] ∼ −0.6) associated to the halo red sequence (see e.g. Haywood et al. 2018). Cluster #3 occupies a region in IoM that is intermediate between accreted halo stars from GE and in situ hot thick disc stars. Despite showing average values of age and metallicity that could be compatible with a thick-disc population (∼10.5 Gyr and [Fe/H] ∼ −0.75), Fig. 12 clearly reveals a very important contribution from metal-poor stars to its MDF. In comparison, substructure A is much more dominated by a thick disc population, with their MDF peaking at [Fe/H] ∼ −0.6, and displaying a more tenuous tail of metal-poor stars. This is also apparent for cluster #12, but a possible interpretation in this case is that it suffers from contamination from substructure A and that the peak at lower metallicity should be considered as an independent subgroup (this is further supported by the CaMD analysis; see Table 12).

We interpret these findings as hints that an accreted, metal-poor population might be buried in the region that overlaps with the hot thick disc, although we cannot rule out the possibility of this being a relic of the old, in situ halo or disc. The chemical analysis presented in Sect. 5.3 supports this interpretation to some extent, but more information is necessary to establish this firmly.

Figure 12 shows that, with our definition of Sequoia, its MDF might hint at the existence of three metallicity peaks. These have similar properties to those reported in Naidu et al. (2020, which led these authors to divide Sequoia into three, namely L’itoi, Sequoia proper, and Arjuna), but the peaks are not statistically significant in our sample (see also the discussion in Paper I) and there is no correspondence between Arjuna, I’itoi, and Sequoia in Naidu et al. (2020) and our clusters #62, #63, and #64. The remaining substructures shown in Fig. 12 display MDFs with peaks at lower metallicities, compatible with an origin in accreted systems of lower mass.

We now focus on the CaMDs of the various substructures. The results of the isochrone-fitting procedure for each of the substructures are given in Table 12. This lists for each substructure its average age and [Fe/H], where the latter is in good agreement with the estimates based on the LAMOST LRS data. As stated above, here we also observe a clear dichotomy between accreted structures (old and metal-poor) and those with more disky dynamics (younger and more metal-rich). Figure 13 shows the CaMDs of all 12 independent substructures with their best-fit isochrone overplotted with an orange solid line. To compare with one of the main constituents making up the inner accreted halo, we also overplot the GE best-fit isochrone. We see that there are only subtle differences in the distribution of stars in the CaMD, which are nevertheless picked up with a proper quantitative analysis (via the fitting of isochrones) as can be seen in Table 12.

Fig. 13. CaMDs for the 12 extended substructures identified in this work. We include the best isochrone from our isochrone-fitting approach (see Sect. 3.1) to each structure (orange solid line) as well as that for the combined Gaia-Enceladus substructure (purple solid line).

5.2. Substructure characterisation in IoM

In this subsection we re-examine the extent of each of the individual substructures in IoM. To do this, we follow a similar approach to that presented in Paper I for individual clusters, and compute the covariance matrix of the members of each of the 12 independent substructures. Here the members we refer to consist of the union of the original members from clusters that are associated to each given substructure. The results are given in Table 13, along with the significance of each substructure, which is calculated by comparing the density of member stars to that from randomised data sets within the PCA ellipsoidal contour encompassing 95.4% of the members (see Sect. 3.3 of Paper I for a full description). To transform the values listed in this table to the corresponding physical quantities, for the means of the ellipsoids we can use that

$$\begin{array}{c}\hfill ⟨{\mu }_i⟩=\frac{2}{{\mathrm{\Delta }}_i}(⟨I_i⟩-I_{i,\min })-1,\end{array}$$(3)

where Δ_i = I_i, max − I_i, min, with i = 0…2 and I_i corresponding to E, L_⊥, and L_z, respectively, with the minimum and maximum values given by E = [ − 170 000, 0] km² s⁻², L_⊥ = [0, 4300] kpc km s⁻¹, and L_z = [ − 4500, 4600] kpc km s⁻¹, which roughly correspond to the minimum and maximum values in the halo set. For the covariance matrix, these definitions lead to

$$\begin{array}{c}\hfill {\sigma }_{\mathit{ij}}=\frac{4}{{\mathrm{\Delta }}_i{\mathrm{\Delta }}_j}{\mathrm{\Sigma }}_{I_{i,j}},\end{array}$$(4)

where Σ_Ii, j is the covariance matrix in IoM space.

Figure 14 shows the ellipsoidal boundaries of such a characterisation of the various substructures but now in the E − L_z plane. These were determined by a PCA analysis in this plane, where the contours indicate the 1.77σ extent along the principal axes of the structure (roughly encompassing 80% of a perfect Gaussian distribution).

Fig. 14. Ellipsoidal boundaries of the 12 independent extended substructures identified in this work, using a PCA. The black dashed ellipse represents the ellipsoid that would be representative of Sequoia if defined as the union of clusters #62, #63, and #64. The original sample of halo stars is shown as a grey background.

We observe in Fig. 14 that the orientation of the ellipsoid of GE is somewhat odd, as it does not conform to the expectation that the ellipsoids’ tilt approximately follows the circularity parameter (see e.g. Helmi & de Zeeuw 2000); that is, typically stars from a given accreted object ought to have higher |L_z| as they have more energetic or less-bound orbits. The tilt is likely in part due to how the random sets to assess the significance of a cluster have been constructed. Whereas in the original dataset, GE stars and the hot thick disc are clearly separated (see e.g. the distribution of grey points in Fig. 14), in our random sets they are not (see also Fig. 3 of Paper I). This is because stars from both these dominant components are being moved by the randomisation process, which consists in reshuffling their velocities to the region in between. A better way to assess the significance would therefore be to model the presence of the two components. However, this can be remedied to a certain extent by considering clusters with lower significance as well (down to 2σ in comparison to the randomised datasets; see Paper I). Furthermore, given the distribution of stars in IoM, especially the increase in density towards higher binding energies and slightly prograde orbits (i.e. positive L_z), the identification of highly significant clusters is more challenging in that area, and could therefore cause part of the tilt.

The analysis presented in Sect. 4.4 prompted us to deem the clusters located in the Sequoia region as independent, mainly based on their different dynamical properties (position in IoM planes as well as velocity distributions). Another argument against their union is that this would require a very massive progenitor (see e.g. Koppelman et al. 2019a). The black dashed ellipse in Fig. 14 represents the extent in the E − L_z plane of an alternative definition of Sequoia including clusters #62, #63, and #64. Such a large system would occupy an area similar to that of GE in the E − L_z plane (translated into a very massive system). This would be inconsistent with the average low metallicity of the three clusters ([Fe/H] ∼ −1.6, with a dispersion of 0.3 dex for clusters #63 and #64; and [Fe/H] ∼ −1.47, with a dispersion of 0.25 dex for cluster #62, and which can be roughly compared to the ⟨[Fe/H]⟩ ∼ −1.3 of GE, with a dispersion of 0.3 dex), if we assume they follow the known mass–metallicity relation (from the local to the high-redshift Universe, e.g. Erb et al. 2006; Kirby et al. 2013; Ma et al. 2016).

Interestingly, cluster #64 (the most retrograde structure in the Sequoia region) presents a distribution of stars in the CaMD that is almost ‘parallel’ to that of GE towards the more metal-poor end (as shown by its lower metallicity compared to GE; see Table 12). This could be consistent with the possibility that this substructure originates in part from the metal-poor tail of GE (indeed they show compatible ages, with cluster #64 being marginally older than GE), an interpretation that would be favoured by the simulations of Koppelman et al. (2019a). However, as recently shown in Matsuno et al. (2022), some stars that would be associated to cluster #64 have chemical properties that appear to be different from those of GE, and more similar to those of Sequoia. As shown in that paper, there is tentative evidence of contamination by GE in that region of IoM, at a level of 10%−20%. To draw firmer conclusions, a more profound characterisation of the chemistry of these clusters is much needed.

Figure 14 is similar to figures presented in Koppelman et al. (2019a) and Naidu et al. (2020) based on Gaia DR2 data, but reveals many more substructures, with an important number in the region populated by the hot thick disc. However, as pointed out in Paper I we also miss some known substructures, such as Wukong/LMS-1 (Naidu et al. 2020; Yuan et al. 2020b), likely because our sample is concentrated on the local halo, and Aleph (Naidu et al. 2020) because of our kinematic selection.

The characterisation of the 3D ellipsoids in IoM listed in Table 13 enables us to identify additional tentative members of the structures. To do this, as in Paper I and following the refinement of individual clusters in Paper I and Sect. 2, we consider stars that are located within a Mahalanobis distance D_cut ≤ 2.13 of the centre of the different independent structures as additional tentative members. There are a total of 10477 stars within the D_cut Mahalanobis distance to one of the 12 substructures, implying that approximately 20.3% of the total sample under analysis can likely be associated to these large substructures. The final number of stars associated to each independent substructure (obtained as described immediately above) is listed in Table 13 in parentheses.

5.3. Substructure chemical characterisation

We now study the behaviour of [Mg/Fe] with [Fe/H] using data from APOGEE (see Sect. 3.2) for the different substructures identified. Following the procedure described in the previous section, we consider members of the different substructures to be those stars with Mahalanobis distances D_cut ≤ 2.13 to the centre of the different substructures. To guide the eye, we fit the chemical evolution trends of substructure A and Gaia-Enceladus (representative of the hot thick disc and of the accreted population, respectively) with quadratic polynomial functions. For the fitting itself, we use the subset of core members of the clusters that were merged to form substructure A or Gaia-Enceladus, respectively (see Sects. 4.1 and 4.3). After a first fit is obtained, we remove 3σ chemical outliers, where σ includes both the uncertainty in the fit (estimated through bootstrapping) and the residual scatter. The results are shown in Fig. 15, where the shaded region denotes the 1-σ width described above. In this figure, we also plot the abundances of the individual stars associated to the various identified substructures.

Fig. 15. [Mg/Fe] and [Fe/H] abundances of the different substructures identified in this work, using data from APOGEEDR16. Filled symbols represent core members of the clusters associated with the substructure considered, while smaller open symbols are for stars that are “added” based on their Mahalanobis distance to each substructure (see Sect. 5.2 and Table 13). Chemical abundance trends of A (thick disc; purple) and GE (yellow) are shown in all panels to guide the eye. Colours are as in Fig. 14.

The bottom right panel of Fig. 15 presents chemical abundances of substructure A, and clusters #12 and #38. Substructure A has a high [Mg/Fe] ratio as expected from the suggested hot thick disc origin. While both clusters #12 and #38 overlap with the track followed by A, cluster #12 includes a few stars on the Gaia-Enceladus track. The additional chemical dimension of Mg abundance enables us to see chemical differences between the two stellar populations more clearly, highlighting the importance of studying abundance ratios.

Compared to substructure A, the substructure associated to GE has a lower [Mg/Fe] (middle right panel of Fig. 15), which is consistent with the well-known chemical abundance differences (first reported as high-Mg and low-Mg sequences respectively, in Nissen & Schuster 2010). In comparison to Fig. 9, which shows only the core members of the clusters associated to GE, we note a broader chemical sequence which is the result of adding extra tentative members likely associated with Gaia-Enceladus based on their Mahalanobis distance to the substructure. For instance, above [Fe/H] > −0.9, there are no contaminant stars (i.e. with high [Mg/Fe] > 0.2) among cluster core members (9 stars in this metallicity range), while this fraction is 19% (5 out of 27 stars) when adding members based on the Mahalanobis distance to the whole structure. This might suggest that the actual distribution of stars from Gaia-Enceladus is better captured by those stars belonging to the clusters associated with it, and that adding members to substructures assuming these follow a Gaussian distribution might lead to a lower purity. We also include substructure #6 in the same panel as GE. The two stars from this cluster seem to follow the same chemical track as GE (but small number statistics prevent us from discerning the differences seen in their MDFs).

The top left panel of Fig. 15 is for the Helmi streams. It is not possible to come to any conclusions on the chemical properties of this substructure as there is only one star with abundance information. Sequoia and cluster #64 are shown in the top right panel of Fig. 15. It is hard to decipher whether or not #64 and Sequoia have conclusively different abundance ratios because of the small number of member stars in APOGEE. It is nevertheless interesting that there is a tendency for Sequoia and #64 stars to have lower [Mg/Fe] at [Fe/H] ∼ −1.5 than Gaia-Enceladus, which is consistent with Matsuno et al. (2019, 2022) and Monty et al. (2020). This extremely low [Mg/Fe] ratio suggests that the progenitor galaxies for the highly retrograde substructures in the nearby halo were less massive and that they have potentially experienced inefficient star formation such that the chemical enrichment by type Ia supernovae becomes apparent at a lower [Fe/H] than for GE (further supporting our division of Sequoia into three identified clusters). Unfortunately, this scenario is based on information from only about a dozen stars with chemical information. More follow-up chemical abundance studies should be carried out in order to confirm or rule out this hypothesis. Finally, and as already noted on the basis of their MDFs and CaMDs, Sequoia and #64 have similar stellar populations.

The two other significantly retrograde substructures (but on lower inclination orbits), Thamnos 1 and Thamnos 2, do not seem to have very low [Mg/Fe] values (middle left panel of Fig. 15). Thamnos 2 has, on average, higher [Mg/Fe] than Gaia-Enceladus when compared at a fixed metallicity, while Thamnos 1 seems to be roughly consistent with the Gaia-Enceladus trend. If Thamnos 1 and Thamnos 2 are debris from small dwarf galaxies, their high [Mg/Fe] might not be unexpected given their high binding energies. This is because, in order to bring the debris from such low-mass galaxies to the inner Milky Way, the accretion event has to happen quite early on (Amorisco 2017). By that time, there would not have been much time for a significant number of type Ia supernovae to contribute to the chemical evolution in the system.

The low-energy cluster #3 contains stars that overlap both with the hot thick disc and the Gaia-Enceladus/accreted population in terms of chemical abundance ratios (bottom left panel of Fig. 15). We note however that a fraction of its stars fall in between the sequences.

5.4. A global view

In this section, we take a step back and present a more global view of the nearby Galactic halo provided by the various substructures we have identified.

Figure 16 plots the different substructures in the E − L_z and L_⊥ − L_z planes colour-coded according to their age and metallicity as presented in Sect. 5.1. This figure shows an interesting trend. Apart from the already known tendency of metal-poor stars to be on more retrograde orbits (e.g. Carollo et al. 2007, or more recently, Koppelman et al. 2019a), we find that the average age of the retrograde substructures tends to be higher than those on more prograde orbits (with the notable exception of the Helmi streams, showing an isochrone average age of ∼12 Gyr). This trend is not the result of an age–metallicity degeneracy (which in fact would go in the opposite sense), because the metallicities from isochrone fitting are fully consistent with those obtained from spectroscopic surveys.

Fig. 16. Distribution of the independent substructures detected in this work in the E − Lz (top) and the L⊥ − Lz (bottom) planes, colour-coded according to their isochrone fitting age (left) and [Fe/H] (right). This information is displayed in Table 12. In all panels, the grey background represents the original sample of halo stars analysed in this work.

The trend seen in Fig. 16 has not yet been reported in the literature, either in observations or in numerical simulations, as far as we are aware. As objects on retrograde orbits show weaker dynamical coupling to the host galaxy, one might expect that such objects would be accreted later (in which case they would have younger populations which is opposite to what is observed). A possible explanation is perhaps that the structures on retrograde orbits were part of a group that fell in together, for example, with GE. In that case, pre-processing in the group would prevent star formation in the smallest objects, and could lead to an outcome such as that observed.

6. Conclusions

In this paper we thoroughly explore the outcome of a data-driven clustering algorithm (see Paper I) applied to a nearby sample of halo stars (up to 2.5 kpc from the Sun) with 6D phase-space information. We performed a detailed analysis of the internal properties (i.e. metallicity distribution functions, CaMDs, and chemical abundances) of the clusters identified by the algorithm in IoM space to robustly assess their independence and to establish and confirm possible links among them. This allowed us to identify 12 large independent substructures encompassing approximately 20% of the halo stars in our sample. The main results can be summarised as follows:

• The large substructures identified include the Helmi streams, Sequoia, Gaia-Enceladus, Thamnos 1, and Thamnos 2. Given our sample definition and its spatial extent, other known substructures such as Wukong or Aleph are missing. Although we find clusters (substructure C2 and cluster #14; see Sect. 4.3) with similar energies and more retrograde orbits than GE, we cannot associate these to the substructure known as Rg5 (Yuan et al. 2020a), because of their substantially different L_⊥. Furthermore, with the information at hand, both substructures are fully compatible with being part of GE.

• We find that the most retrograde substructures in our sample (clusters #62, #63, and #64) display kinematic properties that differ from one to another. We hypothesise that stars populating cluster #64 (the most retrograde one) might stem in part from the metal-poor tail of Gaia-Enceladus and thus may not be fully related to Sequoia. Proper chemical abundance characterisation of many more stars from these substructures is needed to understand this portion of the halo.

• The largest (in number) substructure in our sample is constituted by the single cluster #3. It presents a clear mixture of stars with metallicities characteristic of the (hot) thick disc and a prominent metal-poor component. Our chemical abundance analysis based on APOGEE hints at the presence of an intermediate population between the high-Mg and low-Mg sequences. This suggests the existence of a buried accreted population with high binding energy (but lower than that of Kraken/Heracles, Kruijssen et al. 2019a; Horta et al. 2021, and hence apparently unrelated).

• The large substructure A is clearly dominated by the hot thick disc, but also depicts a small contribution of a metal-poor, halo-like population. This substructure is formed by the union of ten different clusters with compatible stellar populations. The existence of such a large amount of lumpiness within the hot thick disc (also supported by the presence of four additional clusters in a similar region), could be due to the response of a disc to past mergers, as seen in the simulations by Gómez et al. (2012), Jean-Baptiste et al. (2017).

• Most of the independent substructures identified present statistically discernable differences in their MDFs, with some similarities supporting an origin in systems of comparable size (such as Sequoia and the Helmi streams); although in some cases contamination from neighbouring structures is likely also an issue (e.g. Thamnos 1 and 2).

• We find that substructures on more retrograde orbits are, on average, older than those in less-retrograde or prograde orbits, in addition to being more metal-poor.

• We further identify eight small clusters that, based on their stellar populations or their small size, could not be associated to the 12 main independent substructures reported above or identified as an independent substructure on their own. One of these clusters has dynamical and stellar population properties that could support a link with the Sgr streams.

In this paper, we demonstrate how the combined analysis of photometry, astrometry, and chemistry in vast samples of stars in the Milky Way can help us to discern the past history of our Galaxy. The advent of new Gaia data releases, together with upcoming spectroscopic surveys, including chemical abundance follow-ups, will tremendously boost our knowledge of the history of the Milky Way.

Acknowledgments

We thank the anonymous referee for useful comments improving the original version of this manuscript. We gratefully acknowledge financial support from a Spinoza prize. H.H.K. acknowledges financial support from the Martin A. and Helen Chooljian Membership at the Institute for Advanced Study. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research makes use of python (http://www.python.org); vaex (Breddels & Veljanoski 2018), a Python library for visualization and exploration of big tabular data sets; matplotlib (Hunter 2007), a suite of open-source python modules that provide a framework for creating scientific plots; astropy (Astropy Collaboration 2013, 2018), a community-developed core Python package for Astronomy; numpy (van der Walt et al. 2011); and jupyter notebooks (Kluyver et al. 2016). This work also made use of the Third Data Release of the GALAH Survey (Buder et al. 2021). The GALAH Survey is based on data acquired through the Australian Astronomical Observatory. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This paper has made as well use of APOGEE DR16 data part of the SDSS IV scheme. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. We have made use of RAVE data for this work, see the RAVE website at https://www.rave-survey.org. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

References

Appendix A: Individual clusters properties

Table A.1 shows an overview of the characteristics of the clusters detected by the algorithm described in Paper I after refining them using a Mahalanobis distance limit of D_ij <  2.13. We encourage the reader to refer to Paper I for a thorough characterisation of the original clusters. Cluster ID is the cluster number identified by the algorithm; N_⋆ is the number of stars included in this analysis (D_ij <  2.13); ⟨L_z⟩ represents the mean value of L_z; ⟨L_⊥⟩, the mean value of L_⊥; ⟨E⟩, the mean value of E; N_[Fe/H] is the number of stars with [Fe/H] determinations from LAMOST LRS, with ⟨[Fe/H]⟩ its average value and σ([Fe/H]) its standard deviation.

Appendix B: Extended substructure summary plots

In this Appendix we compare the distribution of stars in IoM as well as in velocity planes of the 12 independent substructures identified in this work.

Fig. B.1. Distribution of stars in several projections in IoM space as well as in velocity space for the Helmi streams (substructure E) as identified in this work. In different colours we show the different clusters making up this substructure. For reference, the grey background shows our full halo sample.

Fig. B.2. Cluster #64 in IoM and velocity space. Information as in Fig. B.1.

Fig. B.3. Sequoia (cluster #63) in IoM and velocity space. Information as in Fig. B.1.

Fig. B.4. Cluster #62 in IoM and velocity space. Information as in Fig. B.1.

Fig. B.5. Thamnos 1 (cluster #37, substructure B1) in IoM and velocity space. Information as in Fig. B.1.

Fig. B.6. Thamnos 2 (substructure B2) stars in IoM and velocity space. Information as in Fig. B.1.

Fig. B.7. Gaia-Enceladus (substructure C) in IoM and velocity space. Information as in Fig. B.1.

Fig. B.8. Cluster #6 in IoM and velocity space. Information as in Fig. B.1.

Fig. B.9. Cluster #3 in IoM and velocity space. Information as in Fig. B.1.

Fig. B.10. Substructure A in IoM and velocity space. Information as in Fig. B.1.

Fig. B.11. Cluster #12 in IoM and velocity space. Information as in Fig. B.1.

Fig. B.12. Cluster #38 in IoM and velocity space. Information as in Fig. B.1.

All Tables

All Figures

Fig. 1. Relationship between the significant clusters according to the single linkage algorithm, depicted according to their Mahalanobis distance in the IoM space (Paper I). The different colours represent a tentative classification into different substructures, namely six large groups (named A to F) and a set of individual clusters (in grey; see text for details). The order of the clusters along the x-axis is such that it facilitates the interpretation and readability of the dendrogram, and does not reflect internal properties. In our analyses (e.g. tables, etc.), we broadly follow this same order.

In the text

	Fig. 2. Distribution of stars for the tentative substructures in IoM space: E vs. Lz (top two rows) and L⊥ vs. Lz (bottom two rows). Each panel plots one of the six large groups and the independent clusters from Paper I. Colours are as in Fig. 1 except for the independent clusters panel. The full original sample upon which this study is built is in grey in all panels.
In the text

Fig. 3. Isochrone fitting procedure applied to clusters #60 and #61 together, which are associated with the Helmi streams. Left: age–metallicity plane, where the background colours show the distribution of the d2 parameter used for assessing the quality of the fit obtained. The average age and metallicity, together with the corresponding errors, are represented by the black point with errorbars. To compute these average values, we use isochrones with a d2 parameter within its lowest 10th percentile (small black dots). Horizontal dotted lines show the [Fe/H] values adopted for the isochrones. The histogram of d2 values is displayed in the bottom left for all the different isochrones considered. Right: CaMD for stars in the union of clusters #60 and #61. The best-fit isochrone is depicted in red, and all compatible isochrones (black dots in the left-hand panel) are overplotted in black with transparency. Stars bluer than the grey dashed lines (horizontal branch and blue stragglers, shown with transparency) are avoided in the fit. Every star is colour-coded according to its weight in the fit (see text for details).

In the text

Fig. 4. Analysis of the stellar populations of the Helmi streams. Left: distribution in the CaMD of the stars classified as part of the Helmi streams (clusters #60 and #61) colour-coded according to the colour distance to the best-fitting isochrone (black dashed line, see Sect. 3.1). Center: histograms of the difference in colour for stars in clusters #60 and #61, and showing compatible distributions with a KS p-value of ∼0.6. Right: metallicity distribution function from LAMOST LRS for both clusters (p-value of 0.99). Error bars were computed assuming Poisson statistics.

In the text

Fig. 5. Characterisation of substructure A and the independent clusters #12, #38, #28, and #40. Left: distribution in the L⊥ − Lz plane. The original sample of halo stars is shown as a grey background. Middle: normalised MDFs for the three largest (in number of stars) structures/clusters. As an inset we add the MDF of cluster #12 alone to highlight its apparently double population. See text for details. Right: normalised MDFs for clusters #28 and #40, showing a narrow MDF. We add the MDF of substructure A (dashed, blue line) for comparison. Error bars in all plots are computed assuming Poisson statistics.

In the text

Fig. 6. Comparison of the distribution of stars in CaMD for cluster #12 (blue) and substructure A (red). The black line represent the isochrone that best reproduces the distribution of stars in the CaMD for substructure A (see Sect. 3.1). The inset panel shows the distributions of colours with respect to the best fit isochrone for both structures. We mark low-metallicity stars ([Fe/H] < −1.4) with blue crosses.

In the text

	Fig. 7. Characterisation of the substructures identified within the tentative structure B. Left: distribution in the E − Lz plane. The original sample of halo stars is shown as a grey background. Right: normalised MDFs. We use blue for B1 (Thamnos 1) and cyan for B2 (Thamnos 2; see text for details). Error bars are computed assuming Poisson statistics.
In the text

	Fig. 8. E − Lz (top) and L⊥ − Lz (bottom) distribution of the subgroups making up GE (substructure C). As in Fig. 2, the panels present each of the groups identified and the three independent clusters. The full original sample upon which this study is built is represented as a grey background in all panels. The full extension of substructure C is highlighted in silver.
In the text

	Fig. 9. [Mg/Fe] abundance ratios as a function of [Fe/H] for substructure C (GE) from APOGEE DR16. We use different symbols for the different clusters. Colours for the groups (C1 to C5) follow Fig. 8. The small dots show abundances of all the halo stars within 2.5 kpc in APOGEE DR16.
In the text

Fig. 10. Combined metallicity distribution function of all clusters in substructure C/Gaia-Enceladus. We overlay the [Fe/H] distribution for all 34 532 stars with [Fe/H] from LAMOST LRS in our initial halo sample (containing 51 671 stars). We note that our selection of GE stars based on our clustering algorithm naturally does not suffer from contamination from the hot thick disc (red sequence). Error bars were computed assuming Poisson statistics.

In the text

	Fig. 11. Detailed distribution in the E − Lz (left), the vR − vz (middle), and the L⊥ − Lz (right) planes of the tentative substructure D, which is colour-coded according to the different clusters associated to it. The original sample of halo stars is shown in grey in all panels.
In the text

Fig. 12. Metallicity distribution functions of the 12 extended substructures identified in this work. These MDFs are compared to that of the whole halo population under analysis (blue shaded area) and to that of Gaia-Enceladus (grey dashed histogram). Substructure ID and number of associated stars with a [Fe/H] determination from LAMOST LRS are listed. Error bars were computed assuming Poisson statistics.

In the text

	Fig. 13. CaMDs for the 12 extended substructures identified in this work. We include the best isochrone from our isochrone-fitting approach (see Sect. 3.1) to each structure (orange solid line) as well as that for the combined Gaia-Enceladus substructure (purple solid line).
In the text

	Fig. 14. Ellipsoidal boundaries of the 12 independent extended substructures identified in this work, using a PCA. The black dashed ellipse represents the ellipsoid that would be representative of Sequoia if defined as the union of clusters #62, #63, and #64. The original sample of halo stars is shown as a grey background.
In the text

Fig. 15. [Mg/Fe] and [Fe/H] abundances of the different substructures identified in this work, using data from APOGEEDR16. Filled symbols represent core members of the clusters associated with the substructure considered, while smaller open symbols are for stars that are “added” based on their Mahalanobis distance to each substructure (see Sect. 5.2 and Table 13). Chemical abundance trends of A (thick disc; purple) and GE (yellow) are shown in all panels to guide the eye. Colours are as in Fig. 14.

In the text

	Fig. 16. Distribution of the independent substructures detected in this work in the E − Lz (top) and the L⊥ − Lz (bottom) planes, colour-coded according to their isochrone fitting age (left) and [Fe/H] (right). This information is displayed in Table 12. In all panels, the grey background represents the original sample of halo stars analysed in this work.
In the text

	Fig. B.1. Distribution of stars in several projections in IoM space as well as in velocity space for the Helmi streams (substructure E) as identified in this work. In different colours we show the different clusters making up this substructure. For reference, the grey background shows our full halo sample.
In the text

	Fig. B.2. Cluster #64 in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.3. Sequoia (cluster #63) in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.4. Cluster #62 in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.5. Thamnos 1 (cluster #37, substructure B1) in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.6. Thamnos 2 (substructure B2) stars in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.7. Gaia-Enceladus (substructure C) in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.8. Cluster #6 in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.9. Cluster #3 in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.10. Substructure A in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.11. Cluster #12 in IoM and velocity space. Information as in Fig. B.1.
In the text

	Fig. B.12. Cluster #38 in IoM and velocity space. Information as in Fig. B.1.
In the text