Milky Way's eccentric constituents with $Gaia$, APOGEE $\&$ GALAH

We report the results of an unsupervised decomposition of the local stellar halo in the chemo-dynamical space spanned by the abundance measurements from APOGEE DR17 and GALAH DR3. In our Gaussian Mixture Model, only four independent components dominate the halo in the Solar neighborhood, three previously known $Aurora$, $Splash$ and Gaia-Sausage/Enceladus (GS/E) and one new, $Eos$. Only one of these four is of accreted origin, namely the GS/E, thus supporting the earlier claims that the GS/E is the main progenitor of the Galactic stellar halo. We show that $Aurora$ is entirely consistent with the chemical properties of the so-called Heracles merger. In our analysis in which no predefined chemical selection cuts are applied, $Aurora$ spans a wide range of [Al/Fe] with a metallicity correlation indicative of a fast chemical enrichment in a massive galaxy, the young Milky Way. The new halo component dubbed $Eos$ is classified as $in situ$ given its high mean [Al/Fe]. $Eos$ shows strong evolution as a function of [Fe/H], where it changes from being the closest to GS/E at its lowest [Fe/H] to being indistinguishable from the Galactic low-$\alpha$ population at its highest [Fe/H]. We surmise that at least some of the outer thin disk of the Galaxy started its evolution in the gas polluted by the GS/E, and $Eos$ is an evidence of this process.


