Dynamical substructures of local metal-poor halo

Based on 4,\,098 very metal-poor (VMP) stars with 6D phase-space and chemical information from \textit{Gaia} DR3 and LAMOST DR9 as tracers, we apply an unsupervised machine learning algorithm, Shared Nearest Neighbor (SNN), to identify stellar groups in the action-energy (\textbf{\textit{J}}-$E$) space. We detect seven previously known mergers in local samples, including Helmi Stream, Gaia-Sausage-Enceladus (GSE), Metal-weak Thick Disk (MWTD), Pontus, Wukong, Thamnos, and I'itoi+Sequoia+Arjuna. According to energy, we further divide GSE and Wukong into smaller parts to explore the orbital characteristics of individual fragments. Similarly, the division of Thamnos is based on action. It can be found that the apocentric distances of GSE parts of high and medium energy levels are located at $29.5\pm3.6\,{\rm kpc}$ and $13.0\pm2.7\,{\rm kpc}$, respectively, which suggests that GSE could account for breaks in the density profile of the Galactic halo at both $\approx30$\,kpc and $15\text{-}18$\,kpc. The VMP stars of MWTD move along prograde orbits with larger eccentricities than those of its more metal-rich stars, which indicates that the VMP part of MWTD may be formed by accreting with dwarf galaxies. Finally, we summarize all substructures discovered in our local VMP samples. Our results provide a reference for the formation and evolution of the inner halo of the Milky Way (MW).


