A Walk on the Retrograde Side (WRS) project I. Tidying-up the retrograde halo with high-resolution spectroscopy

Relics of ancient accretion events experienced by the Milky Way are predominantly located within the stellar halo of our Galaxy. However, debris from di ff erent objects display overlapping distributions in dynamical spaces, making it extremely challenging to properly disentangle their contribution to the build-up of the Galaxy. To shed light on this chaotic context, we started a program aimed at the homogeneous chemical tagging of the local halo of the Milky Way, focusing on the component in retrograde motion, since this is expected to host a large fraction of stars accreted from past mergers. The A Walk on the Retrograde Side (WRS) project targets retrograde halo stars in the Solar Neighborhood having accurate 6-D phase space information available, measuring the precise chemical abundance of several chemical elements from high-resolution spectroscopy. In this first paper, we present the project and the analysis of high-resolution spectra obtained with UVES at VLT and PEPSI at LBT for 186 stars. Accurate radial velocity and chemical abundance of several elements have been obtained for all the target stars. In particular we focus on the chemical composition of a specific subset of substructures identified dynamically in the literature. Our study reveals that two among the more recently discovered structures in the retrograde halo, namely Antaeus / L-RL64 and ED-3, have identical chemical patterns and similar integrals of motion, suggesting a common origin. In turn, the abundance patterns of this unified system di ff er from that of Gaia-Enceladus, confirming that it is an independent structure. Finally, Sequoia exhibits a di ff erent chemistry with respect to that of Gaia-Enceladus at [Fe / H] < − 1 . 5 dex, showcasing an excess of stars with lower Mg and Ca in the common metallicity range.