INTRODUCTION
One curious thing to ponder about our Galaxy's continuous evolution is that the guises the Milky Way takes as it transforms all remain in place -albeit somewhat hidden -for us to make sense of. Today, with new data that let us see the Galaxy in phase space (thanks to the ESA's space observatory Gaia (Gaia Collaboration et al. 2021)) and in chemical space (thanks to massive high-resolution surveys like Apache Point Observatory Galactic Evolution Experiment (APOGEE; Abdurro'uf et al. 2022) and Galactic Archaeology with HERMES (GALAH; Buder et al. 2021a)), it is possible to reveal the Milky Way's forgotten states and discarded guises. The Galactic halo specifically has remained the centre of intense attention thanks to its capacity to keep orbital information neatly preserved for every star deposited -this at a clear advantage compared to the studies of the disk where the past state is less straightforward to decipher due to the effects of radial mixing (see e.g., Sellwood & Binney 2002;Schönrich & Binney 2009;Loebman et al. 2011;Minchev et al. 2013;Frankel et al. 2020).
The halo is vast and even though its deep trawl only just started, the haul offered up by Gaia has been copious. For example, a large number of low-mass accreted systems were detected in the Solar neighbourhood (see e.g., Myeong et al. 2018a,b;Koppelman et al. 2019;O'Hare et al. 2020;Borsato et al. 2020;Limberg et al. 2021;Sofie Lövdal et al. 2022) and beyond (see e.g., Malhan et al. 2018;Ibata et al. 2019;Palau & Miralda-Escudé 2019;Naidu et al. 2020;Yuan et al. 2020;Ibata et al. 2021). On the other hand, clear evidence for rem-nants of more massive galaxies consumed by the Milky Way appears rather scarce indicating a modest and restrained accretion course our Galaxy has followed.
So far, only one ancient and massive intruder has been identified unambiguously, with its body now nothing but an enormous cloud of tidal wreckage scattered throughout the inner Milky Way. The existence of a vast debris structure left behind by this merger event, known today as the Gaia Sausage/Enceladus (GS/E), was suggested before Gaia (see Evans 2020, for the historical development). In particular, Deason et al. (2013) argued that the rapid transition in the Galactic stellar halo structural properties at break radius of 20-30 kpc is likely associated with the apocentric pile-up of a relatively early (8-10 Gyr ago) and single massive accretion event. It is hypothesised that the progenitor dwarf galaxy was massive enough for its orbit to shrink and radialise quickly as a result of a complex interplay between dynamical deceleration, host recoiling and self-friction (Vasiliev et al. 2022). Sinking deep in the heart of the Milky Way, the dwarf sprayed the bulk of its stars in a region enclosed by its last apo-centre, i.e. some ∼ 30 kpc (see Deason et al. 2018). As a result, the region of the Galactic halo within this so-called break radius (see Deason et al. 2011) is inundated with the GS/E stars, consistently showing up as the most striking sub-structure even in relatively small volumes around the Sun (see examples of pre-Gaia hints in e.g., Chiba & Yoshii 1998;Brook et al. 2003;Meza et al. 2005;Nissen & Schuster 2010;Hawkins et al. 2015). Subsequently, the Gaia data made crystal clear the unusually strong radial anisotropy of the relatively metalrich GS/E debris and helped to reveal its dominance in the Solar neighbourhood (see Belokurov et al. 2018;Myeong et al. 2018c;Haywood et al. 2018b). The genesis of the GS/E debris cloud was also unambiguously confirmed by its unique chemical fingerprint (see e.g., Helmi et al. 2018;Mackereth et al. 2019) and the large group of globular clusters (GCs) associated with it (see e.g., Myeong et al. 2018dMyeong et al. , 2019Massari et al. 2019). Close to the Sun, the GS/E structure has been thoroughly scrutinised both kinematically and chemically (see Necib et al. 2019;Evans et al. 2019;Sahlholdt et al. 2019;Das et al. 2020;Monty et al. 2020;Kordopatis et al. 2020;Feuillet et al. 2020;Molaro et al. 2020;Carollo & Chiba 2021;Aguado et al. 2021a;Matsuno et al. 2021;Bonifacio et al. 2021;Feuillet et al. 2021;Buder et al. 2021b). Outside of the Solar neighborhood, fewer studies exist, nonetheless the global structure of the GS/E cloud is starting to come into focus (see e.g., Simion et al. 2019;Lancaster et al. 2019;Naidu et al. 2020;Bird et al. 2021;Balbinot & Helmi 2021;Iorio & Belokurov 2021;Wu et al. 2022).
Two lines of enquiry help to cement the unique place the GS/E merger holds in the history of our Galaxy: its mass and the time of its arrival. Using established mass-metallicity relations (see Kirby et al. 2013), the total stellar mass of the progenitor dwarf before disruption can be gleaned from the typical metallicity of its stellar population Feuillet et al. 2019;Naidu et al. 2021). An independent mass constraint can be obtained by scaling up the total number of globular clusters associated with the merger (Myeong et al. 2018d(Myeong et al. , 2019Massari et al. 2019;Forbes 2020;Callingham et al. 2022). A third method of constraining the total stellar mass is via the analysis of the GS/E chemical pattern (e.g., Fernández-Alvar et al. 2018;Helmi et al. 2018;Vincenzo et al. 2019;Mackereth et al. 2019). All of the above estimates agree that the stellar mass of the GS/E progenitor had to be of order of few ×10 8 M . This can be contrasted with the total stellar halo mass measurements of ∼ 10 9 M (see Deason et al. 2019;Mackereth & Bovy 2020). As the comparison shows, the GS/E merger must have contributed a significant fraction of the total stellar halo, perhaps as much as a half (for the inner ∼ 50 kpc region) in agreement with local ) and global  estimates. In other words, a simple budget of the Galactic stellar halo leaves little room for other significant accretion events.
The dominance of the GS/E merger also makes sense when the timeline of the Galactic accretion is considered. Constraining the time of the accretion is far from straightforward. Not many asteroseismological ages exist for relatively metal-poor old stars (see Montalbán et al. 2021;Grunblatt et al. 2021;Alencastro Puls et al. 2022). These direct, albeit scarce, measurements place the epoch of the GS/E merger at 8-10 Gyr ago, in agreement with White Dwarf cooling ages (Fantin et al. 2021). The age estimates for a bulk of GS/E stars based on main-sequence-turn-off stars ) and red giants (Das et al. 2020) also show good agreements with this time frame. Alternatively, the time of accretion can be obtained via the inference of the starformation history (SFH) of the progenitor galaxy (e.g., Helmi et al. 2018;Vincenzo et al. 2019;Mackereth et al. 2019;Sanders et al. 2021). These modelling attempts agree that the dwarf's SFH was shut down between 8 and 10 Gyr ago. Finally both the progenitor mass and the accretion time can be gauged by comparing the observed redshift z = 0 properties of the tidal debris cloud with its numerical counterparts in galaxy formation simulations (see Belokurov et al. 2018;Fattahi et al. 2019;Bignone et al. 2019;Elias et al. 2020;Dillamore et al. 2022). Across different simulation suites, these stud-ies find that the accretion events that lead to creation of a highly radially anisotropic metal-rich stellar halo structure near the Sun typically take place 7-10 Gyr ago and that the progenitor galaxies bring between 10 8 and 10 9 M of stars.
Does this imply that the GS/E merger is the only significant merger our Galaxy has experienced? Not necessarily. Remnants of even older interactions could potentially be hidden in the inner Galaxy, i.e. the area known as the bulge and, indeed, several such claims have been made recently (e.g., Kruijssen et al. 2019;Forbes 2020;Horta et al. 2021). Additionally, more recent or ongoing accretion events may be lurking in the Milky Way's outskirts. There are at least two such relatively recent and relatively massive accretion events: that of the Sagittarius dwarf galaxy (see Ibata et al. 1994;Majewski et al. 2003) and of the Magellanic System (see Wannier & Wrixon 1972;Mathewson et al. 1974). The GS/E, the Magellanic Clouds and the Sgr dwarf are all estimated to have been hosted in dark matter halos with masses not far off 10 11 M (see e.g., Fattahi et al. 2019;Erkal et al. 2019;Gibbons et al. 2017). In addition, tidal debris from several systems with at least an order of magnitude lower mass have also been found, e.g., Sequoia (see Myeong et al. 2018bMyeong et al. , 2019Matsuno et al. 2019).
The analysis reported here is inspired by the recent discoveries of a large number of diverse halo substructures and the need for a systematic tally of past mergers with a focus on the most important events. To our knowledge, not many attempts have been made to build a well-defined census of Galactic accretion remnants. Taking advantage of a rich but relatively small Galactic Globular Cluster dataset, several groups sought to apportion the known GCs to a small number of significant mergers Kruijssen et al. 2019;Forbes 2020;Callingham et al. 2022). However, the choices as to the decision boundaries and the number of independent components were made by hand, and therefore, unfortunately, were neither objective or reproducible (in e.g., a numerical simulation). In a similar vein, Naidu et al. (2020) suggested a break-down of the entire stellar halo accessible to the H3 survey into individual accretion events but their substructure membership, similarly to the works above, was designed by hand. Most similar to our study are the works of Das et al. (2020) and Buder et al. (2021b) who both aimed to extract the footprint of the local GS/E tidal debris in the space spanned by chemical abundances. They proceeded by modelling high-quality APOGEE and GALAH samples with a mixture of Gaussians and identifying the component representing the GS/E against the background of other Galactic populations.
The aim of our study is to build an unbiased and selfconsistent inventory of the stellar halo in the vicinity of the Sun. Our main objective is to identify and characterise all of the significant independent halo components without postulating their number and their behaviour a priori. For the task in hand we will also use an approach based on a mixture of Gaussians. The stars belonging to the in situ populations such as the disk are expected to dominate the halo denizens by almost two orders of magnitude. To avoid wasting model components on known Galactic populations at the expense of the objects of our experiment, we limit the scope of our study solely to stars on orbits with high eccentricity. Thus, the stellar sample considered consists only of stars unambiguously attributed to the Galactic halo. While demanding high eccentricity ensures that the halo sample is pure, it necessarily leads to incompleteness because tidal debris with significant angular momentum will be excised. In terms of the accreted sub-structures, there are plenty of local examples with both prograde (e.g., Helmi and the associated S2 stream, Helmi et al. 1999;Myeong et al. 2018a;Aguado et al. 2021b) or retrograde (e.g., Sequoia, Myeong et al. 2018bMyeong et al. , 2019Matsuno et al. 2019) orbital motion.
Note that thanks to Gaia we are now certain that stars on high-eccentricity, or halo-like, orbits do not necessarily have to be accreted. As first demonstrated by Bonaca et al. (2017) roughly half of the local stellar halo is of high metallicity, i.e. [Fe/H]> −1 dex and is likely of in situ origin. These in situ halo stars had of course been seen before, e.g., as the high-α halo sequence by Nissen & Schuster (2010) and Hawkins et al. (2015) as later confirmed by e.g., Haywood et al. (2018b). The origin of the in situ halo stars with [Fe/H]> −1 dex has been linked to the GS/E accretion event (Gallart et al. 2019;Di Matteo et al. 2019;Belokurov et al. 2020). These stars were born in the ancient high-α disk of the Milky Way before the merger and subsequently "splashed" out of their circular orbits during the interaction with the GS/E progenitor, a scenario ubiquitous in numerical simulations of galaxy formation (see Belokurov et al. 2020;Grand et al. 2020;Dillamore et al. 2022). Not all of the in situ halo is splashed high-α disk. As Belokurov & Kravtsov (2022) demonstrate using APOGEE DR17 data, the Milky Way stars with [Fe/H]< −1.5 dex were likely formed in a hot and messy pre-disk state that later phase-mixed to form a centrally concentrated spheroidal in situ halo component dubbed Aurora. Conroy et al. (2022) also see evidence for an ancient Galactic pre-disk population in the low-metallicity sample of H3 stars. This paper is structured as follows. Section 2 describes the two sets of selection cuts applied to the APOGEE and the GALAH data to produce the samples of high-quality abundance measurements analysed. As the APOGEE DR17 reports abundances with noticeably smaller uncertainties compared to GALAH we consider the APOGEE sample as our primary dataset. Section 3 briefly gives the details of the Gaussian Mixture Model methodology used here. We analyse the components extracted from APOGEE and GALAH in Section 4. Our results are placed in context and summarised in Section 5. The sample is cleaned by rejecting any stars with STAR BAD flag in ASPCAPFLAG, and PROGRAMNAME = magclouds. Only the main red star samples with EXTRATARG = 0 are used. To obtain a red giant star sample, logg < 3.0 is adopted. It is further filtered with x fe flag = 0, x fe err < 0.1 dex for almost all elements x used in our study. The exception is [Ce/Fe], for which ce fe err < 0.15 dex is used instead of 0.1. To focus on the halo or halo-like Galactic components such as the GS/E Haywood et al. 2018b;Myeong et al. 2018c,d) and the Splash , an eccentricity > 0.85 cut is adopted to obtain stars on nearly radial orbits (see e.g., Belokurov et al. 2018;Myeong et al. 2018cMyeong et al. , 2019Mackereth et al. 2019;Buder et al. 2021b). We also require the orbital apocenter > 5 kpc to minimize the contamination from the bulge populations (see e.g., Horta et al. 2021;Buder et al. 2021b). In addition, the distance error < 1.5 kpc, and orbital energy E < 0 km 2 s −2 cuts are applied.
For  [Mg/Mn], and the orbital energy (E) information, forming a 6-dimensional chemodynamical space to study with Gaussian Mixture Modelling (GMM). There are a total of ≈1700 APOGEE stars that pass our quality cuts 1 . The choice of our 6-dimensional space is motivated by the desire to include dynamical properties (E) as well as the "most reliable" (or "reliable" if necessary) elemental production channels for the categories of odd-Z (Al), α-(α, Mg), iron-peak (Mn), and s-process (Ce) elements, following APOGEE's recommendation for giant stars.
As a cross-check, we also use the main and the VAC data of GALAH DR3 (Buder et al. 2021a) which contain 678,423 spectra of 588,571 mostly nearby stars. The catalogues provide stellar parameters, radial velocity, and abundances for up to 30 chemical species. The VAC also provides various additional information including the distance estimates (see also Sharma et al. 2018) and Galactic orbital dynamics calculated using galpy (Bovy 2015) with the McMillan (2017) potential.
We adopt the recommendations from the GALAH Collaboration for the choice of data columns and quality cuts 2 . Our sample is first filtered with snr c3 iraf > 30, flag sp = 0. We adopt similar cuts from our APOGEE sample, with the exception of using e x fe < 0.2 dex for elements x used for our study.
For This includes similar tracers as the APOGEE sample, plus the r-process element (Eu) available in GALAH. Therefore for GALAH, the GMM is applied to a 12dimensional chemo-dynamical space sampled with a total of ≈1100 stars that passed the above quality cuts 3 .