INTRODUCTION
The proto-MW has undergone frequent mergers with small progenitor galaxies.According to the ΛCDM cosmological model, the MW forms through hierarchical processes (White & Frenk 1991), which predicted that the MW involves a series of accretion events.Many studies have shown that the Galactic stellar halo is primarily formed by the merging of numerous progenitor galaxies (Helmi et al. 1999;Majewski et al. 2003;Newberg et al. 2009;Helmi et al. 2018;Belokurov et al. 2018;Myeong et al. 2018b;Koppelman et al. 2019b;Myeong et al. 2019;Yuan et al. 2020b;Naidu et al. 2020Naidu et al. , 2021;;Malhan et al. 2022).The tidal disruption of the merging galaxy occurs slowly enough to form a vast stellar stream in the Galactic halo (e.g., this is the case for the Sagittarius merger; Ibata et al. 2020;Vasiliev & Belokurov 2020).However, for radially accreting galaxies, the stars will quickly get phase-mixed and no clear signature of the stream will be visible (e.g., Gaia-Sausage/Enceladus Belokurov et al. 2018;Helmi et al. 2018).Fortunately, they retain some of the orbital and chemical characteristics of their progenitor galaxies.Therefore, in order to explore the tidal debris from their progenitor galaxies, a common approach is clustering in the integrals of motion space (e.g., actions).
Massive amounts of all-sky high-quality photometry and astrometry have been provided by the Gaia Data Releases (Gaia Collaboration et al. 2016Collaboration et al. , 2023)).For more complete information, its precise proper motion and parallax can be combined with chemical information and radial velocities from other spectroscopic surveys, such ★ E-mail: ducuihua@ucas.ac.cn as LAMOST (Cui et al. 2012), APOGEE (Majewski et al. 2017;Abdurro'uf et al. 2022), RAVE (Steinmetz et al. 2020a,b), SEGUE (Yanny et al. 2009;Alam et al. 2015), GALAH (De Silva et al. 2015;Martell et al. 2017;Buder et al. 2021).Thanks to these large stellar spectroscopic surveys, the combination of motion and chemical abundances has become available for millions of stars in the solar neighbourhood, which allows us to detect substructures in the local halo.
Many studies have focused on identifying and describing the MW's stellar populations, and try to determine whether they originated from the MW (i.e., in-situ halo) or nearby dwarf galaxies.Detection of substructure in phase space has worked successfully for the identification of accretion events.Helmi et al. (1999) found that the Helmi Streams originate from a dwarf galaxy of ∼ 10 8 M ⊙ that was accreted 5-8 Gyr ago (Helmi et al. 2006;Koppelman et al. 2019a).Other prominent discoveries are the core of Sagittarius (Sgr) dwarf galaxy (Ibata et al. 1994(Ibata et al. , 1995;;Yanny et al. 2000;Vasiliev & Belokurov 2020) and it is still forming tidal stream (Ibata et al. 2001;Majewski et al. 2003;Belokurov et al. 2014;Hernitschek et al. 2017;Ramos et al. 2020), which have been used as a paradigm for how dwarf galaxies merger with the MW.The majority of the inner halo comes from the merger of a massive dwarf galaxy with an initial stellar mass of 5 × 10 8 −5 × 10 9 M ⊙ known as Gaia-Sausage-Enceladus (Belokurov et al. 2018;Helmi et al. 2018;Haywood et al. 2018;Mackereth et al. 2019;Vincenzo et al. 2019), it was accreted by the MW at  ≈ 2 (∼ 10 Gyr ago, Di Matteo et al. 2019;Gallart et al. 2019).Furthermore, in recent years, several phase-space substructures believed to be remnants of merged galaxies have been identified in the Galactic stellar halo, including Sequoia (Myeong et al. 2018a(Myeong et al. , 2019)), Tham-nos (Koppelman et al. 2019b), Wukong/LMS-1 (Yuan et al. 2020b;Naidu et al. 2020), I'itoi, Arjuna, Aleph (Naidu et al. 2020), Heracles (Horta et al. 2021), Icarus (Re Fiorentin et al. 2021), Cetus (Newberg et al. 2009;Thomas & Battaglia 2022), and Pontus (Malhan et al. 2022).In addition, Wang et al. (2022) used the friends-of-friends algorithm to identify substructures in five-dimensional space, i.e., eccentricity , semimajor axis , the direction of the orbital pole ( orbit ,  orbit ) and the angle between apocenter and the projection of -axis on the orbital plane    , and found three remaining unknown substructures and one of them has large angular momentum and a mean metallicity -2.13 dex.
The MW's potential presumably evolved adiabatically.Considering that the actions J are conserved for a long time if the potential evolved adiabatically, those stars merged from a progenitor galaxy would be clumped in actions J space.Many recent studies (e.g., Li et al. 2019;Koppelman et al. 2019b;Yuan et al. 2020a;Li et al. 2020;Naidu et al. 2020;Malhan et al. 2022) refer to energy  as an extra quantity to distinguish the mergers even if it is not adiabatic, and their results are highly reminiscent of the substructures from numerical simulations of the disruption of satellite galaxies by the Galactic potential (Helmi & de Zeeuw 2000;Gómez et al. 2013;Orkney et al. 2023).
In this paper, we use the two-stage SNN to identify dynamical relics in the VMP local halo.In Section 2, we describe the selection criteria and dataset.We introduce the clustering method in Section 3, and then analyze all the substructures in Section 4. Finally, conclusions are given in Section 5.

DATA
In this paper, we use proper motions from Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) as well as chemical information (i.e., [Fe/H], [/Fe]) and line-of-sight velocities from LAMOST DR9 (Cui et al. 2012;Zhao et al. 2012;Liu et al. 2020).We apply selection cuts for Gaia DR3: where the last three conditions are for reliable parallaxes zero-point, and then cross-match with the LAMOST DR9 catalogue (< 1 ′′ ) to obtain metallicities and line-of-sight velocities.By comparing line-of-sight velocities between Gaia DR3 and LAMOST DR9, we notice a systematic bias of 5.32 kms −1 with a dispersion  rv = 7.26 km s −1 , which we add to the LAMOST velocities, and remove stars with | los,Gaia −  los,LAMOST − 5.32| > 3 rv that may be mismatched or have unreliable line-of-sight velocities.Finally, we utilize the following selection cuts: where  los,Gaia and  los,LAMOST are the line-of-sight velocities of Gaia and LAMOST.Finally, there are 4098 VMP stars with full 6-D phase-space information.We apply the method in Bailer-Jones et al. (2021) to estimate geometric distances, in which we sample the posterior probability with Goodman & Weare's affine-invariant Markov Chain Monte Carlo (MCMC, Goodman & Weare 2010), using the emcee python package (Foreman-Mackey et al. 2013).When calculating the distance, we take into account the correlation coefficient between parallax and proper motion (e.g., Du et al. 2019;Yan et al. 2020;Li et al. 2020;Zhu et al. 2021;Liao et al. 2023).
We adopt a right-handed Galactocentric frame of reference similar to the one defined in Li et al. (2019): here ,  and  indicate the Cartesian coordinates;  is the cylindrical radius,  is the spherical radius, and  and  represent the azimuthal and zenithal angle, respectively.In this coordinate system, the Sun is located at (8.2, 0, 0.025) kpc (Jurić et al. 2008;Bland-Hawthorn & Gerhard 2016).We take  LSR = 232.8km s −1 (McMillan 2017) for the local standard of rest (LSR) and ( ⊙ ,  ⊙ ,  ⊙ ) = (10, 11, 7) km s −1 (Tian et al. 2015;Bland-Hawthorn & Gerhard 2016;McMillan 2017) for the Sun's proper motion with respect to the LSR.
In this study, we use the McMillan (2017) model because the predicted velocity curve of this model is more consistent with the measurements of the Milky Way (e.g., Bovy 2020;Nitschai et al. 2021).To compute actions (J), energy (), and other orbital parameters, we make use of the galpy module (Bovy 2015).We analyze actions (  ,   ,   ) in cylindrical coordinates system, where   is equal to the  component angular momentum (  ) under McMillan (2017) model, and a negative   represents prograde motion.  and   describe the extent of oscillations in cylindrical radius and  directions, respectively.In addition, we obtain other orbital parameters, including eccentricity, apocentric distance, pericentric distance and the maximum vertical height (,  apo ,  peri ,  max ).For each star, we construct a set of 100 initial conditions using a Monte Carlo technique considering the observational uncertainties in heliocentric distances, proper motions, and line-of-sight velocities.The final dynamical parameters are taken as the means of the derived distributions, and associated uncertainties are the corresponding standard deviations.Note that we calculate the 4 × 4 covariance matrices about energy and actions by biweight midcovariance in the Python module astropy (Astropy Collaboration et al. 2013) because it can remove spurious outliers.