Introduction
The Λ cold dark matter (ΛCDM) cosmological paradigm predicts that galaxies grow their masses hierarchically through merging events with smaller systems (White & Frenk 1991).The discovery of the ongoing accretion event of the Sagittarius dwarf galaxy (Ibata et al. 1994) has demonstrated for the first time that the Milky Way (MW) is not an exception.As extensively shown by N-body simulations of galaxy formation (Newton et al. 2018), the stellar halo of the surviving galaxy mainly comprises stars that have been accreted from smaller, disrupted satellite galaxies.Nevertheless, it has been established that the stellar halo is not exclusively composed of accreted objects (Carollo et al. 2007;Pillepich et al. 2015;Monachesi et al. 2019;Font et al. 2020;Belokurov & Kravtsov 2022;Khoperskov et al. 2023).Detecting stellar substructures in the Galactic halo yields valuable insights into the formation and evolution of the MW, since the overall configuration of the halo Full Tables 1, 3, 4, 6, and 7 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/684/A37has been deeply affected by the merger history of our Galaxy.Hence, comprehending the characteristics and mechanisms of formation of these substructures plays a pivotal role in reconstructing the evolutionary history of the MW.The search for substructures in the phase space is currently the most common and successful method used to unravel ongoing or recent accretion events (Ibata et al. 1994;Helmi et al. 1999;Bonaca et al. 2012;Helmi 2020;Belokurov et al. 2020).However, phase-mixing has proven some limitations to this approach, making it challenging to identify merging events that took place in the early stages of the MW's formation.Instead, N-body simulations have shown that debris from the same progenitor maintain a dynamical coherency in the space defined by the integrals of motion (IoM, e.g., energy and angular momentum) in a slowly evolving potential (Johnston et al. 1996;Helmi & de Zeeuw 2000;Gómez et al. 2013).In this context, the advent of the ESA/Gaia mission (Gaia Collaboration 2016) has provided a deeper understanding of the MW overall structure, shedding light on its formation and evolution, thanks to the extremely accurate measurements of the 3D position and velocity for ≈33.8 million of stars brighter than G = 14.0.According to the widely accepted scenario (see Helmi 2020 for a review), a substantial fraction of the MW stellar halo is believed to have formed as a consequence of a merging event with a dwarf galaxy, known as Gaia-Enceladus Sausage (GES, Belokurov et al. 2018;Helmi et al. 2018).
Recent studies revealed that the stellar halo contains an increasing number of other smaller substructures, including streams, clumps, and overdensities that are believed to be the remnants of past accretion events (Koppelman et al. 2019;Massari et al. 2019;Myeong et al. 2019Myeong et al. , 2022;;Kruijssen et al. 2020;Naidu et al. 2020;Yuan et al. 2020;Horta et al. 2021;Ibata et al. 2021;Malhan et al. 2022;Oria et al. 2022;Ruiz-Lara et al. 2022;Shank et al. 2022;Tenachi et al. 2022;Belokurov et al. 2023;Dodd et al. 2023;Mikkola et al. 2023;Cabrera Garcia et al. 2024).However, the interpretation of these numerous recent discoveries can be extremely challenging.In fact, numerical simulations of a major merger between a MW-like galaxy and a GES-like progenitor have revealed that remnants from the same merging event can result in overdensities at various locations in the space defined by the IoM, which could be mistaken for independent substructures (Koppelman et al. 2020;Amarante et al. 2022;Belokurov et al. 2023;Davies et al. 2023).On the other hand, relics with different origin may significantly overlap in the IoM space (Helmi & de Zeeuw 2000;Naidu et al. 2020;Lövdal et al. 2022;Dodd et al. 2023) and unrelated sources can drastically contaminate samples of candidate members of a given substructure selected kinematically (Buder et al. 2022;Rey et al. 2023).
A possible solution to this dynamical degeneracy can be achieved by adding an entirely independent dimension to the parameter space, that is, the chemical composition of stars (chemical tagging, Freeman & Bland-Hawthorn 2002).Indeed, the chemical abundance patterns of individual stars exhibit characteristics that reflect the star formation and chemical enrichment histories of the site where they were born.Therefore, they can greatly help in understanding the origin of a given substructure, namely, by distinguishing between accreted and in situ formation, and/or to infer the properties of the (now dissolved) progenitor satellite (Gallart et al. 2005;Tolstoy et al. 2009;Nissen & Schuster 2010;Mackereth et al. 2019;Kruijssen et al. 2020;Minelli et al. 2021;Montalbán et al. 2021;Mucciarelli et al. 2021b).
Large spectroscopic surveys, such as APOGEE (Abdurro'uf et al. 2022), GALAH (Buder et al. 2021), and H3 (Conroy et al. 2019), have made it possible to conduct in-depth studies in this field, providing detailed elemental abundances for an enormous number of stars.Many recent works have taken advantage of these large chemical datasets, focusing primarily on the chemical patterns of GES (Aguado et al. 2021;Bonifacio et al. 2021;Feuillet et al. 2021;Hasselquist et al. 2021;Limberg et al. 2022).Still, there are interesting regions of the parameter space where the sampling of chemical abundances is sparse and occasional, as, for example, the highly retrograde halo (RH) that is expected to host a very large fraction of accreted stars (Naidu et al. 2020;Myeong et al. 2022;Horta et al. 2023).Also, the merging and comparison between chemical abundances from different sources has been shown to be somehow problematic due to the lack of homogeneity in the atmospheric parameters and abundances scales of the various surveys (see, e.g., Helmi 2020).
Within this framework, we present a small spectroscopic survey aimed at delivering detailed chemical tagging of RH stars in the solar neighborhood with accurate kinematics, based on the precise measurements of several elements (e.g., Fe, Na, Mg, Al, Ca, Sc, Ti, Mn, Ni, Zn, Y).The A Walk on the Retrograde Side (WRS) project is aimed at obtaining chemical abundances from high-resolution spectroscopy for at least a few hundred stars, using accessible observational facilities.Spectroscopic observations of relatively bright stars (G 14.0) can be an excellent use of nights with sub-optimal seeing and lunar illumination at large telescopes.This is the niche that is currently being used for WRS observations.
The multi-dimensional dataset we plan to obtain will provide deep insight on the origin of each target star, by looking at similarities in the orbits and abundance patterns among them along with those of known relics of merging events.The approach we propose is complementary to that of largescale general surveys such as H3 (Naidu et al. 2020(Naidu et al. , 2022) ) and APOGEE (Horta et al. 2023).Indeed, the WRS project is highly focused, targeting only RH stars in the solar neighborhood, within D ≤ 1.0 kpc.Retrograde stars are particularly interesting in this context, as it is more likely that they come from the accretion of ancient satellites, with respect to their prograde counterparts (Carollo et al. 2007).Moreover, the WRS targets are selected among those having very accurate astrometry from the Gaia mission (Gaia Collaboration 2016).In addition to the constraint in distance and to the accurate radial velocity measure obtained within WRS, this ensures the most reliable determination of orbital parameters and IoM for the considered stars.In the end, the stellar spectra obtained for WRS have higher spectral resolution (R ≈ 40 000) than large spectroscopic surveys.This enhanced spectral resolution facilitates the precise determination of chemical abundances and offers finer insights into stellar characteristics.
In this first paper of the series, we describe our sample selection process and present the results of the analysis of a first batch of 186 stars.In particular, we conduct a high-precision and homogeneous abundance analysis on RH stars that are associated with different dynamical substructures, according to Dodd et al. (2023).We compare the results obtained for these substructures with the aim of uncovering their true nature and understanding whether they are all independently accreted components or whether some have originated from the same progenitor.
The paper is structured as follows.In Sect.2, we describe the selection criteria used to identify RH stars and we present the high-resolution spectroscopic dataset.In Sect.3, we present the dynamical characteristics of our sample of stars.In Sect.4, we provide a detailed description of the criteria that were employed to define each individual substructure in the RH.In Sect.5, we provide a comprehensive account of the detailed procedure employed to conduct the chemical analysis.In Sect.6, we show the comparison with literature results.In Sect.7, we showcase an extensive examination of the halo substructures in different chemical planes, focusing in particular on odd-Z elements (Sect.7.2), α-elements (Sect.7.3), iron-peak elements (Sect.7.4), and neutron capture elements (Sect.7.5).In Sect.8, we summarize and discuss the main results of our work.

Data selection
At the epoch when the WRS project was conceived, we had to face the problem of selecting retrograde stars from catalogues, like the Gaia DR2 (Gaia Collaboration 2018b) and A37, page 2 of 21 EDR3 (Gaia Collaboration 2021) databases 1 , where most of the stars were lacking a line of sight velocity (V los ) measurement in the range of apparent magnitudes of interest (G ≤ 14.0, for the time being).We tackled this issue by realizing that a selection in tangential velocity V T = V 2 ra + V 2 dec , that is, the linear velocity in the plane of the sky, can be used to effectively select very likely retrograde stars.
This technique is illustrated in Fig. 1, where the wellmeasured Gaia EDR3 stars within 1.0 kpc from the Sun with the complete 6D phase space information available are plotted in a Toomre diagram that is color-coded according to V T 2 .It is quite clear that in the vicinity of the Sun, high V T values are displayed only by stars on retrograde orbits or on slightly prograde orbits but high V 2 R + V 2 Z , the latter being equally interesting in our view, as they are typically associated with GES.Hence, we adopted V T > 400 km s −1 as our criterion to select likely retrograde stars.The criterion proved to be very effective as only 16 out of the 186 stars considered here and selected in this way are not retrograde, all of them having L z < 800 kpc km s −1 , while 123 of them have L z ≤ −500 kpc km s −1 .In future observational campaigns, we plan to select retrograde targets directly in V Φ among stars having V los measurements in the Gaia DR3 catalogue or catalogues from future Gaia data releases.This will allow us to sample also the mildly retrograde stars at relatively low total velocity (−100 km s −1 V Φ < 0.0 km s −1 and ) that are now excluded by our selection in V T , as well as to target specific substructures.
1 That have the information on the V los for 7.2 × 10 6 stars with G 13.0, while the most recent release, DR3 (Gaia Collaboration 2023), includes V los for 33.8 × 10 6 stars down to G 14.0. 2 The adopted solar motion is from Schönrich et al. (2010), while the distance between the Sun and the Galactic Center and the rotational velocity of the Sun are the same used in McMillan (2017).
The additional selection criteria we adopted were intended to select well-measured stars in the solar neighborhood.In particular we selected our targets among V T > 400 km s −1 stars having G 14.0.This suitable magnitude range was selected in order to get high signal to noise ratio (S/N), high-resolution spectra with 8 m class telescopes with relatively short exposures (typically ≤3000 s).Then, we imposed an error on the parallax lower than 10%.This implies small errors in the distance and in the proper motion.We selected stars with a distance from the Sun of D ≤ 1.0 kpc.We adopted the following criteria in order to ensure that the Gaia astrometry and photometry of the target stars are of the highest quality: RUWE ≤ 1.3 (Gaia Collaboration 2021) and C < 3σ (Riello et al. 2021).The first criterion is also effective in excluding astrometric binaries (see Andrew et al. 2022, and references therein).Finally, by inspecting the colormagnitude diagram, we excluded by hand a few dwarf stars that may be unresolved binaries, which are located slightly above the main sequence.
The most relevant outcome of the selections above, complemented by uncertainties in the V los measured from high-resolution spectra lower than 0.9 km s −1 , is that the typical uncertainty on the IoM is just a few percent (see Sect. 3), hence, the noise in the phase space is minimal in our sample.Basic information on our target stars is listed in Table 1.

High-resolution spectroscopy for selected stars
At the present stage, we have collected high-resolution spectra for 98 stars (Program IT-2021B-004 and Program IT-2022B-006, P.I.: M. Bellazzini) with the optical spectrograph PEPSI (Strassmeier et al. 2015) mounted at the Large Binocluar Telescope in Arizona.The observations have been performed with the CD3 grating for the Blue Arm, which cover the spectral range between 4800-5440 Å, and the CD6 grating for the Red Arm with a spectral coverage between 7410-9140 Å, both with a spectral resolution of 40 000.The choice of this set-up for the two arms was dictated by the trade-off between the access to measurable lines, allowing to derive the abundance of scientifically relevant chemical elements, and the efficiency of the spectrograph with the goal of achieving S /N > 30 for a reliable measure of the chemical abundances we are interested in.In order to maximize the number of observed stars, we aimed for short exposures.This translates into a magnitude limit of G < 14.0.CD3 is less efficient than CD6, therefore it drives the exposure times.The short exposures allow us to reach typically S /N 40 for the CD3 grating and S /N 60 for CD6.The spectra have been reduced with the dedicated PEPSI pipeline, which contains bias subtraction, flat-fielding, spectral extraction, wavelength calibration and normalization.The observations of close sky region in the same exposure times of the targets have been executed in order to subtract the sky contribution to the stellar spectra.
We also observed 100 stars (Program 0109.B-0522, P.I.: A. Mucciarelli) with the optical spectrograph UVES (Dekker et al. 2000) mounted at the Very Large Telescope of the European Southern Observatory.Observations were performed with UVES in Dichroic mode adopting the standard settings Dic 1 Blue Arm CD2 390 (3600-4800 Å) and Dic 1 Red Arm CD3 564 (4800-6600 Å) and with the 1 × 12 slit, thus yielding a resolution of R = 40 000.We are able to obtain, on average, S /N 20 for the Blue Arm and S /N 40 for the Red Arm.All the observed spectra have been reduced with the ESO pipeline3 .The spectra of 12 targets were collected with both instruments as a reference point for the purpose of calibration and direct comparison.Detailed information on some representative stellar spectra are listed in Table 2, bracketing the range of exposure times adopted and reporting the achieved S/N in several spectral region of interest.As an example of the quality of the spectroscopic dataset, we plot in Fig. 2 the spectra of four target stars with very similar atmospheric parameters and a difference in [Fe/H] of ≈1.2 dex.

Orbital parameters
Position, proper motion, and parallax measurements were taken from Gaia EDR3 (Gaia Collaboration 2021).We corrected each star's parallax for zero-point offsets following the prescriptions described by Lindegren et al. (2021).Also, V los was measured from the high-resolution spectra by cross-correlation with template spectra using the fxcor task of IRAF4 .For each star, as the template spectrum, we exploited a synthetic spectrum calculated using the correct atmospheric parameters (see Sect. 5.2).The V los is measured at the maximum peak of the cross-correlation function.Uncertainties on the line-of-sight velocity were computed following the method described in Tonry & Davis (1979).
As shown in the central panel of Fig. 3, they are of the order of 0.1−0.9km s −1 , which significantly reduce the uncertainties arising from the Gaia measurements, especially for faint stars.Furthermore, we corrected the measured V los for the gravitational redshift effect following the prescription described by Zwitter et al. (2018).Since most of our stars are dwarfs and the precision of our V los measures is so high, this correction (albeit small) is non-negligible, ranging from 0.1−0.7 km s −1 .The distribution of the final V los is displayed in the left panel of Fig. 3 and the values are listed in Table 3.The difference between our measurements and those from Gaia DR3 (Gaia Collaboration 2023; Katz et al. 2023), defined as ∆V los = V los,spectrum − V los,GaiaDR3 , is consistent within 3σ for all except five stars, where σ is calculated as the squared sum of the uncertainties arising from high-resolution spectroscopy and Gaia DR3 (see right panel of Fig. 3).The discrepancy observed in these five stars could hint at a binary nature for these objects.This implies the possibility of their having undergone a distinct and unconventional chemical evolution or either having the spectrum contaminated by the companion star.Consequently, we flag these stars with a distinct marking in Table 3, and they will be omitted from the following chemical analysis.The final V los values have been coupled with the parallax and proper motions provided by the EDR3 to draw the complete 6D phase space information for each target.
We then transformed the observed kinematic information to the Galactocentric reference frame.To do so, we assume the Sun to be located R 0 = 8.122 kpc away from the Galactic Center (GRAVITY Collaboration 2018) and z = 20.8pc above the midplane (Bennett & Bovy 2019).We obtain the solar  velocity by combining the definition of the local standard of rest by Schönrich et al. (2010) with the proper motion of Sgr A * (Reid & Brunthaler 2004).In this frame, the Sun moves with (U , V , W ) = (12.9,245.6, 7.78) km s −1 (Drimmel & Poggio 2018).
We computed stars orbital parameters with the software AGAMA (Vasiliev 2019) using a McMillan (2017) potential for the MW.We ran 100 Monte Carlo simulations of the orbit for each star assuming Gaussian distributions for the uncertainties in distance, proper motion, and radial velocity.The final values of the dynamical parameters are computed as the median of their derived distributions, with associated uncertainties at the 16th and 84th percentiles.The typical error in E, L z and L ⊥ are 0.8%, 2.6%, and 1.8%, respectively.The distribution of our sample in these dynamical spaces is shown in Sect. 4.

Identification of substructures in the RH
Recent studies have revealed that the RH exhibits plenty of substructures that significantly contribute to the total mass budget of the stellar halo (Naidu et al. 2020).In most cases, they have been identified as overdensities in the space of IoM, such as their actions, the total energy (E), the angular momentum along the z-direction (L z ), and the angular momentum perpendicular to L z (L ⊥ ).We note that L ⊥ has been extensively used to identify stars formed in the same progenitor, even if it is not completely A37, page 5 of 21 Table 3. Kinematic information for the stars: ID from Gaia EDR3, orbital energy, angular momentum along the z-axis, perpendicular angular momentum, and line-of-sight velocity with uncertainties.conserved in an axisymmetric potential such as the one we assumed for our galaxy (Helmi et al. 1999;Helmi & de Zeeuw 2000;Bonaca et al. 2021).
As a first application of this first bunch of WRS data, we attempted to carry out a chemical characterization of some of the substructures recently listed and discussed by Dodd et al. (2023, D23 hereafter), who built on the findings by Lövdal et al. (2022) and Ruiz-Lara et al. (2022).These authors used a specifically designed clustering algorithm based on the IoM computed using Gaia DR3 data to identify groups of stars sharing similar properties in the phase space.
In the region of the E, L z , and L ⊥ space covered by our sample, the most widely represented among the D23 groups is GES.We sample the most retrograde and high energy slice of GES, based on our selection window.In addition, there are two mid-sized slightly retrograde substructures associated to Sequoia (Myeong et al. 2019) and Thamnos (Koppelman et al. 2019), as well as three newly found, small, high-energy, and extremely retrograde clumps, namely L-RL64, also known as Antaeus (Oria et al. 2022), ED-2, and ED-3, plus other minor groups that do not appear to have significant overlap with our sample (see more details below).Indeed, we find that 179 of our 186 stars are included in the D23 sample.55 of them are associated with one of the D23 groups.However, there are five groups that have just one member star included in our sample, namely Thamnos, ED-5, ED-6, Group 13, and Group 14.Since no meaningful conclusion can be drawn from an individual star, we do not consider these substructures in the following.Still, we do make the chemical abundances available since they might be of interest for upcoming investigations.Of the remaining 50 stars, 12 are associated to GES, 18 to Sequoia, 11 to Antaeus, 4 to ED-2, and 5 to ED-3.
Starting from these 50 "core members" and with the aim to increase the number of stars associated to a substructure, we performed a clustering analysis with the algorithm DBSCAN (Ester et al. 1996).We used E, L z , and L ⊥ as the input parameters for the clustering algorithm, for the purpose of having a fair comparison with D23.We scaled the space of these parameters with the RobustScaler tool in scikit-learn (Pedregosa et al. 2011).Due to the significant overlap between halo substructures in the IoM spaces, it is difficult for any clustering algorithm to identify each individual overdensity.Furthermore, some of the newly discovered substructures include a low number of stars.For this reasons, the algorithm used in this study was fine-tuned to detect all the small-scale clumps identified by D23.We find the best matching results with D23 by using DBSCAN default hyperparameters and by setting eps = 0.34 and min_samples = 5. Figure 4 shows the position of the final candidate members of each substructure in the IoM space.With our final selection, we were able to relate 119 stars in our parent sample to the known retrograde substructures.In the following sections, we discuss the sample defined for each substructure in detail.
Before starting the discussion of the various substructures, it is important to note that the identification and classification of completely disrupted relics of ancient merging events is still in its infancy and, as made clear in the following discussion, it is subject to considerable uncertainty.On the one hand, there is no consensus on the best selection criteria to isolate the purest samples of a given substructure and any criterion will imply some degree of contamination from unrelated sources (see, e.g., Buder et al. 2022, for references and discussion).On the other hand, as stated above, there are substructures that may be erroneously considered as an individual entity when they may, in fact, be the sum of independent substructures or the unrecognized part of a larger one (Koppelman et al. 2020;Naidu et al. 2020;Belokurov et al. 2023;Dodd et al. 2023).Here, we aim to get a clearer insight into the nature of some of the retrograde substructures by adding the chemical dimensions, as done by other authors (e.g., Feuillet et al. 2020Feuillet et al. , 2021;;Buder et al. 2022;Horta et al. 2023).Overall, this is not an easy task in the present phase and it is worthwhile to proceed cautiously, taking care of every step along our exploratory path.

Gaia-Enceladus Sausage
The dominant accreted component in the nearby stellar halo is the remnant of the last major merger event experienced by the MW with an ancient (8-11 Gyr ago) and relatively massive (M ∼ 10 8 −10 9 M ) dwarf galaxy, known as GES (Belokurov et al. 2018;Helmi et al. 2018;Gallart et al. 2019;Mackereth et al. 2019;Massari et al. 2019;Vincenzo et al. 2019;Kruijssen et al. 2020;Hasselquist et al. 2021).One of the defining characteristics of this population is the stars' highly eccentric orbits and slightly retrograde motion.The metallicity distribution function of GES peaks from [Fe/H] −1.5 to −1.1 dex, depending on the sample selection and the spectroscopic survey (Naidu et al. 2020; A37, page 6 of 21  Due to the criterion employed to select retrograde stars (see Sect. 2), which significantly cuts the IoM region typically populated by GES stars, we were able to retain only 12 stars in common with D23 belonging to GES.By using the clustering algorithm, all of the original members of this cluster were retained as GES members and an additional 41 stars were added to it, making the total number of stars in this cluster equal to 53 (see the green-colored symbols in Fig. 4 and thereafter).

Sequoia
The first evidence of the remnant of an accreted dwarf galaxy in the RH was the identification of a clump of stars moving on highly retrograde orbits, with large orbital energy and forming an arc-like structure in the Toomre diagram, now dubbed Sequoia (Myeong et al. 2018(Myeong et al. , 2019;;Koppelman et al. 2018;Matsuno et al. 2019).Sequoia has been proposed to be the remnant of an accreted independent small dwarf galaxy due to its position and smaller extent in the IoM space with respect to GES (Koppelman et al. 2019;Myeong et al. 2019).However, the actual complexity of this substructure is far greater than what was initially perceived.A different scenario has been proposed to explain the Sequoia overdensity, possibly making it tightly related to the GES accretion event.In fact, it has been shown that a GES-like progenitor, depending on its orbital configuration and morphology, is able to display its debris on orbits that reproduce the same arc-like feature asso-ciated with Sequoia (Koppelman et al. 2020;Amarante et al. 2022;Belokurov et al. 2023).The chemical composition of Sequoia remains a subject of ongoing debate, as different studies have yielded conflicting results regarding its abundance of α-elements compared to the GES substructure.Several works have reported either an enhancement (Mackereth et al. 2019;Myeong et al. 2019) or a depletion (Monty et al. 2020;Matsuno et al. 2022) of α-elements compared to GES abundances.Moreover, Naidu et al. (2020) found a multiple-peaked metallicity distribution for Sequoia, and claimed that two additional sub-components, namely Arjuna and I'Itoi, occupy the same region in the IoM space (see also, Bellazzini et al. 2023).In addition, Ruiz-Lara et al. ( 2022) discovered three disconnected clusters in the region of the IoM space previously associated to Sequoia that might be related to different progenitors.
The cross-match provides 18 common members between our targets and the Sequoia sample of stars identified by D23.When using the DBSCAN algorithm, only 8 stars were found to belong to this clump; additionally, 2 stars from the original sample were associated with ED-2 (see Sect. 4.4).This association was primarily based on the alignment of these two stars with the ED-2 clump in the L z − L ⊥ space.However, they were found to move on higher energy orbits with respect to ED-2.Consequently, this distinctive orbital behavior led us to consider the possibility of these stars being affiliated with the Sequoia substructure.In the end, since all the 8 stars are part of the 18 original members, we decided to stick with the original D23 classification for Sequoia.Thus, our final Sequoia sample contains 18 stars (brown-colored symbols in Fig. 4).

Antaeus/L-RL64
L-RL64 is a newly discovered cluster located in the most retrograde region of the IoM space (Lövdal et al. 2022;Ruiz-Lara et al. 2022).L-RL64 is one of the three clusters A37, page 7 of 21 associated to the Sequoia merging event and it was initially believed to be its extreme retrograde tail.Subsequently, Oria et al. (2022) independently rediscovered this substructure in the action space, giving it the name Antaeus.The peculiar kinematics of this clump (which we henceforth refer to as Antaeus) combined with the mass estimate derived from its mean metallicity, indicate a different origin for this substructure with respect to Sequoia.Indeed, it is now thought that Antaeus is an independent accretion event, not related to Sequoia (Oria et al. 2022;Ruiz-Lara et al. 2022;Dodd et al. 2023).
The analysis based on D23 data shows that 11 stars in our entire sample are associated with Antaeus.However, the clustering algorithm identifies two separate clumps composed of 15 stars moving on extremely retrograde orbits.Here, we treat them as one single substructure (blue-colored symbols in Fig. 4), since their position in the IoM spaces is totally consistent with the one of Antaeus and 13 out of these 15 stars are labeled by Oria et al. (2022) as members of this group.Thus, with our final selection we define a fair sample of this newly discovered substructure and we will be ableto provide its detailed chemical characterization for the first time.

ED-2 and ED-3
D23 spotted two new small highly retrograde clusters near Sequoia, namely ED-2 and ED-3, which are extremely tight in the IoM and velocity spaces (see their Fig. 3).We thus have a total of 4 and 5 stars in common with the ED-2 and ED-3 clusters, respectively.The results provided by DBSCAN show that we are able to retain all the 5 stars in the ED-3 clump and also add 1 star to the ED-2 substructure (orange-and magenta-colored symbols in Fig. 4, respectively).Notably, as remarked above, the clustering algorithm spotted 2 additional stars that were previously associated with Sequoia as members of ED-2.4.5.Thamnos, Clump-6, and Clump-7 Koppelman et al. (2019) revealed the presence of a less prominent substructure with typically lower orbital energy in the retrograde region of the E − L z space with respect to GES, which they named Thamnos.Given its small extent in the E − L z plane, they estimate a stellar mass <5 × 10 6 M .Thamnos has been suggested to be the remnant of an old merging event with the MW, due to its low mass and peculiar location in the IoM (Kruijssen et al. 2020), but it is still unclear if the stars associated to this group are from a unique remnant or two (Koppelman et al. 2019;Bellazzini et al. 2023).
In our entire sample, only 1 star is linked to Thamnos according to D23.As is in the case with GES, the selection method used does not allow us to include stars with orbital energy as low as Thamnos, which means we are unable to sample this particular substructure.Interestingly, DBSCAN spots two clusters composed of 17 and 6 stars in the low energy regime of our sample at slightly higher orbital energy with respect to Thamnos.These two groups have very similar E and L z , but they differ in L ⊥ .It is worth noting that these clumps may be artificial overdensities identified by the algorithm and produced by our selection cut.In the following, we refer to them as Clump-6 and Clump-7 (cyanand red-colored symbols in Fig. 4).

Unassociated stars
A significant fraction of our sample, 67 out of 186 stars, does not appear to be associated with any known structure according to D23 and the clustering algorithm.In the rest of the paper, we do not focus on the chemical features of these unassociated stars, as they will be part of future contributions to the WRS project, exploiting their chemistry to determine whether they are accreted.

Stellar parameters
Establishing precise atmospheric parameters (i.e., effective temperature, T eff , surface gravity, log g, and microturbulent velocity, v t ) is a crucial step to obtaining reliable chemical abundances.
We derived both T eff and log g relying on the EDR3 photometry (see Table 1).In particular, T eff has been calculated with the (BP − RP) 0 −T eff relation provided by Mucciarelli et al. (2021a).The color (BP-RP) has been corrected for extinction following the iterative procedure described in Gaia Collaboration (2018a).The color excess E(B − V) adopted for each star is provided by the online tool EXPLORE5 .Since the adopted color -T eff relation depends on the star metallicity ([Fe/H]), which is not known a priori, we calculated T eff while assuming a fixed value of [Fe/H] = −1.5 dex for any star.The uncertainty on T eff is due to the errors in photometric data, parallax, color excess, and (BP − RP) 0 − T eff relation.Given that the observed targets are bright, very nearby and affected by minimal extinction, the main source of error in the effective temperature arises from the color and T eff relation.Thus, we assumed the dispersion of the (BP − RP) 0 − T eff relation, which is 61 K (Mucciarelli et al. 2021a), as the typical internal error on T eff .Surface gravities have been estimated assuming the photometric T eff , a stellar mass derived from a BaSTI isochrone (Pietrinferni et al. 2021) with 12 Gyr, [Fe/H] = −1.5 dex and [α/Fe] = +0.4dex (according to the typical values for the accreted component of the stellar halo derived in previous studies, Nissen & Schuster 2010, Kruijssen et al. 2020, Montalbán et al. 2021, Horta et al. 2023) and the G-band bolometric corrections described by Andrae et al. (2018).We note that employing an alternative isochrone varying the parameters has a minimal impact on the resulting surface gravity of the star.For instance, when using an isochrone with 10 Gyr age, [Fe/H] = −1.1 dex and [α/Fe] = +0.4dex, the discrepancy with the first results remains consistently below 0.04 dex.The uncertainties on log g are derived through the propagation of the errors on the effective temperature, photometry, and distance.The obtained values are systematically lower than 0.1 dex.Finally, microturbulent velocities were calculated by requiring no trend between iron abundances and reduced equivalent widths, defined as log 10 (EW/λ).The uncertainties on v t are estimated as described in Mucciarelli et al. (2013) and they are of the order of 0.2 km s −1 .
Subsequently, we improved the estimate of the stellar parameters by using the correct [Fe/H] value derived from the chemical analysis.The final effective temperature value deviates from the initial estimate by a margin persistently lower than 60 K.We excluded from the chemical analysis the three hottest stars of the sample (T eff > 6650 K).This decision stems from the fact that the combination of their atmospheric parameters leads to the formation of extremely weak absorption lines.As a consequence, this renders the measurement of abundances unreliable.In Fig. 5, we show the color-magnitude diagram (CMD) for all the analyzed stars, color-coded based on the metallicity.The final values of the atmospheric parameters are listed in Table 4.

Line selection
To obtain the chemical abundances of each star, we made a list of unblended lines to analyze, automatically pre-selected with our own code AUTOKUR and consequently refined with a visual comparison between the observed and the synthetic spectra.The latter were calculated with the code SYNTHE (Kurucz 2005), adopting the stellar parameters of the observed targets (see Sect. 5.1), assuming ATLAS9 (Kurucz 2005) model atmospheres, and using the linelists containing atomic and molecular data from the Kurucz/Castelli 6 database.The synthetic spectra have been convolved with a Gaussian profile to reproduce the observed resolution.
To address the fact that the degree of blending for a particular transition is dependent on the metallicity, which is not known a priori in this case, we utilized an iterative approach to establish the list of spectral lines to analyze for each star.The synthetic spectra have been calculated adopting an αenhanced mixture ([α/Fe] = +0.4dex) with an initial metallicity of [Fe/H] = −1.5 dex.Following the initial chemical analysis, a revised set of unblended spectral lines has been identified for each star using a synthetic spectrum that was computed using the appropriate chemical composition.
During the visual inspection, we noticed that the spectra of three low-main sequence stars (T eff < 5250 K and log g > 4.60) were highly contaminated by molecular transitions from MgH.Thus, we opted to exclude them from the abundance analysis.In the end, when considering also the exclusion of the 5 binary stars, a total of 11 stars were removed from the initial sample.We chose to exclude them from further consideration in our study.Among these 11 excluded stars, 4 were associated with GES, 1 with ED-2, 1 with Sequoia, 1 with Clump-6 and 4 stars remained unassociated with any specific substructure.

Chemical analysis
The chemical abundances of Na, Ca, Ti, Fe, Ni, and Zn were derived through a comparison between the observed equivalent widths (EWs) and the theoretical line strengths using the code 6 http://wwwuser.oats.inaf.it/castelli/linelists.html GALA (Mucciarelli et al. 2013).The EWs were measured with the code DAOSPEC (Stetson & Pancino 2008) via the automatic tool 4DAO (Mucciarelli 2013).
The chemical abundances of elements which show saturated core lines, but with wings that are still sensitive to the abundance or hyperfine/isotopic splitting transitions (Mg, Al, Sc, Mn, and Y) were derived through spectral synthesis with the proprietary code SALVADOR.This program runs a χ 2 -minimization between the observed line and a grid of synthetic spectra computed onthe-fly by the code SYNTHE, keeping the stellar parameters fixed and varying only the abundance of the matching element.
All the abundances were scaled to the solar composition, taking as reference the values from Grevesse & Sauval (1998).We adopted these solar abundances since the ATLAS9 model atmospheres used in this study are computed based on a chemical mixture consistent with the one outlined in Grevesse & Sauval (1998;see Castelli & Kurucz 2003).For all the elements that were examined, in cases where the spectral lines were not clearly detectable, we determined upper limits by comparing the observed spectra to the synthetic spectra.
To put the chemical abundances measured for the stars in the UVES and PEPSI samples onto the same scale, we made a direct comparison between the final abundances ratios of the 12 stars that have been observed with both instruments.For each element, we determined a median abundance difference (∆ abu = [X/Fe] PEPSI − [X/Fe] UVES ) and we applied an offset correction by adopting UVES as the reference scale (due to the largest number of available lines) in the cases where: and where σ is the standard deviation of the abundance differences for a given element and N is the number of stars with a measured abundance with both instruments for that specific element.For Fe, Ca, Sc, Ni, Zn, and Y, we corrected the abundances for the PEPSI sample by subtracting the offsets listed in Table 5, thus bringing them onto the UVES scale.We focus in the next paragraphs on the methodology we employed specifically to derive chemical abundances for problematic species.
Sodium.The available two lines of Na in the spectral range of PEPSI are centered at 8183.3 and 8194.8Å.We highlight that these atomic lines are situated within a spectral region significantly affected by telluric absorption lines.As a precautionary measure to ensure the robustness of our analysis, we visually inspected the possible telluric contamination of the Na lines using appropriate synthetic spectra of the Earth atmosphere calculated with TAPAS (Bertaux et al. 2014).Subsequently, we excluded any Na line exhibiting signs of contamination attributable to telluric features.We applied the correction for departures from LTE conditions from Lind et al. (2011).All the corrections for Na are available through the INSPECT database7 .The values listed in Table 6 are corrected for NLTE.
Aluminum.Within the UVES spectral range, the presence of Al is characterized by solely two strong lines located at 3944.0 and 3961.5 Å.However, these lines become saturated at [Fe/H] > −1.3 dex.In contrast, the PEPSI sample features only three weak aluminum lines suitable for the chemical analysis: one line at 7836.1 Å and the doublet at 8772-73 Å.Unfortunately, these weak lines are not detected at metallicities lower than ∼−1.5 dex.The corrections for non-LTE effects for each line are taken from Lind et al. (2022).Notes.These values were employed to ensure a consistent scale alignment between the two sets of measurements.

Abundance uncertainties
When inferring the uncertainties on the abundance ratios, it is essential to consider two primary sources of error: internal errors related to EW measurements and errors arising from the adopted atmospheric parameters.Uncertainties due to EW measurements have been estimated as the dispersion around the mean of the individual line measurements, divided by the root mean square of the number of used lines.For elements whose abundances have been measured using spectral synthesis, we exploited Monte Carlo simulations to estimate the internal error.Specifically, we generated a sample of 500 synthetic spectra with Poissonian noise to reproduce the S/N of the observed spectra and then we repeated the analysis to derive the abundance.The internal error is estimated as the standard deviation of the abundance distribution of the 500 noisy synthetic spectra.
To determine the uncertainties arising from atmospheric parameters, we replicated the chemical analysis by varying only one stellar parameter at a time and keeping the others fixed (see Sect. 5.1 for the details about the derivation of the uncertainties on the atmospheric parameters).Since we refer to chemical abundances with abundance ratios ([X/Fe]), the uncertainties on the iron abundance [Fe/H] have been considered as well.
The intrinsic and systematic components have been summed in quadrature to compute the total uncertainty.Therefore, the final errors have been estimated as: where σ X,Fe is the dispersion around the mean of chemical abundances, N X,Fe is the number of used lines and δ i X,Fe are the abundance differences obtained by varying the parameter i.All the abundance ratios for individual stars are listed in Tables 7 and 6, together with the corresponding uncertainties, as described in Sect.5.4.

Comparison with the literature
To put our results in a more general context, in this section we compare chemical abundances of our sample with external datasets, starting with the sample of metal-poor dwarf and subgiant stars associated with the halo population analyzed by Nissen & Schuster (2010;NS10, hereafter).This choice was guided by the extremely high-quality stellar spectra used by NS10 (S/N ranging from 250 to 500) that guarantees accurate measurements of chemical abundances.By using a sample of stars observed with UVES at VLT and FIES at NOT, these authors revealed for the first time the existence of two distinct populations within the nearby MW halo, composed of stars having different kinematics and chemical patterns.NS10 proposed that the Galactic halo population ([Fe/H] −0.8 dex) following the low-α sequence in the classical [Fe/H] vs. [α/Fe] plane (see the red squares in Fig. 6) is composed of stars accreted from disrupted dwarf galaxies, whereas the high-α population (green triangles in Fig. 6) is made up of stars that have formed in situ, a view that that is generally accepted nowadays (Bonaca et al. 2017;Hayes et al. 2018;Haywood et al. 2018;Helmi et al. 2018;Mackereth et al. 2019;Helmi 2020).To avoid the introduction of systematic errors due to the different chemical analysis and ensure a fair comparison between homogeneous A37, page 10 of 21   results, we defined our reference targets based on a sub-sample of 41 NS10 stars that have been observed with UVES using the same instrument setup as our targets and covering a range of atmospheric parameters similar to our sample.We re-analyzed the NS10 spectra following the same procedure described in Sect. 5 and using the same linelist.We found a median discrepancy of [Fe/H] = −0.05dex and σ [Fe/H] = 0.08; given that these values adhere to the criteria described in Eq. ( 1), we made the decision to adjust the abundances to ensure congruence within the NS10 reference frame.The estimates of [Mg/Fe] display a median difference of 0.01 dex and σ [Mg/Fe] = 0.08.We show the final comparison in Fig. 6: as expected from the selection we imposed to our targets, which should predominantly include retrograde, accreted objects, stars from our sample (blue symbols) are located in the same region of the α-elements chemical planes as the low-α population from NS10, with abundances that are systematically lower (≈0.2 dex on average) than the high-α population composed of in situ stars.It is also interesting to note that the low-α sequence of NS10 targets basically defines the upper limit of the distribution of our targets, which are mostly at lower [α/Fe] abundances.This might reflect the fact that we are sampling regions of the RH populated also by several different low-mass mergers, for which theoretical models predict low (and different) [α/Fe] levels already at fairly low metallicity (Matsuno et al. 2022).On the other hand, the accreted population of NS10 likely belongs to the most-massive merger, GES.
As a further and independent check of our results, we also made a comparison with chemical abundances from the large and widely used APOGEE DR17 (Abdurro'uf et al. 2022) dataset.We limited the comparison to well-measured APOGEE abundances by selecting only stars having ASPCAPFLAG = 0, A37, page 11 of 21 Interestingly, the stars from our sample that move on more extreme retrograde orbits (L z < −2500 kpc km s −1 ) appear to lay on a different pattern with respect to the low-α population at [Fe/H] < −1.5 dex, with lower [Mg/Fe] at fixed metallicity.In the right panel of Fig. 7, we show the [Mg/Mn] vs. [Al/Fe] chemical space, which serves as a diagnostic to distinguish between accreted and in situ populations (Das et al. 2020;Horta et al. 2021).Contrary to the APOGEE control sample, nearly all of our targets are located in the area of the plot where we expect to find accreted stars (Horta et al. 2021), emphasizing the fact that we are observing a poorly explored area of the solar neighborhood.Indeed, only less than 1% of the stars within the APOGEE control sample are located at D < 1 kpc and move on retrograde orbits.It is worth noting that only 11 stars among our sample are listed in the entire APOGEE DR17 catalogue, with only one of them aligning with the selection criteria expounded above.Moreover, through our program, we have the capability to provide detailed chemical abundance data for elements that are notably absent in APOGEE (e.g., Sc, Zn, and Y) and for elements whose measurements in APOGEE may diverge from theoretical expectations (Na and Ti).The comparisons presented above show that the selection we adopted was very effective in selecting stars that very likely originated in ancient and now-dissolved MW satellites.In addition to the location of our stars in the [Mg/Mn] vs. [Al/Fe] plane, it is worth noting that: (a) our kinematic selection picked up only metal-poor stars ([Fe/H] −0.8 dex), the majority of which (65%) having [Fe/H] ≤ −1.5 dex, with a significant very metalpoor tail (47 stars with −2.0 ≤ [Fe/H] < −3.0 dex, 26% of the sample), a component that is poorly probed by the APOGEE control sample; and (b) despite the higher resolution, our stars show a larger scatter in [Mg/Fe] with respect to APOGEE, increasing with decreasing metallicity, suggesting again that we are more effective in sampling stars from a variety of disrupted satellites, with different chemical evolutionary paths.In summary, the results of the adopted selection are fully compliant with the main rationale of WRS.
In Fig. 8, we compare the IoM of our targets with those of the NS10 sample and APOGEE.We computed the IoM for the NS10 stars using the 6D phase space information provided by the Gaia EDR3.Instead, for the APOGEE control sample we used the phase space values from Gaia listed in the APOGEE DR17 catalogue.First, Fig. 8  of our selection criterion in picking up halo stars in retrograde motion, as discussed in Sect. 2. However, this criterion does introduce a clear bias, since it prevents us from picking stars with low orbital energy (E ≤ 1.5 × 10 5 km 2 s −2 ) moving on slightly retrograde orbits (L z ≥ −0.5 × 10 3 kpc km s −1 ).At the same time, it is clear that high-quality APOGEE measures primarily targets disc stars with a quite sparse sampling of the retrograde component of the halo.For what it concerns the NS10 sample, although some high-α stars are interestingly located in the retrograde region (they might be part of the dynamically heated disc, sometimes referred to as the Splash, Belokurov et al. 2020), the bulk of this population moves on disc-like orbits, whereas the low-α moves on high eccentricity orbits, as expected for GES stars.

Chemical abundances of RH substructures
Here, we present the chemical abundances of the substructures in the RH of the MW.Having inspected the [Fe/H] vs. [X/Fe] diagrams for all the elements we measured, we discuss only those that appear to provide the clearest view of the substructures we consider here.In particular, we focus on the abundances of Na, Mg, Al, Ca, Sc, TiI, TiII, Mn, Ni, Zn, and Y.It is of utmost importance to acknowledge that due to impact of the combination of the atmospheric parameters onto the shape of absorption lines, our capacity to measure the chemical abundances for every depicted chemical element across all stars is restricted.Consequently, certain stars may not feature in specific figures.

Metallicity distribution functions
We go on to examine the metallicity distribution functions (MDFs) of the six substructures that were built exploiting only the IoM.In Fig. 9 we display the results, which suggest that the global MDF of our sample is bimodal, with a primary peak at [Fe/H] −1.4 dex and a secondary peak at [Fe/H] −2.4 dex.The high-metallicity end of the RH stars distribution is dominated by GES, with a median metallicity of [Fe/H] = −1.47dex and σ [Fe/H] = 0.32.Our result is consistent with Bonifacio et al.
Once again, we wish to emphasize the importance of our adopted selection criterion for GES stars.This criterion is structured in such a way that prevents us from including stars from the solar neighborhood with [Fe/H] > −0.80 dex, while exclusively observing the higher energy populations of GES.As shown in Koppelman et al. (2020), the high-energy debris of a merger event typically consists of the first stars to have been lost by a progenitor, meaning those in its peripheries.As commonly observed in dwarf spheroidal galaxies orbiting the Milky Way (see, e.g., Tolstoy et al. 2023), it is reasonable to expect a negative metallicity gradient that would explain why our targets are slightly more metal-poor than lower-energy GES samples.Stars assigned to Sequoia are, on average, slightly more metal-poor with respect to GES and display a broad MDF, with a median metallicity of [Fe/H] = −1.66dex and σ [Fe/H] = 0.37, which is consistent with the result by D23.The two most metal-poor stars in Sequoia have a metallicity that overlaps with that of ED-2.However, these 2 stars occupy a place in the IoM space that is significantly different with respect to the ED-2 clump.Moreover, they do not exhibit the distinctive high [Ti/Fe] ratios observed in the stars of ED-2 at a comparable metallicity.As a result, we exclude a possible association of these stars with the ED-2 substructure based on their chemical features.To obtain a clearer view, we explored the space defined by the orbital eccentricity, defined as ecc = (apo − peri)/(apo + peri), where apo and peri are the orbital apocenter and pericenter, respectively; this dynamical parameter can help to efficiently differentiate the various substructures.Each substructure comprises stars that predominantly follow orbits characterized by very similar eccentricity (see Fig. 10).The exceptions are Clump-6 and Clump-7, whose stars are spread across the range of ecc 0.1−0.8.
The broad spread in the eccentricity values for the stars in these substructures, coupled with the shape of their MDFs, compels us to attribute the existence of these cluster to a combination of our selection criterion and the choice of hyperparameters in the clustering algorithm.As a result, we do not consider Clump-6 and Clump-7 to be genuine substructures or independent accretion events.Thus, we excluded stars associated to these two clusters from the rest of the analysis.
Based on the MDF and the eccentricity parameter, three substructures emerge as notably more coherent.GES distinguishes itself as the most metal-rich component, characterized by a peak in the MDF that closely mirrors the primary peak observed in the overall sample.Furthermore, GES comprises stars that exhibit the highest eccentricities in their orbital motion.Then, Antaeus is identified as another well-defined substructure, with its constituent stars primarily having orbits characterized by a relatively high eccentricity ( 0.7).Finally, ED-2 exhibits distinctive features in its MDF and ecc, both being characterized by a remarkably narrow distribution.It is interesting to note that the 2 stars that were incorporated into the ED-2 clump by the clustering algorithm have [Fe/H] −1.90 dex and also showcase a behavior that is more akin to that of Sequoia in terms of eccentricty (ecc > 0.6) and [TiI/Fe].These values support our choice to incorporate them as part of Sequoia.

Odd-Z elements
Odd-Z elements, such as Na and Al, are believed to be mainly synthesized in massive stars exploding as core-collapse supernovae (CC-SNe) at low metallicity and in asymptotic giant branch (AGB) stars at intermediate metallicity (Nomoto et al. 2013).Since they are formed starting from a neutron excess produced in the CNO cycle, their yield depends on the metallicity of the progenitor (Kobayashi et al. 2020).
Recent studies have demonstrated that stars that have formed in situ are enhanced in odd-Z elements abundances with respect to accreted populations and stars in dwarf satellite galaxies of the MW (Tolstoy et al. 2009;Feuillet et al. 2021;Hasselquist et al. 2021;Belokurov & Kravtsov 2022).
We show in Fig. 12 that the stars from our sample display sub-solar values of [Na/Fe], in agreement with the findings by Nissen & Schuster (2010).Moreover, it is evident that Antaeus displays a depletion in Na relative to GES.This distinction is discernible through the positioning of its stars on a distinct sequence, consistently registering values ∼0.2 dex lower.The sequence defined by Antaeus closely resembles the pattern found by Matsuno et al. (2022) for the stars belonging to Sequoia, whereas our stars associated with Sequoia exhibit a behavior that is more akin to that of GES.However, Matsuno et al. (2022) employed a relatively wide retrograde region within the IoM space to identify members of Sequoia (with L z values ranging from −1600 to −3100 kpc km s −1 ), making very likely that some Antaeus stars inadvertently contaminate their Sequoia sample, thus potentially contributing to the observed similarities.Figure 12

α-elements
With the notable exception of oxygen, the α-elements are mainly synthesized in the core of massive stars and diffused into the interstellar medium through CC-SNe, while a small fraction is also produced in SN Ia (Kobayashi et al. 2020).Thus, the chemical patterns of these elements can be used as a proxy of the different timescales of star formation.A typical indicator of the efficiency of a galaxy's star formation is the metallicity of the knee in the [α/Fe] vs. [Fe/H] plane, which denotes the beginning A37, page 15 of 21  (2003,2006), Barklem et al. (2005), Bensby et al. (2005Bensby et al. ( , 2014)), Roederer et al. (2014) and Reggiani et al. (2017).
of a sizable contribution to the chemical enrichment by SNe Ia (Tinsley 1979;Matteucci & Greggio 1986).Furthermore, Mg is synthesized through the hydrostatic burning of He, C, and Ne in massive stars (30-35 M ) with a negligible contribution by SN Ia (Woosley & Weaver 1995;Arnett 1996).The explosive α-elements (Ca and Ti) are mainly produced in less massive stars (15-25 M ), with a considerable fraction formed also in SN Ia (Kobayashi et al. 2020).
As shown in Fig. 13, the interpretation of our findings is rather challenging.The spread of the measured [α/Fe] is significantly larger than our uncertainties.The fact that our sample is systematically more metal-poor than typical disc stars makes the comparison with the APOGEE control sample not trivial (see Mg and Ca in Fig. 13).Indeed, the [Mg/Fe] of our sample seems systematically lower at [Fe/H] > −1.5 dex, when the patterns expected for a dwarf galaxy start to part ways from that expected for the MW (Venn et al. 2004;Tolstoy et al. 2009;Nissen & Schuster 2010;Hayes et al. 2018;Horta et al. 2023).In the most metal-poor regime, the two populations mix, but this is not surprising given that SN Ia still have to start exploding.Conversely, for elements like Ti, where the reliability of APOGEE derived abundances is uncertain, we have opted to juxtapose our findings with established high-resolution spectroscopy literature data.As demonstrated in Fig. 13, our results exhibit a strong concordance with those obtained MW halo.With a knowledge of this general behavior, we can now focus in detail on the α-elements trends of the six substructures and discuss the possible origin of the distinct chemical patterns, described below.
Gaia-Enceladus.The Mg abundance ratios of the GES stars overlap the metal-poor tail of the distribution of the low-α MW population, as already observed in previous works (Helmi et al. 2018;Mackereth et al. 2019).The reason for the difference in abundance between GES and high-α stars is commonly believed to be the result of chemical enrichment by SN Ia.In fact, in situ stars have not undergone significant enrichment by SN Ia at [Fe/H] −1.0 dex, whereas a dwarf galaxy like GES has already experienced such enrichment due to a less efficient star formation history.
Sequoia.The stars associated to Sequoia exhibit flat [α/Fe] patterns across the entire metallicity range.In particular, Sequoia occupies the same region of GES at [Fe/H] > −1.5 dex, especially in the chemical planes of TiI.At lower metallicity, however, the behavior exhibited by the Sequoia substructure more closely resembles that of the smaller retrograde substructures.At the same time, Sequoia seems to display [Mg/Fe] and [Ca/Fe] ratios that are 0.1-0.2dex lower than GES, in agreement with the results by Matsuno et al. (2022).Based on these findings, we favor the interpretation of Sequoia as an independent merger A37, page 16 of 21 event with respect to GES.Lower [α/Fe] values cannot be explained by the scenario that considers Sequoia as the outskirt of GES that were first lost during the tidal disruption of the galaxy, as this would mean that the outer regions were polluted by SN Ia ejecta before the inner ones, which is a feature yet not observed in existing dwarf MW satellites.
Antaeus.The abundances of the α-elements suggest that Antaeus follows a distinctive chemical sequence in the all the [α/Fe] planes at lower values at fixed metallicity with respect to the dominant components in the RH (GES and Sequoia).At [Fe/H] −1.5 dex, the [α/Fe] of Antaeus is consistently 0.1-0.2dex lower than those observed in GES.We interpret this feature as evidence for the existence of Antaeus as an independent merging event, characterized by a progenitor with lower SFR compared to GES and Sequoia.
ED-3.This substructure is located in a region of the [α/Fe] chemical plane that is consistent with that of Antaeus.Constituted of only five stars, ED-3 overlaps with the Antaeus clump across its entire metallicity range.The two systems display a high degree of proximity in terms of their E and L z , resulting in significant overlap within this parameter space.Based on this new chemo-dynamical evidence, we conclude that Antaeus/L-RL64 and ED-3 are actually part of the same relic.This evidence is confirmed in all of the chemical planes, as we will show in the next paragraphs.
ED-2.The four stars associated with ED-2 occupy the very metal-poor end of the chemical planes, being all located at −2.6 < [Fe/H] < −2.4 dex.This system exhibits a markedly narrow metallicity distribution and its chemical composition is generally consistent with that of the other RH substructures, with the exception of Ti, where it seems to be enhanced.Its position in the [α/Fe] plane seems to recall the one of I'itoi (Naidu et al. 2022).Given the small sample, any conclusions inferred about the chemical pattern of this system would be premature and require further investigation on a much extended sample of stars.

Iron-peak elements
The production of iron-peak elements is sourced from a variety of nucleosynthesis formation pathways, with SN Ia, CC-SNe, and hypernovae all making a sizeable contribution to the overall production of these chemical elements (Romano et al. 2010;Leung & Nomoto 2018;Kobayashi et al. 2020;Lach et al. 2020).
Sc is mainly produced in CC-SNe in the burning of C and Ne (Woosley et al. 2002), with a negligible contribution by SN Ia (Kobayashi et al. 2020), while Mn is almost entirely produced by SN Ia, with a relative contribution by CC-SNe and hypernovae even lower than Fe (Kobayashi & Nomoto 2009;Kobayashi et al. 2020).SN Ia are also the predominant channel of formation of Ni, although it can be produced in CC-SNe, as for Fe.Overall, Zn is synthesized in the highest energy hypernovae.
As shown in Fig. 14, all the RH substructures exhibit abundances that are remarkably similar to those observed in accreted populations and low-mass satellite galaxies of the MW (Venn et al. 2004;Nissen & Schuster 1997;Tolstoy et al. 2009).The interpretation of these results is complex, as the small differences between the various substructures can be explained by both the diverse production sites of these elements and the metallicity dependencies of their yields at low metallicities (Weinberg et al. 2022).The chemical patterns of the iron-peak elements for the different substructures seem to reveal that GES and Sequoia display a similar distribution.Moreover, Antaeus and ED-3 exhibit, once again, an overlap in their respective loci within these chemical plane, providing additional support for the hypothesis that they are associated with the same accretion event.

Neutron-capture elements
Of all neutron-capture elements, we focus only on yttrium (Y) in this work.This element is predominantly formed by the slow neutron capture process (s-process, i.e., neutron capture that is slower than the β-decay).
Y is predominantly formed by s-processes in the envelope of low-mass (1-4 M ) AGB stars and, to a lesser extent, in massive stars (Busso et al. 1999;Pignatari et al. 2010;Cristallo et al. 2015).At low metallicities, magneto-rotational SNe and mergers of neutron stars contribute significantly to the production of these neutron-capture elements via r-processes, thus explaining the observed scatter (Kasen et al. 2017;Molero et al. 2023).
Despite the huge scatter between abundances of stars born in the same progenitor, the results in Fig. 15 show that the average abundances for the Sequoia substructure is comparable to the one of GES, as previously hinted also by Aguado et al. (2021).Within the common metallicity range, Antaeus and ED-3 exhibit a pronounced degree of overlapping, consistently displaying average abundances that are lower (∼0.1-0.2 dex) than those observed in GES.

Discussion and conclusions
In this first paper of a series, we introduce the WRS project, describing the methodology employed and highlighting the scientific goals it aims to pursue.To demonstrate the potential of this project, we have narrowed our focus to the subset composed of 90 stars associated by the DBSCAN clustering algorithm with various RH substructures within the IoM space (see Dodd et al. 2023) and we have analyzed their chemical properties.The interpretation of the chemical patterns of unassociated stars will be left to future contributions to the series.We present a concise summary of the principal findings of this paper below.1.A selection criterion based on tangential velocity only (thus not requiring the V los information) is able to select stars that move on retrograde motions with a 91.4% efficiency (by requiring V T > 400 km s −1 ).All of the stars included in our sample consistently exhibit chemical pattern across all chemical element planes that aligns with the characteristic abundance patterns displayed by accreted systems and dwarf galaxies.2. In spite of the selection criteria leaving out the bulk of GES stars at higher L z and lower energy from our sample, GES stands out as the largest and most metal-rich component within the RH.This supports the hypothesis that it represents the remnant of the most massive galaxy that was accreted by the MW. 3. The chemical patterns of GES and Sequoia display similar trends at [Fe/H] > −1.5 dex.However, significant distinctions become evident at lower metallicity, specifically characterized by a discernible difference in the α-elements.In particular (as shown in the top panel of Fig. 16), Sequoia exhibit a pronounced excess of stars with lower [Mg/Fe] and [Ca/Fe] ratios with respect to GES.Consequently, this discrepancy leads us to interpret Sequoia as an independent accretion event distinct from GES.   Gratton et al. (2003), Reddy et al. (2003Reddy et al. ( , 2006)), Bensby et al. (2005Bensby et al. ( , 2014)), Roederer et al. (2014) and Reggiani et al. (2017).Edvardsson et al. (1993), Fulbright (2000), Stephens & Boesgaard (2002), Reddy et al. (2003Reddy et al. ( , 2006)), Bensby et al. (2005) and Reggiani et al. (2017).
4. For the first time, we have achieved a comprehensive chemical characterization of Antaeus, encompassing detailed chemical abundances for numerous elements.Our analysis reveals distinct behavior in the chemical space compared to GES within the same metallicity range (see bottom panel of Fig. 16), reinforcing the assertion that it constitutes the remnant of an accretion event independent of GES. 5. Antaeus and ED-3 exhibit analogous positions in the IoM spaces, albeit with some distinctions in L ⊥ .Additionally, they share similarities in terms of their metallicity distribution function and orbital eccentricity.Furthermore, the chemical patterns observed in Antaeus and ED-3 are consistent across all chemical elements.In the end, we integrated backwards in time (for 2.5 Gyr) the orbits of the stars associated to these two substructures in a McMillan (2017) potential for the MW.As depicted in Fig. 17, stars from the two substructures appear to move on very similar orbits.Collectively, these lines of evidence compel us to interpret them as components of the same merging event.Henceforth, we refer to this merging event as Antaeus.6. ED-2 is the most metal-poor component, with a remarkably narrow MDF peaked at [Fe/H] = −2.43dex, in agreement with the findings of Balbinot et al. (2023).The limited number of stars within the ED-2 substructure, coupled with the specific combination of their stellar parameters (T eff > 6000 K and [Fe/H] < −2.4 dex) leading to the manifestation of weak spectral lines, restricts our capacity to attain detailed chemical abundances across all elements.As a consequence, we are unable to definitively assert whether ED-2 constitutes A37, page 18 of 21

Fig. 1 .
Fig. 1.Toomre diagram of the 2.4 million stars within 1.0 kpc from the Sun, having V los values measured in Gaia EDR3 and passing all the quality criteria we adopted for the WRS sample.Stars are color-coded according to their tangential velocity and stars with V T 400 km s −1 are highlighted with a black empty square.

Fig. 3 .
Fig. 3. Line-of-sight velocity distribution for the entire stellar sample, shown on the left.Central panel shows comparison between the uncertainty on the line-of-sight velocity provided by Gaia (blue filled circles) and the one measured with high-resolution spectroscopy (red filled circles) as a function of observed G magnitude.The blue arrows point out the outliers in the distribution, specifically representing target stars with errors in the Gaia velocity measurement larger than 10 km s −1 .Right panel shows the difference between the two measurements of the line-of-sight velocity.The red dashed lines mark the ±3σ threshold.

Fig. 4 .
Fig. 4. Distribution of the observed stars in IoM and velocity space (top and bottom panels, respectively).We highlight stars associated with various substructures with different colors (green for GES, brown for Sequoia, blue for Antaeus, magenta for ED-2, orange for ED-3, cyan for Clump-6 and red for Clump-7).

Fig. 5 .
Fig. 5. Gaia EDR3 distance-corrected color-magnitude diagram of the 175 analyzed stars.Targets are color-coded based on the metallicity.

Fig. 6 .
Fig. 6.Comparison of Mg abundance ratios of the stars from our selection observed with UVES (blue filled circles) and the NS10 reference sub-sample (red filled squares and green filled triangles for low-α and high-α sequences, respectively).The errorbar in the top-right corner indicates the typical uncertainties.

Fig. 7 .Fig. 8 .
Fig. 7. [Mg/Fe] abundance ratios as a function of [Fe/H] for our sample (color-coded by the angular momentum along the z-direction, L z ) and for the APOGEE control sample (greyscale 2D histogram) shown on the left.The errorbar in the top-right corner displays the typical uncertainties.[Mg/Mn] abundance ratios as a function of [Al/Fe] for the same samples, with the same arrangement shown on the right.Red arrows indicate upper limits.

Fig. 9 .
Fig. 9. Metallicity distribution functions for the RH substructures labeled at the top of each panel compared to the MDF of the entire sample of stars.The median metallicity value and the standard deviation of the distribution are also reported.The color coding is the same as in Fig. 4.

Fig. 10 .
Fig. 10.Eccentricity distribution functions for the RH substructures labeled at the top of each panel compared to the distribution of the entire sample of stars.The median eccentricity value and the standard deviation of the distribution are also reported.The color coding is the same as in Fig. 4.

Fig. 11 .
Fig. 11.Distribution of the RH substructures in the eccentricity vs. [Fe/H] space (left panel) and in the eccentricity vs. Z max (right panel).The color coding is the same as in Fig. 4.
also shows that the [Al/Fe] trend of our targets is almost flat at [Al/Fe] ∼ −0.3 dex for [Fe/H] < −1.0 dex.The [Al/Fe] abundances observed in the majority of the target stars show a distinctively lower level compared to the typical values observed in the MW disc, thereby indicating a consistency with the position occupied by accreted systems.Interestingly, we find a spread of about 0.5 dex in the values of [Na/Fe] and [Al/Fe] for the associated stars, which might hint at a different fraction of massive stars exploding as hypernovae in the progenitors(Smiljanic et al. 2016).

Table 1 .
Main target information: ID, coordinates, apparent magnitudes and parallax from Gaia EDR3, color excess from the EXPLORE tool, and the instrument used to collect the spectrum.

Table 2 .
Magnitude, exposure times, airmass, and S/Ns for some of the observed target spectra with the spectrographs UVES and PEPSI.
Comparison between the spectra of two stars observed with PEPSI at LBT, shown on the left.Right column shows the same comparison, but for two stars observed with UVES at VLT.These stars have very similar T eff and log g but they differ in the iron content.Black arrows indicate the position of some atomic lines of interest.

Table 4 .
Stellar parameters for the selected targets: ID from Gaia EDR3, derived effective temperature, surface gravity, microturbulent velocity, and iron abundance ratio.The complete table is available at the CDS. Notes.

Table 5 .
Differences in chemical abundances obtained from the UVES and PEPSI samples, including the calculated offsets, corresponding standard deviations, and the number of stars used for each element.

Table 6 .
Chemical abundances for the odd-Z, iron-peak and neutron capture elements for the stars associated to known substructures.Notes.The values of [Na/Fe] are corrected for NLTE effects.The entire table is available at the CDS.

Table 7 .
Chemical abundances for the α-elements for the stars associated to known substructures.
Notes.In the last column we report our association with the former progenitor and the result provided by D23 (inside the brackets).The entire table is available at the CDS.