METHOD
Gaussian Mixture Modelling is an unsupervised probabilistic modelling technique. It assumes that the overall distribution of some (unlabelled) data is the sum of a finite number of Gaussian distributions with unknown parameters. The optimal (maximum likelihood) model parameters can be searched for using Expectation-Maximization (EM) iterations. For our modelling, we use the Extreme Deconvolution (Bovy et al. 2011) (XD) algorithm. One attractive feature of XD is the ability to include uncertainties associated with the observed data. The maximum number of EM iterations for each XD run is 10 9 by default We run XD multiple times with different numbers of Gaussian components, from N = 1 to 10. For each model with N Gaussians, we repeat the fit 100 times with different initial parameters randomly drawn from the parameter space. This is to insure that the algo-rithm obtains the true maximum likelihood, not a local maximum. The whole session of 100 random runs is also repeated multiple times to verify convergence to the best-fitting case. We check the results of the XD fitting with a different GMM package from Scikit-learn (Pedregosa et al. 2011) to demonstrate consistency.
The APOGEE sample provides abundances with noticeably smaller uncertainties as compared to GALAH. Therefore, the APOGEE-Gaia data serves as the primary sample for our study. The difference in the quality of the data is apparent in the overall level of scatter seen in each dataset 4 .
To avoid potential over-fitting, we compute the Bayesian Information Criterion (Schwarz 1978) (BIC) and the Akaike Information Criterion (Akaike 1974) (AIC) for each XD run. We use the BIC score to assess the necessary number of components for each model. The lowest BIC score (highest maximum log-likelihood among all XD runs) indicates the GMM with 7 components as the optimal choice. Although the BIC scores for 5 to 8 Gaussian components cases appear relatively comparable, the 7 Gaussians case always scores the lowest BIC consistently from our repeated trials of XD fitting sessions. Since we repeat the XD fitting 100 times during each trial (with randomized initial parameters) for each model with N Gaussians, we check the overall distribution of the score in each case. When compared across all models (from N = 1 to 10), the change in the BIC distribution shows that the median and the highest score are also the lowest at N = 7 case in addition to the fact that it achieves the overall lowest (preferred) BIC value. Figure 1 shows the scores from the BIC and AIC tests. For each model with N Gaussians, the solid line gives the lowest score out of 100 XD runs. The dashed line is the median, and the dotted line is the highest score for each N Gaussian case. Note that, although the BIC criterion provides a clear minimum, the AIC does not. The differences between AIC and BIC are discussed in standard texts on Bayesian model fitting using information criteria (e.g., Gelman et al. 2015). Here, we merely note that model selection is a difficult task, without a clear-cut answer as to which information criteria works best. For each number of components, the GMM fit is performed a large number of times starting from random initial conditions to ensure convergence and to investigate scatter in obtained solutions. For each number of components probed, solid line gives the lowest (best) AIC/BIC score, dashed line corresponds to the median, and dotted to the highest score. The BIC test indicates the GMM with 7 Gaussian components as the optimal case. The zoom-in panel shows the BIC result around the best score.