METHOD
In this work, we search for structure among the VMP stars from LAMOST DR9 in the energy-action space.Energy is conserved as long as the potential of the Milky Way is static, and the actions are insensitive to the slow, adiabatic time-dependence of the potential.Halo stars belonging to the same structure, even when they are scattered across the sky, retain similar coordinates in the energy-action space, so they could be revealed through clustering algorithms.
In this section, we summarize the method of two-stage SNN clustering; further details can be found in Chen et al. (2020).The SNN (Ertöz et al. 2003) can be seen as a modified version of DBSCAN, short for density-based spatial clustering of applications with noise (Ester et al. 1996).The SNN adopts the Jaccard distance metric of the nearest neighbours in J space to make the density threshold more flexible.The Jaccard distance is defined as: where   and   represent the sets of the nearest neighbors in J space for two stars,  and .|  ∩   | and |  ∪   | represent the sizes of the intersection and union of the nearest neighbours of two stars, respectively.Thus, for those sets that are completely overlapping,  Jaccard equals 0, whereas for those sets that are entirely disjoint,  Jaccard equals 1.
The first stage consists of three steps: retrieve the same number (N J ) of nearest neighbours for each star in actions space, and then remove those with dissimilar energies (>   ) in order to keep neighbours that share both similar actions and energies; compute the Jaccard distance matrix; cluster by DBSCAN implemented in Python module scikit-learn1 (Abraham et al. 2014).Two free parameters are involved in first-stage clustering: one is the number (N J ) of the nearest neighbours for each star in the action space, and the other is the maximum energy difference (d  ) between a star and its nearest neighbours in actions space.For a large N J , the elements of the Jaccard matrix tend to be larger, which, in turn, increases the sizes of groups composed of stars with similar actions.However, such a small set of the nearest neighbour leads to a weak correlation among most of the stars, consequently resulting in very few or even no clusters.In addition,   is a secondary parameter, as it imposes energy constraints on the nearest neighbours composed of N J members.Excessive   values loosen these constraints, while extremely small values significantly decrease the size of the nearest neighbour.
Therefore, adjusting these parameters can indeed be a formidable challenge.
We apply a Monte Carlo approach to accomplish parameter tuning and calculate the frequentist stability of groups (the number of times a group appears) and the probability of members belonging to an assigned group (the percentage of times a star appears in the duplicates of an assigned group).First of all, the maximum possible value of   should not be set too high initially, because some adjacent substructures in the   -  space exhibit small energy differences, such as Wukong and GSE.Additionally, in cases where N J is too large, it usually leads to stars that are not strongly correlated in action space being assigned to a common group.Similar to Chen et al. (2020), we set the ranges of N J and d  to 20 − 220 (each star has a chance to connect to at most about 5% of the sample closest to it in the action space) and 0 − 6600 km 2 s −2 (twice the average error of energy), and then redraw 10,000 parameter values from uniform distributions, respectively.Taking into measurement errors, we resample the energies and actions of our sample stars according to their 4 × 4 covariance matrices about energy and actions.We ran the SNN algorithm 10,000 times with randomly generated samples and parameter values, and we obtained 41,229 stellar groups.Of course, not all of these ∼ 40, 000 groups are reliable, so we need to cluster these groups, regarding each group as an individual, in order to identify high-frequency or stable groups.
In the second stage of clustering, we calculate the Jaccard matrix using these stellar groups instead of the nearest neighbour of each star in J space, and then apply DBSCAN to cluster again.In order to remove stars assigned to, with low probabilities, groups as well as to exclude groups that appear by chance due to specific N J and   values, we must set probability and stability thresholds.To avoid confusion from fortuitous overlaps, we select stars with stability (the number of times a group appears) above 30 and a probability (the percentage of times a star appears in the duplicates of an assigned group) > 30% so that only a few stars are assigned to multiple groups.Even for the group with a stability of 31, the probability of its members (∼ 73%) is much larger than the probability threshold (30%), indicating that the number of its members is not large enough to be significantly recognized multiple times but it is highly compact.For a group composed of members with a probability of ∼ 30%, its stability is outstanding.Finally, we summarize 26 groups with 1, 515 stars, as listed in Table 1, simultaneously removing 2, 583 stars that probabilistically ( probability < 30% ) belong to poor-stability ( stability < 30 ) multiple groups and represent smooth, diffuse VMP background halo without clumping in our 4D phase space.Smooth stellar halo is defined as background halo obtained by removing obvious substructures.It is utilized for measuring the anisotropy profile (Bird et al. 2021), fraction of GSE (Wu et al. 2022b), density shape (Wu et al. 2022a), and Milky Way mass (Bird et al. 2022).Smooth VMP background halo stars also involve member stars from ancient mergers, but further determination of their origins requires detailed chemical information or model fitting.
Considering that for some mergers, their member stars could be multiply stacked in the 4D phase space, such as energy wrinkles and phase-space folds from GSE (Belokurov et al. 2023;Wu et al. 2023), we cannot simply think that each group corresponds to a distinct and individual merger.Based on orbital properties, we divided the 26 groups into seven substructures, including Gaia-Sausage-Enceladus, Thamnos, Metal-weak Thick Disc, Helmi Streams, Wukong/LMS-1, I'itoi+Sequoia+Arjuna, Pontus, as follows: For the subdivided substructures (e.g., GSE, Thamnos, Wukong), each part is represented in a different colour but with the same shape.
Note that we select substructures by constraining the mean orbital parameters in each group, so a few stars belonging to a substructure cannot satisfy its selection criteria.In this work, we did not identify previously discovered five substructures, including Sagittarius (Sgr), Heracles, Aleph, Icarus, and Cetus due to: (i) the heliocentric distances of Sgr stars generally are beyond 10 kpc (Hayes et al. 2020)  .Seven substructures VMP stars described in the (J, E) space, using the same colour and shape scheme as Figure 2. Note that the retrograde GSE part (royal blue circle) coincides with the polar Thamnos (blue square) in the -  space, but does not happen in the projected action space (bottom-left).

RESULTS
In this section, for each substructure discussed below, we analyze its dynamical properties (J, ) and orbital parameters.In particular, we show all the mergers in "projected action space" represented by a di- The reason for using the projected action space is that this plot is effective in separating objects that lie along circular, radial, and inplane orbits, and it is considered to be superior to other kinematic spaces (e.g., Lane et al. 2022).

Gaia Sausage/Encelaus
It has been suggested that the structure known as the Gaia-Sausage-Enceladus (GSE) is the remnant of the Galaxy's last major accretion event (Belokurov et al. 2018;Helmi et al. 2018;Myeong et al. 2018b;Haywood et al. 2018).Wu et al. (2022b) recently found that the GSE contributed about 41% − 71% of the inner ( < 30 kpc) stellar halo by fitting their K giant sample with the Gaussian mixture model.The GSE members can be distinguished in velocity space, as they form an elongated distribution in   around a close-to-zero azimuthal velocity (Koppelman et al. 2018).In addition, due to their high eccentricities, the GSE stars can be selected using their orbital properties (e.g., Myeong et al. 2018b;Naidu et al. 2020;Bonaca et al. 2020).Naidu et al. (2020) defined the GSE as the highly radial population and selected the GSE stars by requiring  > 0.7.We preliminarily selected 17 groups by  > 0.7, and then removed a slightly retrograde group belonging to Pontus (Malhan et  space, which they interpret as folds of GSE tidal debris as it stretches and winds up due to phase-mixing in the MW gravitational potential. In addition, they found that the (,   ) distribution appears more prograde at high energies but is roughly symmetric with respect to   at low energies.Further numerical simulation shows that the least bound population of GSE-like satellite have on average slightly positive   , inheriting it from the satellite's orbital angular momentum, and more tightly bound debris are mostly located closer to   = 0.
In the top panel of Figure 1, based on the GSE VMP stars, the chevrons can also be seen, indicating that they are strongly related to energy.The angle between two wings of chevron decreases with energy, which is consistent with that in Belokurov et al. (2023).In the bottom panel of Figure 1, we additionally find that the (  , ) distribution exhibits the short and wide chevrons, and their angles increase with energy.The GSE stars are symmetrically distributed relative to   ∼ 0 between  ∼ −1.To sum up, we believe that energy is a highly significant quantity for understanding the internal structure of GSE, so we divide GSE into three parts according to energy.Firstly, a prominent high-energy GSE group (hot pink) is separately regarded as one part because it is also isolated in the ( apo , ) space, as shown in the top-left panel of Figure 2. Around  rapo ∼ 10 kpc, the GSE groups assemble into a large clump in the ( apo , ) space, corresponding to a major clump with  < −1.5 [×10 5 km 2 s −2 ] in the (,   ) space that would transform into a smooth asymmetric peak found by Belokurov et al. (2023) if the dispersed and smooth GSE halo stars in this space are considered.In order to ensure sufficient stars in the medium and low energy parts and to understand the distinct properties on both sides of the peak, we divide the clump or peak around  apo ∼ 10 kpc or  ∼ −1.7 [×10 5 km 2 s −2 ] into two parts (pink and cyan circles represent the medium and low energies, respectively).A slightly retrograde group (royal blue) has a mean eccentricity greater than 0.7 and lies on the boundary between GSE and Thamnos in the 4D space, so we consider it as an individual GSE part for comparison with the major GSE and Thamnos.
In Figure 2, the eccentric GSE groups are located at the bottom of (  −   )/ tot and   / tot space, indicating that GSE stars move on extremely radial orbits.Their pericentric distances aggregate around  peri ≈ 1 kpc, but the GSE stars definitely do not build up here because their radial velocities reach their maximal values when the stars on eccentric orbits approach their pericentres.Naidu et al. (2021) and Belokurov et al. (2023) have both successfully reproduced the generic features of GSE in phase space and spatial position through simulations of radial accretions of massive galaxies, despite discrepancies concerning stripped GSE's outer disc.While the overall characteristics of GSE are widely recognized, it is also worthwhile to study its internal structure, as we analyze below.In Figure 2, we found that the high-energy GSE (hot pink) stars constitute an obvious apocenter pile-up in the range of  apo = 29.5 ± 3.6 kpc and move on orbits with strongly radial action ((  −   )/ tot = −0.89± 0.11,   = 2199.86± 324.72 kpc km s −1 ) and high eccentricity ( = 0.93 ± 0.04).Therefore, the high-energy GSE stars could be responsible for the break in a broken power-law density profile at  ≈ 28−30 kpc (Naidu et al. 2021;Han et al. 2022;Ye et al. 2023) that also is predicted by chevrons in the -  space (Belokurov et al. 2023;Wu et al. 2023).Instead, the apocentric distance of the low-energy GSE (cyan) is the smallest, and its orbit has relatively low eccentricity ( = 0.84 ± 0.07) and radial action ((  −   )/ tot = −0.71± 0.17,   = 569.81± 83.42 kpc km s −1 ) compared to the high/mediumenergy GSE.The medium-energy GSE (pink) part with the largest proportion has an intermediate apocentric distance, eccentricity and radial action ( apo = 13.0±2.7 kpc,  = 0.88±0.06,(  −   )/ tot = −0.84±0.12,  = 921.27±272.08 kpc km s −1 ), so they could be associated with the apocenter pile-up that creates the break at 15-18 kpc (Naidu et al. 2021).For the three GSE parts with different energy levels, we suggest that the part with higher energy has a farther apocentric distance and larger eccentricity but a smaller (  −   )/ tot and radial action   , which indicates that its members move on more radial orbits, because chevrons in the -  space show high radial velocities for the high-energy GSE stars.
In Figure 3, we give all the mergers in dynamical space for easy comparison of their dynamical properties.The maximum vertical height of high-energy GSE part (hot pink circle) has very large uncertainty, because the high-energy GSE group (hot pink circle) is relatively scattered in the action-energy space, as shown in Figure 3, especially   and   , but it is significantly higher than those of medium and low due to much farther apocentre.The slightly retrograde GSE part has a somewhat larger pericentric distance than those of the other parts but smaller than that of Thamnos, so it could belong to the GSE merger.Figure 4 shows the 3D kinematic properties of GSE (hot pink, pink, cyan, and royal blue circles).The major difference of distributions of velocity for three distinct energy GSE (hot pink, pink, cyan) is about radial velocity   .Note that the cylindrical coordinate system is used here to compare with action J.The slightly retrograde GSE part (royal blue) is positioned between the other three GSE parts and Thamnos in the action-energy space or (  ,   ) space, so it could be a mixture of the GSE and Thamnos.However, the GSE members could be dominated due to their high eccentricities ( = 0.73 ± 0.06).
In order to understand the correlation between the eccentricity and velocity for the four GSE parts, we also show the GSE member stars in the velocity-eccentricity space in Figures 5 and 6.Note that the spherical coordinate system is used here to understand the dependence of radial (  ) and tangential ( √︃  2  +  2  ) velocities on eccentricity.We found that the eccentricities of GSE stars with high radial velocities slowly decreases with |  |, because, for a star moving with high radial velocity   at low zenithal angles, its orbit tends to be radial, while fast rotation (  ,   ) makes its shape more circular under the assumption that they orbit the Galactic centre.Since there is no significant difference in the zenithal velocity for the four parts, we can only determine the influence of   and   on eccentricity for a given level of   in the solar neighborhood.We conclude that, for the local GSE, the energy level is positively correlated with eccentricity or radial degree of orbit, and vice versa for   .Koppelman et al. (2019b) suggested the existence of a significant, retrograde substructure, which could be divided into two parts on the basis of metallicity and azimuthal velocity, named Thamnos 1 and Thamnos 2. We select four groups associated with Thamnos from the residual 10 groups.In order to study the fine dynamical structures of Thamnos, we split it into two parts: (  −   ) < 0 (orange square) and (  −   ) > 0 (blue square), as shown in Figures 2-6.When comparing the two parts of Thamnos, we found no peculiar difference in energy, pericentric distance, apocentric distance and eccentricity, while significant differences in azimuthal velocity (  ) and vertical velocity (  ), leading to their separation in the   -  space as shown in Figure 3.The distribution of the members of the Thamnos part (blue square) with (  −   )/ tot > 0 in (  −   )/ tot vs.   / tot space is consistent with that of Naidu et al. (2020).Koppelman et al. (2019b) found that Thamnos 2 embraces Thamnos 1 in the   -  space but Thamnos 1 has higher energy and azimuthal velocity than Thamnos 2. In the top-right panel of Figure 4, we found that the Thamnos part with (  −   )/ tot > 0 cannot fully embrace another Thamnos part with (  −   )/ tot < 0 in the   −   space, and that their azimuthal velocities (  = 84.84± 24.88 km s −1 , 120.07 ± 34.14 km s −1 ) are lower than those (  ∼ 150 km s −1 , 200 km s −1 ) of Koppelman et al. (2019b), which might be caused by their clustering in a different 4D space (i.e., ,   ,  and [Fe/H]) to probe the high-density regions of Thamnos and other differences.Koppelman et al. (2019b) divide Thamnos into two parts with the different metallicities and azimuthal velocities in the -  space, and do not clearly classify stars at low energies ( < −1.5 [×10 5 km 2 s −2 ]) and |  |, and we split it into more polar and radial parts in the projected action space, with distinct vertical and azimuthal velocities.In addition, they centred the clustering results from HDBSCAN (McInnes et al. 2017) and set boundary lines in the -  space to select Thamnos stars, while we directly used the clustering results from SNN as Thamnos stars because we are not sure whether the smooth stellar halo stars in these areas belong to Thamnos, maybe they are the slightly retrograde GSE stars.Finally, the Thamnos stars in the range of   ∼ [1200, 1600] kpc km −1 is actually Thamnos 1 proposed by Koppelman et al. (2019b).