GMM view of the local halo
GMM is an unsupervised probabilistic modelling technique. Naturally enough, it recovers the chemical signatures of the three known eccentric components of the Galaxy, namely the GS/E, the Splash and Aurora. This is in itself valuable as it confirms earlier results derived by extracting samples of these populations by handmade cuts. Equally, this confirmation adds extra confidence to the GMM's results when the method breaks completely new ground.
For each GMM component in the APOGEE-Gaia sample, the parameters describing the weight (fractional contribution), the mean and the standard deviation σ of the model are presented in Table 1. The value for the combined GS/E population is also marked as GS/E sum . The sum of the two low-amplitude physically implausible components, associated with the background, (back 1 , back 2 ) are also marked as back sum for a reference. Figure 2 shows several 2-dimensional projections of the original 6-dimensional chemo-dynamical    APOGEE-Gaia space used for the GMM. The best-fit model contains 7 Gaussian components shown with 2σ ellipses. As a result of the fit, each star in the sample has a probability of belonging to an individual Gaussian component. By comparing these probabilities, the component membership of each star is established. There can be cases that a star has comparable probabilities between multiple components, making the classification less certain. Although this can result in slight deviation between the GMM estimated weight and the number of stars assigned to each component, we note that this deviation is within a small fraction. For our samples, there are < 1% in APOGEE and < 2% in GALAH that the difference between the highest and the second-highest component probabilities for each star is less than 0.05 (probability ranges from 0 to 1). While we use the stars with high orbital eccentricity (e > 0.85) for our GMM study, the density histogram of the whole APOGEE-Gaia dataset is also shown as a 2-dimensional greyscale background for each projection to put the properties of the model Gaussian components into context. We also apply the XD algorithm to GALAH-Gaia, which include additional tracers. For the odd-Z and iron-peak elements, Na and Cu are included in addition to the abundances used in the fitting of the APOGEE-Gaia sample. For the s-process, Y and Ba are used instead of Ce. In the case of the r-process, Eu is included. This is particularly desirable, as no r-process tracer is available in the APOGEE-Gaia sample. Similar to the APOGEE-Gaia sample, the parameters of each of the components identified by the GMM in GALAH are listed in Table 2. Figures 3 and 4 shows sets of 2-dimensional projections among the 12-dimensional GALAH sample space.
The GMM for both samples finds strong evidence for 4 independent, physically plausible components. These comprise the three known constituents, together with an entirely new population called Eos. Only one of these components is of accreted origin, namely the GS/E. This dominates the eccentric sample locally. Importantly, there is no clear chemical evidence for another major accretion event in the Solar neighborhood. This makes sense -only the most massive accreted objects can fall deep into the potential well of the Galaxy. Innumerable minor accretion events have indeed occurred (e.g., Myeong et al. 2018a,b;Koppelman et al. 2019;O'Hare et al. 2020;Borsato et al. 2020;Limberg et al. 2021;Sofie Lövdal et al. 2022), but their residues are -relatively speaking -inconsequential for the local sample.
There are other reassuring properties of the GMM results that add credibility to its decomposition. First, there is consistency between the results when applied to APOGEE-Gaia and GALAH-Gaia in terms of the number and properties of the components. Both datasets provide evidence for 4 physically-motivated components. Secondly, whenever there are elements in common between the two samples, then there is consistency in the results for the 4 components, as shown for example in Figs. 2 and 3. Thirdly, there is also consistency between APOGEE-Gaia and GALAH-Gaia for the same groups The only substantial difference is in the fractional contributions of the 4 components. For example, the GMM judges that GS/E comprises a majority of the local halo when applied to APOGEE-Gaia, but only a third in GALAH-Gaia. Here, we readily concede that the relative weightings of the populations are due to the surveys' distinct selection functions. For example, the GALAH-Gaia sample is undernourished in stars below [Fe/H] < −1.3 dex, so the GMM over-weights the Splash, the Aurora and the Eos, while under-weighting the GS/E, as compared to APOGEE-Gaia.
Another interesting difference is that, in APOGEE-Gaia, the GMM is powerful enough to split the GS/E into two separate, but closely-related, Gaussians. These correspond to the α plateau and the knee of the progenitor dwarf. In APOGEE-Gaia, the GS/E population is so well represented, especially at low metallicities, that the GMM can even break down the population into its internal chemical structure.
Finally, in both the APOGEE-Gaia and GALAH-Gaia, the GMM finds low-weight "background populations". These are consistent with an outlier contribution, as they do not show physically plausible chemodynamical trends. For example, they have very large and thus physically implausible dispersion in most elements encompassing giant swaths of the abundance space. Since there will clearly be a discrepancy between the true distribution of each genuine Galactic population and the Gaussian model assumed by the GMM, residuals between the data and the model are expected. Such residuals across the sample space can easily be absorbed into a broad, low-weight and featureless background (Gaussian) density GMM component.
We now discuss the GMM results for the 4 physically plausible components in turn. For the GS/E and the Splash, this is largely a confirmation and amplification of earlier results, though now obtained by an unsupervised learning method. For the Aurora, the GMM identifies completely new chemical trends, e.g. in Eu, Y and Ba, whilst Eos is here identified and described for the very first time.

GS/E
The GS/E takes up ≈ 51% of eccentric stars in the APOGEE-Gaia sample (≈ 33% in the GALAH-Gaia sample). GS/E has an ex situ origin, having formed as a result of the disruption of a relatively massive galaxy. The mass of the progenitor dwarf has to be high to shrink and radialize the orbit (Vasiliev et al. 2022), so that the bulk of its stars are deposited into the inner regions of the Milky Way's halo. Only the most massive objects are able to sink sufficiently deep in the host's potential (e.g., dynamical friction being stronger for massive acquisitions) and contribute a sizeable fraction of their phase-mixed (and thus covering a significant volume) tidal debris in a given small region. Of all the components, the GS/E shows the highest mean orbital energy with the widest spread, as well as the lowest metallicity. This is also in accord with its ex situ origin and with its debris being spread during the interaction (see e.g., panels b) and c) of Figures 2 and 3.
With a progenitor mass of at least 1/10 of the Milky Way (e.g., Belokurov et al. 2018) at the time of merger, GS/E shows a metal-poor α-knee at [Fe/H] ≈ −1.3 (e.g., Helmi et al. 2018;Mackereth et al. 2019). Its inefficient star-formation is linked to its meagre enrichment in the odd-Z elements, which is reflected in the GS/E's low [Al/Fe] and [Na/Fe] level.
Since our APOGEE-Gaia sample extends down to lower metallicity, we can trace the GS/E's behavior at low metallicity better than in GALAH-Gaia. In Figure 2, we see that the GS/E population splits into two. The two components have similar mean energy. However, the low metallicity portion (GS/E 2 , yellow) shows a wider spread across the orbital energy range and also reaches down to lower orbital energies compared to its higher metallicity counterpart (GS/E 1 , green). Note that even in the absence of dynamical friction, the dwarf progenitor's stars could cover a wide range of energies with leading and trailing debris occupying preferentially lower and higher energies correspondingly (see also, D' Onghia et al. 2009, for complications in the case of a massive and rotating progenitor). This effect will blur any simple metallicity-energy correlation and could also result in the lower metallicity debris showing a wider spread, including the higher and lower ends of the orbital energy range. GS/E 1 and GS/E 2 correspond to (i) the knee/decline in α and in [Al/Fe] with increasing contribution from SN Ia and (ii) high-α plateau and slowly rising [Al/Fe]. The two GS/E components also show different behaviour in Ni, another element sensitive to SN Ia explosions.
In the APOGEE-Gaia data, the GS/E population covers a wide range of metallicities. So, we can expect its chemical distribution (e.g.,  Figure 2). GS/E's [Al/Fe] shows a mild increase along the metallicity range probed, reaching its peak around [Fe/H] ≈ −1.3, and turns to a decreasing trend beyond this metallicity. This pattern is in good agreement with the earlier analysis of APOGEE data reported by e.g., Hasselquist et al. (2021) and Belokurov & Kravtsov (2022).
The unbiased GMM decomposition also recovers the fact that the GS/E is high in [Eu/Fe] and [Eu/Mg], consistent with previous studies (Aguado et al. 2021a;Matsuno et al. 2021). Eu and Mg are both good tracers for Core-Collapse Supernovae (CCSNe) contribution in star formation history as they are believed to be produced from massive stars most exclusively. But, in case of Eu (r-process element), Neutron Star Mergers (NSMs) are another considerable source of enrichment (see e.g., Kobayashi et al. 2020). These two r-process channels have distinct timescales -CCSNe (NSMs) are known for prompt (delayed) r-process production, respectively. This means for a system with longer duration of star formation the [Eu/Mg] level will increase, as it will give more chance for the delayed Eu enrichment of NSMs to be reflected. Thus Eu trends mentioned above are a consequence of the longer star formation duration for the GS/E than the Milky Way. The Milky Way progenitor's has a much more vigorous and efficient star formation than GS/E progenitor.

Splash
The Splash is another major component in our higheccentricity sample. Compared to GS/E, the orbital energy of the Splash stars is also low. This is consistent with the idea that it is an in situ population kicked out of the Galactic disk by the GS/E merger.
The Splash contributes ≈ 25% in APOGEE-Gaia and ≈ 31% in GALAH-Gaia sample. The differences between the two datasets in the fractional contribution of this population are likely due to the differences in the surveys' selection functions, particularly the metallicity range probed. For example, our GALAH-Gaia sample has very few stars below [Fe/H]< −1.3 dex. Since the GS/E population is dominant below this metallicity, GS/E's total weight in the GALAH sample appears to be lower compared to the APOGEE sample which extends down to lower metallicity.
In The main advance in this paper is here the chemical properties of the Splash have been "discovered" by an unbiased GMM decomposition rather than extrapolated from hand-made cuts.
There is remarkable consistency across all chemical dimensions between the Splash and the high-α disc. This is of course only to be expected -the Splash is born out of the evolved high-α, in situ component of the Galaxy.

Aurora
The local eccentric halo component with the lowest fractional contribution (marked in red, at ≈ 6% in APOGEE and ≈ 11% in GALAH) corresponds to the recently discovered Aurora population (see Belokurov & Kravtsov 2022). This in situ population is identified by the GMM primarily based on the chemical information, yet it shows behaviour entirely consistent with that described in Belokurov & Kravtsov (2022) where it was picked up based on kinematic evidence. For example, according to the GMM, in APOGEE DR17 Aurora is limited to the metallicity range of −1.5 <[Fe/H]< −0.9 dex in line with the earlier study. Aurora has one of the lowest mean orbital energy, similar to Splash. However, given the apo-centre selection cuts imposed in this analysis, the main body of Aurora likely has even lower energy and probably resides mostly in the inner Galaxy.
The unbiased GMM decomposition of both APOGEE and GALAH data offers a new, expanded chemical view of Aurora, demonstrating clearly where and how this population stands out from the rest of the Milky Way. For example, at its high metallicity end, i. Curiously, enhancement in the s-process elements has been seen in a number of currently existing GCs. For example, Mészáros et al. (2020) demonstrated the existence of s-process rich stars in a number of Galactic GCs, together with a strong correlation between Ce and metallicity. Similarly, Marino et al. (2009) reported bimodality in the s-process elements including Y and Ba, also accompanied by a strong correlation with metallicity, which they relate to the multiple cluster populations. These studies suggested that the s-process enhancement takes place in low-and intermediate-mass Asymptotic Giant Branch (AGB) stars. Belokurov & Kravtsov (2022) argued that Aurora's high dispersion and anomalous abundance in N, Al, O, Si could be linked to a significant contribution from massive star clusters (some of which may even survive to become Galactic GCs today). Here we provide additional evidence to support this hypothesis: the anomalous enrichment in odd-z and s-process species together with metallicity correlations.
As far as the iron-peak elements are concerned, in particular those sensitive to SN Ia enrichment, as shown in Figure 5,  Comparing Aurora to Heracles (c.f. this work and those by Naidu et al. 2022;Horta et al. 2022), it appears that in all chemical dimensions these two populations are close to indistinguishable. Yet, the proposed origins of the two are very different. While Aurora is created in situ, Heracles is hypothesised to be an ancient merger event that occurred before the Milky Way-GS/E interaction. Heracles' stellar debris typically occupy higheccentricity orbits (eccentricity > 0.6) but at a considerably lower energy range than the GS/E stars. Horta et al. (2021) used the apparent energy gap seen in the APOGEE dataset (see e.g., their Figure 4) to separate the GS/E and the Heracles stars. The Heracles' chemical signature is also different from that of the GS/E, for example, its [Mg/Fe] level remains consistently high while at the same range of metallicity, the GS/E shows a clear α-knee, i.e. a decreasing trend with [Fe/H]. Additionally, there have been two separate studies proposing the existence of another major merger in our Galaxy, which supposedly took place, like Heracles, before the GS/E event. Based on the age-metallicity and the number of globular clusters residing in the very centre of the Milky Way, Kruijssen et al. (2019) argued for a major merger with a massive galaxy they named Kraken; a similar merger (Koala) was suggested by Forbes (2020). As these studies all suggest an older massive merger event (occurring before GS/E) with its debris mostly populating the inner Galaxy, a consensus is starting to emerge in the literature that Kraken, Koala and Heracles are in fact the same event (see e.g., Naidu et al. 2022).
A recent study by Lane et al. (2022) (see e.g., their Section 6.3) pointed out that the dynamical characteristics of Heracles (e.g., the energy gap used to identify it as a distinct low-energy overdensity), is similar to the artificial density features induced by the APOGEE selection function. Indeed, no evidence for such an energy gap exists in the GALAH data which has a different (less bi-modal) selection function.  (Horta et al. 2021;Naidu et al. 2022;Horta et al. 2022). A more in-depth detailed analysis can be carried out by comparing the Aurora chemical properties as displayed in Figures 2, 3, 4 and 5 here to Heracles' chemical fingerprint as shown in Horta et al. (2022). The similarities are striking. Despite being nearly identical chemically, Heracles and Aurora appear to occupy very different regions of the orbital energy (and consequently Galactocentric radius) distribution. Heracles is selected to have low energy and reside in the bulge, while Aurora has been identified -both here and in Belokurov & Kravtsov (2022) -outside of the bulge, amongst stars on high energy orbits. A comparison of the chemo-dynamical properties of Aurora and Heracles seem to indicate that this is the same population that is not bound to either inner or outer Galaxy. Instead, probably, it extends beyond the Solar neighborhood but is mostly concentrated in the bulge, similar to the numerical examples shown in Figure 12 of Belokurov & Kravtsov (2022).

Eos
The GMM reconstruction of the local halo populations has revealed a new eccentric component that we dub Eos. This previously unseen population occupies the metallicity range similar to that of Splash, i.e. −1 <[Fe/H]< −0.3 dex. Eos, which is shown in blue in the Figures in this paper, contributes ≈ 10% in APOGEE sample and ≈ 15% in GALAH sample 5 . What sets Eos apart is its relatively low α-abundance level, it appears to be intermediate between the high-α sequence (in our sample represented by Splash) and the GS/E at similar metallicity. Note also the difference in the orbital energy between Eos and the high-α components (Aurora, Splash). Eos stars have relatively higher orbital energy implying that it populates larger (outer) Galactic radii.
At its highest metallicity, across all elements Eos connects seamlessly to the low-α (thin disk) sequence, yet unlike the disk it is highly eccentric. On the other hand, at its lowest metallicity Eos resembles more GS/E than any other Galactic stellar populations. This is particularly obvious in [ (2010) and Hawkins et al. (2015) and conclude that Eos was created in situ. As shown by Hasselquist et al. (2021) 2021) as a transition region connecting the old thick (high-α) and young thin (low-α) disks, with a possible age gradient (older at higher [α/Fe], younger at lower [α/Fe]). This Bridge region doesn't exactly overlap with Eos's possible higher-α subbranch (Bridge region covers slightly higher metallicity and spans a very broad range of metallicity; see e.g., Figure 6 of Ciucȃ et al. 2021). Yet, the general idea of the transition signature running from higher-[α/Fe] with higher-[Fe/H] to lower-[α/Fe] with lower-5 Lists of member stars are available electronically from the authors [Fe/H] could be in line with Eos's possible higher-[α/Fe] subbranch feature, as a trace of gradual gas mixing (polluting) between the leftover gas from high-α (think) disk and the gas from GS/E merger (more discussion on Eos's origin in Section 5.1). We should also consider the possible blending from GMM fitting between Eos and the high-α disk populations, as these higher-[α/Fe] Eos stars are relatively closer to the high-α disk population in the chemical space.
Eos's low-α abundance combined with its intermediate odd-Z abundance level also places Eos stars in a unique region of [Mg/Mn]-[Al/Fe] and [Mg/Cu]-[Na/Fe] planes (panel f) of Figures 2 and 3, panel c) of Figure 4). These planes have been successfully used to discern the less-evolved (and/or ex situ) and more-evolved (in situ) stellar populations. As a less evolved portion of the low-α in situ population, Eos appears between the region occupied by GS/E and the bulk of the low-α in situ stars (shown as a 2D density histogram in the Figures). Eos also sits below the region where Aurora is located in this plane emphasising a clear difference in their α abundance levels. Although both Eos and Aurora are likely to be the Galaxy's in situ components, clear differences between them in a number of key elemental trends such as metallicity, α, and odd-Z suggests different formation channels for Eos and Aurora. In several of the chemical projections Eos also displays a correlation with metallicity, i. The chemical trajectory of Eos is unique, its chemical evolution appears to have started from the gas polluted by the GS/E and evolved into a typical low-alpha (thin disk) population.