Thamnos
Below we analyze the differences between polar and radial Thamnos parts.We know that the maximum vertical height ( max ) is related to vertical velocity (  ), that is, a star with large   usually can reach a high  max , so the polar Thamnos (blue square) stars move along orbits reaching higher  max .In the -  space, for polar Thamnos, it is evident that  decreases as |  |, whereas for radial Thamnos, there appears to be no significant correlation with zenithal velocity, somewhat resembling the behaviour observed in the GSE.However, in the -  , they appear similar descending characteristics, with the -  slope of radial Thamnos (orange square) being comparable to that of the retrograde GSE (royal blue circle), while that of polar Thamnos (blue square) is significantly steeper.Within -  space, radial and polar Thamnos parts have similar distributions, their radial velocities differ significantly from GSE because they are better suited to a single Gaussian rather than a mixture of two Gaussians describing radial anisotropic GSE (Lancaster et al. 2019;Necib et al. 2019;Iorio & Belokurov 2021;Wu et al. 2022b) for radial velocity (  ).Compared with the GSE, Thamnos is a retrograde component with low radial velocity and action, so its member stars move along orbits with smaller eccentricity  = 0.51 ± 0.15.Furthermore, we notice that the slightly retrograde GSE part is very close to Thamnos in phase spaces, recent studies (e.g., Koppelman et al. 2020;Naidu et al. 2021) show that the GSE stars lost early have large retrograde motions and a subset of these stars have low eccentricities, so these retrograde populations, such as Thamnos, Arjuna and Sequoia, may come from the same progenitor galaxy as GSE itself, which will require observational details and high-resolution simulations to be further confirmed in the future.

Metal-weak Thick Disk and Helmi Streams
Here, we select the Metal-weak Thick Disk (MWTD, Carollo et al. 2019) group based on its dynamical properties described by Naidu et al. (2020) in (  −   )/ tot vs.   / tot space.However, we found that the VMP stars of the MWTD have larger eccentricities ( = 0.61±0.07> 0.47) than those of more metal-rich ([Fe/H] ∼ −1.12) MWTD members in Naidu et al. (2020) even if we add the criterion of [/Fe] in Naidu et al. (2020), namely 0.25 < [/Fe] < 0.45, which suggests that the orbital properties of the MWTD may be strongly correlated with chemical abundance, and the VMP part of the MWTD may originate from accreted galaxies, which require high-completeness samples for further verification.After removing the MWTD (  = −142.30± 28.51 km s −1 ) from the two extremely prograde groups, the only left group (  = −143.95± 27.86 km s −1 ) could belong to the Helmi Streams, because its actions and energy, as listed in Table 2, are similar to the Helmi Streams obtained by Yuan et al. (2020a) using VMP stars.
The discovery of the Helmi Streams as accreted substructures in the halo was one of the first instances achieved via integrals of motion due to their high vertical velocities (Helmi et al. 1999).In order to further make sure that the group belongs to the Helmi Streams, we use the orthogonal and  components of the angular momentum (i.e.,  ⊥ = √︃  2  +  2  ,   ) to identify the Helmi Streams, similar to Koppelman et al. (2019b) and Naidu et al. (2020).The Helmi Streams selected here are satisfied with the selection criteria of Naidu et al. (2020), namely −1.7 <   /[10 3 kpc km s −1 ] < −0.75 and 1.6 <  ⊥ /[10 3 kpc km s −1 ] < 3.2, as shown in the bottom-right panel of Figure 4, and are of higher energy than the MWTD.Although both the Helmi Streams and MWTD are highly prograde and are adjacent in the (,   ) and (  ,   ) spaces, the vertical velocities and actions of the Helmi Streams are much larger than those of MWTD.Therefore, the MWTD stars move along radial and relatively eccentric orbits near the Galactic disc, and the Helmi Streams appear to extend out of the Galactic disc, which indicates that the progenitor galaxy of Helmi streams, with a low initial tangential velocity, might accrete with the Milky Way in a direction nearly perpendicular to the Galactic disc.In addition, the MWTD has a richer [/Fe] than that of Helmi Streams ([/Fe] MWTD = 0.29 ± 0.17, [/Fe] Helmi Streams = 0.21 ± 0.16), indicating they do not appear to originate from a common progenitor galaxy.