Clues to Eos's origin
Deciphering the Eos's origin may help us to solve the conundrum of the (trans)formation of the Milky Way's disc. A part of the puzzle is the existence of a striking chemical bi-modality: two nearly parallel but well separated sequences are detected in the disc's α-[Fe/H] plane at a wide range of stellar metallicities and Galactocentric distances (e.g. Fuhrmann 1998;Feltzing et al. 2003;Adibekyan et al. 2013;Bensby et al. 2014;Hayden et al. 2015). The high-α sequence and the low-α sequence have been linked to the so-called "thick" and "thin" discs respectively, although given the recent structural analysis it may be more appropriate to label the discs "inner" and "outer" (see e.g. Haywood et al. 2013;Rix & Bovy 2013). Chronologically, the high-α and the lowα discs appear distinct: the former built up most of its mass before 8-10 Gyr ago, the epoch around which the first cinders of star-formation have been detected in the thin disc which has been growing steadily ever since (see e.g. Haywood et al. 2018a;Fantin et al. 2019).
Several scenarios have been proposed to explain the α-dichotomy in the Galactic disc. For example, the "two-infall" model requires that in the distant past an episode of gas accretion diluted the Milky Way's reservoir, resetting its metallicity so that the subsequent star-formation could proceed along the lower α-sequence starting from a lowered value of [Fe/H] (see Chiappini et al. 1997;Chiappini 2009;Grisoni et al. 2017). Alternatively, the observed α-bimodality can be produced via radial migration and "churning" of the stars created along the two α-sequences in distinct (broadly speaking, inner and outer) parts of the disc (see Sellwood & Binney 2002;Roškar et al. 2008;Schönrich & Binney 2009;Minchev et al. 2013). A new idea was proposed recently where the two spatially overlapping α-sequences can start forming contemporaneously: massive dense gas clumps dominate the production of stars with high αabundance, while the smooth distributed star-formation is responsible for the low α-abundance (Clarke et al. 2019). Finally, it is now also possible to form galaxies with α-bimodal discs in Cosmological hydrodynamical simulations. Curiously, in the simulated discs, the α-dichotomy is generally linked to ancient massive gasrich accretion events (see e.g. Brook et al. 2012;Grand et al. 2018;Mackereth et al. 2018;Buck 2020), however the details might vary. For example, in Agertz et al. (2021) around redshift z=1.5 a gas-rich accretion creates an outer low-α disc misaligned with respect to the earlier-formed inner high-α one; after a while, the two discs line up, leading to spatially overlapping low-and high-α stellar populations (see Renaud et al. 2021). In Grand et al. (2020), an analogue of the thick disc exists prior to the GS/E event but this violent interaction contributes to its additional heating; moreover, the gas-rich GS/E merger induces a star-burst in the central MW which adds between 20% and 60% stellar mass to the thick disc. As hinted by these recent studies based on cosmological hydrodynamical simulations, the Galactic discs and the GS/E event are linked in terms of their kinematics, chemical composition, and chronology.
Curiously, both Renaud et al. (2021) and Grand et al. (2020) discuss the emergence of an in situ halo-like population distinct from Splash and not much dissimilar to Eos around the time of the GS/E accretion. In Grand et al. (2020), a GS/E-like merger both splashes the preexisting disk and induces a starburst in the centre of the Galaxy which leads to creation of a puffed-up halo-like in situ component. The simulated Splash would be expected to have a wide range of ages truncated around the time of the GS/E arrival, while the starburst population's SFH would peak at the epoch of the interaction. In Renaud et al. (2021), the GS/E-like merger not only triggers star formation in the mis-aligned outer gaseous disk but also kick-starts its re-orienting with the preexisting inner one. The two gas disks align rapidly but the first stars that were born in the outer disk are left behind in a strongly tilted configuration. The early outer disk stars subsequently phase-mix, and at the present day lack any appreciable angular momentum -the memory of their particular initial condition. The scenario of Renaud et al. (2021) is remarkable in that it predicts the existence of an in situ halo-like population with a narrow age distribution and the chemistry similar to the metal-poor low-α disk.
In the two-infall scenario, a variation of which is played out in the simulations discussed above, the resetting of the disk's chemical sequence is facilitated by the accretion of fresh, poorly-enriched gas from the medium around the host. It is also possible to directly invoke the accretion and recycling of the gas from the GS/E progenitor (see Grand et al. 2020;Palla et al. 2020). The dwarf is not required to donate its gas (which may not be enough to make a huge difference), it merely needs to enrich the gas surrounding the host. The exact amount of the heavy elements the dwarf can infuse into the condensing circumgalactic medium (CGM) will depend on the state of the gas and the parameters of the host-satellite interaction. Importantly, however, the mass loading factor of the pre-enriched outflows are predicted to be higher in the dwarf galaxy compared to the more massive host (Mitchell et al. 2020;Pandya et al. 2021). This is because the energetics at the feedback site do not depend on the host's global properties, but the restoring force (gravity) does. This idea may have some support in the chemical properties of Eos. As Figures 2, 3 and 4 demonstrate, at its lowest metallicity, Eos is closest to the metal-rich portion of the GS/E sequence.