Wukong
It is found by Naidu et al. (2020) and Yuan et al. (2020b) that Wukong/LMS-1 is a less prograde and more metal-poor merger compared to the Helmi Streams, and its member stars also move along polar orbits due to their large vertical actions.Malhan et al. (2022) apply globular clusters, dwarf galaxies, and streams as tracers to probe accretion events and also determine that Wukong has a slight prograde motion and its objects have very polar orbits, and that there are three most metal-poor streams of our galaxy belonging to Wukong.Therefore, we regard the only two left groups with slightly prograde rotation as Wukong, and the distributions of their member stars in (  −   )/ tot vs.   / tot space and (,   ) spaces are very consistent with those in Naidu et al. (2020).Similar to the GSE, we divide Wukong into two parts with an energy difference of ∼ 0.2 × 10 5 km 2 s −2 to understand its internal structure under distinct energies.Two Wukong parts are not satisfied with the selection criteria of the Helmi Streams (Naidu et al. 2020), due to their low  component of angular momentums, and Wukong can reach a vertical height much higher than that for MWTD.Therefore, we can exclude the possibility that they belong to the Helmi Streams and MWTD.
Below we analyze the respective features and relationship between high (purple up triangle) and low (dark cyan up triangle) energy Wukong parts.Their   -  distributions almost overlap, but a significant discrepancy is in vertical action due to distinct absolute vertical or polar velocity (|  | or |  |).Wukong overlaps with the prograde GSE in -  space, but their prominent  max and   as well as unremarkable  all imply that, observationally, it is very different from GSE.However, the massive GSE progenitor merged with the Milky Way motivates the significant motion of the host galaxy, which could disrupt the precession of the angular momentum of the host disc, and after the merger, it may continue to precess.Therefore, we cannot entirely rule out the probability that the Wukong progenitor comes from a metal-poor area of GSE.
Finally, we combine the aforementioned GSE and Thamnos to discuss the correlation between pericentric distance and rotation velocity.We have found that the slightly retrograde GSE part (royal blue) possesses a farther pericentric distance compared to other GSE parts of distinct energy levels.Here, the Wukong part (purple up triangle) with a larger |  | also shows a larger pericentric distance (i.e.,  peri,purple = 6.05 ± 1.44 kpc,  peri,darkcyan = 2.83 ± 0.79 kpc).The stars of the Thamnos part (orange) with smaller   but larger   possess pericentric distances as large as another part (blue square), so pericentric distance usually increases with rotation velocity.

Highly retrograde I'itoi+Sequoia+Arjuna and Pontus
I'itoi (Naidu et al. 2020), Sequoia and Arjuna (Myeong et al. 2019) (ISA) are three highly retrograde and high-energy populations with distinct chemical abundances (i.e., [/Fe], [Fe/H]).Therefore, we take the only extremely retrograde and high-energy group as ISA.Although ISA and Thamnos overlap in (  −  )/ tot and   / tot space, the energy, retrograde rotation and apocentric distance are much larger than those of Thamnos.Its remarkable energy is consistent with that of Sequoia in Figure 2 of Koppelman et al. (2019b) as shown in the bottom-right panel of Figure 3 here.Based on three metallicity levels, Naidu et al. (2020) clearly defined the Arjuna stars as those of [Fe/H] > −1.5, the Sequoia stars of −2 < [Fe/H] < 1.5, and the I'itoi stars of [Fe/H] < −2.The bottom-left panel of Figure 3 shows that a few ISA members satisfy (  −   ) > 0, while Naidu et al. (2020) 2020) for details), the ISA group is composed of a lot of ultra-metal-poor stars from I'itoi and a small number of Sequoia stars as well as very few Arjuna stars.Here the ISA stars could be more retrograde than Sequoia in other works (e.g., Feuillet et al. 2021) because we pursue purity and discover its member stars through clustering, while previous studies wish to get high-completeness member stars by actions J cuts, which actually overlap with our ISA members in projected action space.In addition, the retrograde halo is more metal-poor, as advocated by Koppelman et al. (2019b).
Pontus recently discovered by Malhan et al. (2022) is a slightly retrograde merger, and it is indistinguishable from the GSE in (,   ) space or (  ,   ) space.Since stars with the same origin have a common age-metallicity relationship (Massari et al. 2019), Malhan et al. (2022) ruled it out as a fragment of the GSE by comparing its age-metallicity relationship with that of the GSE.Its dynamical properties are  = −172842 ± 3477 km 2 s −1 ,   = 385 ± 89 km s −1 kpc,   = 177 ± 214 km s −1 kpc,   = 423 ± 128 km s −1 kpc,  ⊥ = 896 ± 154 km s −1 kpc,  = 0.74 ± 0.10,  apo = 8.8 ± 0.7 kpc,  peri = 1.4 ± 0.6 kpc and  max = 7.2 ± 1.0 kpc. Figure 3 shows that Pontus is located at a distinct position from the GSE and Thamnos in the projected action space, which indicates that the members of the GSE move along more in-plane and radial orbits, but Thamnos is more retrograde compared to Pontus.