Caveats
Our study uses GMM to analyse unlabelled data in hyper-dimensional space consisting primarily of chemical information. Our modelling is based on the assumption that any sub-population found in the data follows a Gaussian distribution in the data space of choice. Unfortunately there is no guarantee that this is actually true. If a real physical trend in the data deviates severely from Gaussian distribution (e.g., a stellar population in [α/Fe]-[Fe/H] with changing α abundance slope), this single stellar population can be segmented into multiple Gaussian components in the GMM fit. Any hard boundary cut applied to the data can also degrade the GMM performance. Therefore, we avoid the use of any dimension in our GMM if this dimension is directly used for any selection cuts adopted.
Populations seen in our sample are the higheccentricity tails of their actual orbital distribution. For each population, this will bias the weight of the population assigned by the GMM depending on the degree of radial anisotropy. For example, Aurora is shown to be isotropic with its velocity dispersion, i.e., σ R ≈ σ z ≈ σ φ (see Belokurov & Kravtsov 2022, for more detail). In contrast, GS/E or Splash are known for their strong radial anisotropy (e.g., Belokurov et al. 2018;Myeong et al. 2018c,d;Belokurov et al. 2020). Such a difference is reflected in their orbital eccentricity distribution, and different proportions of each population will be included in our high-eccentricity sample. While we have avoided making a cut on energy, the apocenter > 5 kpc cut puts a slightly softer lower limit on the orbital energy range in our sample. This can bias the components' properties if they are found near this lower energy boundary. For example, Aurora and Splash appear to extend down to the lowest energies, close to the limit imposed.
The reliability of chemical measurements varies with all of the atmospheric parameters and generally deteriorates towards the lower end of the metallicity range. To minimize the effect of varying uncertainty, we adopted the XD algorithm for our GMM which can take account of measurement errors associated with the observed data. However, of course, the systematic (unreported) errors will still affect the results of our modelling. The metallicity probed by both APOGEE and GALAH is truncated earlier as both surveys struggle in the low-[Fe/H] regime. The metallicity limit can also confuse the recovery of the true stellar populations in the data, especially for those components that extend far down in [Fe/H], e.g. GS/E and Aurora.