DISCUSSION AND CONCLUSIONS
Based on the precise parallaxes and proper motions from Gaia DR3 in combination with the radial velocities, [/Fe] and [Fe/H] from LAMOST DR9, we select the VMP stars (i.e., [Fe/H] < −1.8) to explore the dynamical relics of halo objects in the solar neighbourhood.Since small or low-mass systems are primarily populated by very metal-poor stars, the VMP sample can help us to trace fine structures that would be overwhelmed during clustering if all samples were considered.In this work, we apply the improved two-stage SNN clustering algorithm to identify substructures with similar actions and energy (J, ).Not only the observation errors are considered, but also the 10, 000 SNN clustering results are further clustered to obtain stable groups in order to overcome the parameter adjustment problem.
In order to obtain the characteristics of each merger in more detail, we further investigate the parts of distinct dynamical properties (i.e.,  and (  −   )/ tot ) for the GSE, Thamnos, and Wukong.In energy and angular momentum (,   ) space, the GSE stars appear slightly prograde at high energies but symmetric with respect to |  | = 0 at low energies.The  distribution is also rather wrinkly with several overdensities gathered in low |  |.Similarly, it corresponds to over-dense chevrons in the -  phase-space.The high-energy and medium-energy GSE parts, with apocenters in the ranges of 29.5 ± 3.6 kpc and 13.0±2.7 kpc, presumably contribute to the two breaks at 15−18 kpc and 30 kpc in the "double-break" density profile predicted by Naidu et al. (2021), respectively.We found that the difference between two slightly retrograde Thamnos parts, with energy in the same range of −165280 ± 4171 km 2 s −2 , is primarily due to their zenithal or vertical and azimuthal velocities, which results in two clumps (i.e., (  −   )/ tot < 0 and (  −   )/ tot > 0) in the projected action space.The difference between two Wukong parts in vertical action or velocity leads to their slightly distinct orbital parameters such as  max and  apo .The VMP members of MWTD move along the orbits with larger eccentricities than those of more metal-rich ([Fe/H] ∼ −1.12) the MWTD in Naidu et al. (2020), so we suggest that the VMP part of the MWTD may originate from accreted galaxies.Orbits of the highly prograde Helmi Streams in local samples are circular and rise to high vertical heights.We also analyze the dynamical properties of highly retrograde ISA with high energy, and Pontus hidden in the GSE, but we found that the members of Pontus move on more polar orbits than those of the GSE stars.In Table 2 we summarize all the mergers explored in this study that constitute the VMP halo in the solar neighbourhood.
Our work reveals the origin of local metal-poor stellar halos, which may provide significant insights into probing the origins of exotic stars, such as high-velocity stars.Study in the origin of high-velocity stars, such as the tidal debris of disrupted dwarf galaxies (Du et al. 2018a,b) and LMC (Erkal et al. 2019), is very interesting, recent results (e.g., Huang et al. 2021;Li et al. 2022Li et al. , 2023) ) show some candidate high-velocity stars originating from the Sgr dwarf spheroid galaxy.It is known that the energies of the Sgr member stars are so high that some member stars satisfy  > 0 or  >  escape (escape velocity).In the 26 groups we identified here, two special groups, namely high-energy GSE (hot pink circle) and ISA (red diamond), were discussed further due to their high velocities and energies, with their velocities mostly exceeding 300 km s −1 .As we expected, stars at very high energies are so sparse that they can hardly be identified as clumps in the action-energy space and the two high-energy groups are only about  = −1.2[×10 5 km 2 s −2 ].The retrograde ISA stars may come from the GSE's outer disc due to the host galaxy moving significantly in the merger.Therefore, we conjecture that some VMP high-velocity stars, especially those with large radial or extremely retrograde velocities, may originate from the massive GSE satellite, and its orbital angular momentum tilted by 30 • from that of the host disc when merging (Naidu et al. 2021).

Figure 1 .
Figure 1.The member stars in a massive substructure, namely GSE, are in the -  space (top panel) and   - space (bottom panel) and colour-coded according to energy.

Figure 2 .
Figure2. 26 groups assigned to seven substructures plotted by distinct shapes and colour where error bars in each object are the standard deviations of the xand y-axis quantities.For the subdivided substructures (e.g., GSE, Thamnos, Wukong), each part is represented in a different colour but with the same shape.All groups are shown in  apo - space (top-left), projected action space (top-right),  peri - max space (bottom-left), and -  space (bottom-right).
Figure3.Seven substructures VMP stars described in the (J, E) space, using the same colour and shape scheme as Figure2.Note that the retrograde GSE part (royal blue circle) coincides with the polar Thamnos (blue square) in the -  space, but does not happen in the projected action space (bottom-left).

Figure 4 .
Figure 4. 3D kinematics in cylindrical coordinates for all substructures, using the same colour and shape scheme as Figure2.
6 and ∼ −1.8 [×10 5 km 2 s −2 ], while the GSE stars are inclined to prograde (  < 0) between  ∼ −1.1 and ∼ −1.6 [×10 5 km 2 s −2 ].In the bottom-right of Figure2, for GSE (pink, hot pink, cyan, and royal blue circles), the (,   ) distribution appears stratified, corresponding to the energy wrinkles found byBelokurov et al. (2023).In their work, energy wrinkles are only distributed at prograde and retrograde edges, while a continuous, smooth and asymmetric peak appears at low |  |.Here, energy wrinkles are tightly presented at low |  |, probably because what we obtained are only high-stability and compact substructures, leaving out those relatively scattered GSE stars.

Figure 5 .Figure 6 .
Figure 5. Left: GSE, Thamnos and ISA in the   - space.Right: Wukong, Pontus, Helmi Streams and MWTD in the   - space.Note that we present all substructures in two separate panels solely to avoid visual clutter.

Table 1 .
Summary of VMP groups.
Group represents the name of cluster derived from the two-stage SNN.n represents the size of cluster. it is the mean and standard deviation of probability for stars associated with each group.The contents in parentheses indicate colours coded in Figures2-6.

Table 2 .
Summary of existing substructures in local VMP halo.Substructure (  sub )  [ /Fe]   apo  peri  max (   −   )/ tot   / tot   sub represents the number of stars in substructure.