Summary
Our study provides one of the first attempt to assemble an unbiased census of the main halo components in the vicinity of the Sun. To carry out this experiment, we use the Extreme Deconvolution implementation of the Gaussian Mixture Model. To help the GMM concentrate on the stellar populations of interest, we apply an eccentricity cut thus limiting the sample to the halo stars only. The GMM model provides results that show consistency between APOGEE and GALAH in terms of the number and the properties of the components recovered. Moreover, when comparing the detailed properties of the individual components in APOGEE and GALAH we see clear consistency for the same chemical elements and for the same groups of elements (e.g. the [Al/Fe] behaviour is consistent with the [Na/Fe] behavior). Small differences in fractional contributions of the components between APOGEE and GALAH are likely due to the surveys' distinct selection functions.
The reliability of the GMM decomposition reported here can be established by comparing the chemical properties of the known and well-studied halo components to those reported in the literature. For example, our view of the GS/E chemical fingerprint is in full agreement with the recent studies (see e.g. Das et al. 2020;Hasselquist et al. 2021;Buder et al. 2021b;Aguado et al. 2021a;Matsuno et al. 2021;Naidu et al. 2022;Horta et al. 2022;Carrillo et al. 2022). Similarly, we demonstrate that across all chemical dimensions investigated, Splash is indistinguishable from the high-α (thick disk) of the Milky Way. This is of course in agreement with the proposed origin of this in situ halo population (see e.g Bonaca et al. 2017;Haywood et al. 2018b;Gallart et al. 2019;Di Matteo et al. 2019;Belokurov et al. 2020).
While it is reassuring to see agreement between the earlier attempts to decipher the chemical behavior of the halo denizens and our unsupervised analysis, the results of this work are not limited to mere confirmation. Let us summarise the most interesting new findings here.
• The GMM finds only four independent halo components dominating the Solar neighborhood, three of these are previously known, i.e. Aurora, Splash and GS/E, and one is new. There is strong evidence that the newly discovered population Eos is of in situ origin. Therefore, in the Solar neighbourhood, our analysis detects no evidence for tidal debris from any additional large accretion event.
• As reported by the GMM, the multi-dimensional chemical view of the ancient in situ component Aurora (see Belokurov & Kravtsov 2022) agrees remarkably well with the chemical properties of the alleged merger remnant known as Heracles (Horta et al. 2021;Naidu et al. 2022;Horta et al. 2022). The only noticeable difference appears to stem from the fact that in the literature the Heracles stars are identified using a hard-cut in [Al/Fe] while here we dispense with selections in the chemical space. Our study preserves the integrity of the stellar populations in the elemental abundance space and demonstrates that in [Al/Fe]-[Fe/H] dimension Aurora evolves rapidly from low [Al/Fe] to anomalously high [Al/Fe]. This behaviour is predicted for a rapidly star-forming self-enriching galaxy with a mass similar to that of the Milky Way (see Horta et al. 2021) and thus obviates the necessity for introducing an additional merger event to explain this stellar population.
• The chemical properties of the discovered in situ born structure Eos evolve significantly as a function of metallicity. We speculate that Eos appears to originate from the gas polluted by the GS/E and evolved to resemble the (outer) thin disk of the Milky Way.
The unsupervised view of the local halo the GMM provides also allows us to decipher many new and fine details of the four components detected. For example, in APOGEE, the GS/E is split into two components that primarily differ in metallicity but trace distinct epochs in the life of the progenitor galaxy, i.e. one corresponding to the high-α plateau and slowly rising [Al/Fe] and another to the knee/decline in α and in [Al/Fe] due to the increased SN Ia contribution. We also see that the two GS/E components behave differently in Ni (ironpeak element particularly sensitive to variations in SN Ia yields  Naidu et al. (2021).
Deciphering the origin of Eos with the help of numerical simulations, both those currently available such as FIRE (Hopkins et al. 2014), Auriga (Grand et al. 2017), VINTERGATAN ) and ARTEMIS (Font et al. 2020), as well as those upcoming, is the obvious next step. Currently the feedback physics remains somewhat unconstrained but there may be routes to test the proposed formation mechanism of Eos, in particular, the role of the GS/E progenitor galaxy in the pollution of the circum-host gas and the triggering of the in situ star formation. The (relative) stellar ages (see e.g. Conroy et al. 2022;Xiang & Rix 2022) also provide valuable information and will help to expand the current analysis to a chrono-chemodynamical study. Undoubtedly, detailed information from the upcoming large spectroscopic surveys (e.g., DESI, 4MOST, WEAVE, SDSS-V) and astrometric surveys (e.g., Gaia, LSST) will open up new arenas for decoding the genealogy of the Galaxy.