Chronology of the chemical enrichment of the old Galactic stellar populations

The Milky Way accreted several smaller satellite galaxies in its history. These mergers added stars and gas to the Galaxy and affected the properties of the pre-existing stellar populations. Stellar chemical abundances and ages are needed to establish the chronological order of events that occur before, during, and after such mergers. We report precise ages ($\sim$6.5%) and chemical abundances for the Titans, a sample of old metal-poor dwarfs and subgiants with accurate atmospheric parameters. We also obtain ages with an average precision of 10% for a selected sample of dwarf stars from the GALAH survey. We used these stars, located within $\sim$1 kiloparsec of the Sun, to analyse the chronology of the chemical evolution of in-situ and accreted metal-poor stellar populations. We determined ages by isochrone fitting. For the Titans, we determined abundances of Mg, Si, Ca, Ti, Ni, Ba, and Eu using spectrum synthesis. The [Mg/Fe] abundances of the GALAH stars were re-scaled to be consistent with the abundances of the Titans. We separated stellar populations by primarily employing chemical abundances and orbits. We find that star formation in the so-called Gaia-Enceladus or Gaia-Sausage galaxy, the last major system to merge with the Milky Way, lasted at least 3 billion years and got truncated 9.6 $\pm$ 0.2 billion years ago. This marks with very high precision the last stage of its merging process. We also identified stars of a heated metal-poor in-situ population with virtually null net rotation, probably disturbed by several of the early Milky Way mergers. We show that this population is more metal rich than Gaia-Enceladus at any time. The sequence of events uncovered in our analysis supports the hypothesis that Gaia-Enceladus truncated the formation of the high-$\alpha$ disc and caused the gas infall that forms the low-$\alpha$ disc, in agreement with theoretical predictions.


Introduction
According to standard cosmology, large galaxies form hierarchically by merging smaller stellar systems (Blumenthal et al. 1984;Springel et al. 2005).The identification of Milky Way merger remnants has enjoyed great progress in recent years thanks to the combination of the astrometric data from the Gaia mission (Gaia Collaboration 2016) with photometric and spectroscopic data from large stellar surveys such as the Sloan Digital Sky Survey (SDSS, Abazajian et al. 2009), the Apache Point Observatory Galactic Evolution Experiment (APOGEE, Majewski et al. 2017), and the Galactic Archaeology with HERMES survey (GALAH, De Silva et al. 2015).Such investigations provide evidence of the early assembly of the Milky Way halo from satellite galaxies (e.g., Naidu et al. 2020;An & Beers 2021;Shank et al. 2022) and demonstrate that the Milky Way protodisc formed very early on (Di Matteo et al. 2020;Sestito et al. 2020;Venn et al. 2020;Kielty et al. 2021).Current evidence indicates that the last major merger to have occurred in Chemical abundances and updated ages of the Titans, as well as the data used in Sect.7 are only 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/673/A18the Milky Way was with the so-called Gaia-Enceladus or Gaia-Sausage galaxy (Helmi et al. 2018;Belokurov et al. 2018).
Early mergers between the Milky Way and smaller galaxies not only entailed the accretion of stars, but also of gas.Models have shown that this stripped gas has an influence on the star formation history of the accreting galaxy (Vincenzo et al. 2019;Belokurov et al. 2020;Grand et al. 2020).For the Milky Way, it has been shown that a model with two phases of gas infall (Chiappini et al. 1997;Micali et al. 2013) is able to explain the distinct chemical evolution of the thick and thin discs.Thindisc stars are generally younger and have lower [α/Fe]1 ratios compared to thick-disc stars (Fuhrmann 2011;Haywood et al. 2013;Recio-Blanco et al. 2014).The stellar age-metallicity distribution of the two populations splits at an age of about 8 Gyr (Xiang & Rix 2022).
The delay between the infall that forms the thick disc and the next that forms the thin disc has been estimated to be between 4.5 and 5.5 Gyr (Spitoni et al. 2019(Spitoni et al. , 2020)).It has been conjectured that gas from a merger, such as that with Gaia-Enceladus, could have enabled the formation of the low-α thin disc (Brook et al. 2007;Bonaca et al. 2020), thus being the source of the second infall gas.In this sense, the delay in star formation between the discs can be understood within a gas shock-heating theory (Noguchi 2018).Quenching of the previous episode of star formation could also be related to the merger of Gaia-Enceladus (Vincenzo et al. 2019).Furthermore, the Gaia-Enceladus merger also heated up a proto-disc stellar population, creating part of the population that we now call the inner stellar halo (Brook et al. 2003;Bonaca et al. 2017;Di Matteo et al. 2019;Belokurov et al. 2020).
The evidence mentioned above indicates that the Gaia-Enceladus merger played a major role in shaping the configuration of the current stellar populations of the Milky Way, although some uncertainty remains with regard to the chronology of all related events (see, however, Ciucȃ et al. 2021Ciucȃ et al. , 2023, for significant improvements in this regard).To determine the relationship between all of these events, it is particularly critical to obtain a precise timing of the last stage of the merger.Previous attempts have obtained estimates between 8 and 10 billion years ago2 (Helmi et al. 2018;Gallart et al. 2019;Belokurov et al. 2020;Montalbán et al. 2021;Borre et al. 2022).
In this work, we precisely date the chemical enrichment of old Galactic stellar populations taking advantage of the Titans (Giribaldi et al. 2021).This is a sample of 48 relatively nearby dwarf and subgiant stars with metallicities, [Fe/H], between −3.1 and −0.8 dex.More importantly, we have previously determined values of their effective temperatures (T eff ) and surface gravities (log g) with accuracy below 1% (∼50 K and ∼0.04 dex, in T eff and log g, respectively).This accuracy is unmatched by any other sample of metal-poor stars and this is what enables our determination of highly precise ages.To complement the analysis, we selected a second sample of similarly unevolved stars from the Data Release 3 (DR3) of the GALactic Archaeology with HERMES (GALAH) survey (Buder et al. 2021).
Stellar ages were determined by isochrone fitting.Thanks to the accurate parameters, the derived ages have an average precision of 6.5% (900 Myr or 0.9 Gyr) for 45 of the 48 Titans.For the GALAH turnoff stars, the ages have an average precision of 10% (1 Gyr).The precision of these ages enables a robust discussion of the chronology of the events in the Galaxy that can be related to the Gaia-Enceladus merger.
This paper is organised as follows.The data, sample, and stellar orbit computations are described in Sect. 2. Section 3 describes how the parameter scale of the GALAH sample was brought in agreement with the accurate scale defined by the Titans.The estimation of stellar ages is described in Sect. 4. Section 5 describes the determination of the chemical abundances for the Titans and the correction of the Mg abundances of the GALAH sample.Section 6 describes the procedure followed to separate stellar populations.Finally, in Sect.7 we interpret our results and discuss their implications in the context of the Galaxy evolution.

Data, sample, and stellar orbits
The spectroscopic data used to determine the atmospheric parameters and abundances of the Titans are described in Giribaldi et al. (2021).In short, the observations were obtained with the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000), at the Very Large Telescope of the European Southern Observatory (ESO), and are publicly available at the ESO data archive3 .The spectra include the Hα line at 6562.797 Å, have resolving power R ≥ 40 000, and a signalto-noise ratio (S/N) greater than 100.The sample selection described in Giribaldi et al. (2021) resulted in 41 stars, for which the values of T eff and log g were estimated with 1% accuracy.In addition to those stars, in Giribaldi et al. (2021) we also provided accurate parameters for seven metal-poor Gaia benchmark stars (Hawkins et al. 2016;Jofré et al. 2018) which were analysed using the same methods.The combination of these two samples is what we discuss in this work and is the sample that we generically refer to as the 48 Titans.
To increase the sample of stars for age determination, additional metal-poor stars were selected from data release three (DR3) of the GALAH survey (Buder et al. 2021).We first selected stars with [Fe/H] < −0.5, 2 < log g < 5, T eff > 4000 K, and high quality-atmospheric parameters (cannon_flag = 0).We excluded stars that are too cool and/or too bright because their analysis is usually more uncertain.More metal-rich stars were excluded to avoid strong contamination of the sample by thindisc stars.The sample was restricted to stars with Gaia parallax ≥ 0.2 mas that have errors smaller than 20% (Lindegren et al. 2021b), resulting in 12 405 stars.We further removed stars with a total4 Galactic velocity below +210 km s −1 , to avoid stars with disc kinematics, reducing the sample to 6250 stars.This initial sample was used to train the algorithm that separates stellar populations (Sect.6).For the main discussion, ages and chemical abundances are used only for dwarf and turnoff stars with log g in the range of 3.5 < log g < 4.5.This is the range where precise ages can be determined.
To find the optimal way to separate the populations, in the first step of the analysis, we kept in the sample the stars with unreliable abundances of α elements from the GALAH sample (i.e., flag_alpha_fe > 0); this is the sample of 6250 stars mentioned above.This choice was motivated by the need to increase the statistics and the density of points in the multidimensional space where the populations are separated.Applying the flag on α abundances reduces the sample to about 1500 stars.In any case, as we discuss below, we verify that the final separation in the chemical diagram of [α/Fe] vs. [Fe/H] is consistent with the usual findings in the literature.Accreted stars are found to have low-[α/Fe] ratios down to a metallicity around [Fe/H] = −1.5.For lower metallicities, high-and low-α sequences merge.Furthermore, when we select only stars with the best ages (Sect.4), the sample remains only with those stars that have reliable abundances, that is, those with flag_alpha_fe = 0.This final sample has 208 stars: 25 Titans and 183 from GALAH DR3.These stars are presented in Fig. 1, where the in situ populations, Splash (defined in Belokurov et al. 2020) and Erebus5 (defined in this work, see Sect.6), and the accreted population of Gaia-Enceladus are distinguished.
For the GALAH sample, we adopted the kinematic and dynamic parameters provided in one of the value-added catalogues part of DR3.We refer to Sect. 7.3.3 of Buder et al. (2021) for details.For the Titans, we followed the same method to compute the orbits.We used the heliocentric velocities and distances given in Giribaldi et al. (2021).Orbits were computed with the galpy package (Bovy 2015)   potential of McMillan (2017).The orbits were integrated for 13.8 Gyr.For the solar Galactocentric distance and the circular speed at the Sun's position, we adopted the values determined in McMillan (2017).For the solar peculiar velocity, we adopted the values given from Schönrich et al. (2010).

Accuracy and consistency of the parameter scales
Atmospheric parameters (T eff , log g, and [Fe/H]) of the Titans are those determined in Giribaldi et al. (2021).For the determination of T eff , the method relies on fitting the observations with synthetic Hα line profiles computed using hydrodynamic three-dimensional (3D) model atmospheres and taking into account departures from the local thermodynamic equilibrium (Amarsi et al. 2018).As demonstrated in Giribaldi et al. (2021), an accuracy of the order of 1% was reached for the T eff values.To determine log g, Gaia parallaxes of the EDR3 (Lindegren et al. 2021b), corrected for the biases discussed in Lindegren et al. (2021a), were used.A similar accuracy level of ∼1% was achieved.
In Giribaldi et al. (2021), it was shown that the accurate T eff scale of the Titans is compatible with the scale set by the InfraRed Flux Method (IRFM) implementation in Casagrande et al. (2010).New IRFM relations using photometry in the Gaia system were recently published by Casagrande et al. (2021), on a scale that was made compatible with that of the previous work.Therefore, for our selected sample of GALAH stars, Furthermore, the values of log g of the GALAH stars were redetermined here, using the new IRFM T eff values, to be consistent with the method used in Giribaldi et al. (2021), namely log g was determined by fitting Yonsey-Yale isochrones (Kim et al. 2002;Yi et al. 2003) using the code q 2 (Ramírez et al. 2013).As input, the values of T eff , [Fe/H], V magnitude, and parallax ( ) are required.In this case, [Fe/H] and V values compiled in the GALAH DR3 catalogue were used.Metallicity values were not adjusted because they appear to have negligible offsets compared to the scale of the Gaia benchmark stars (Jofré et al. 2018;see Fig. 6 in Buder et al. 2021).The Gaia benchmark stars metallicity scale is already consistent with the metallicity values of the Titans (see Giribaldi et al. 2021).Our stellar ages and the element abundances required for the separation of the populations are based on the scales of T eff and log g determined here.

Stellar ages
Stellar ages for the Titans and the GALAH sample were determined using Yonsey-Yale isochrones (Kim et al. 2002;Yi et al. 2003) and the code q 2 (Ramírez et al. 2013).As mentioned above, the code requires T eff , [Fe/H], V, and , and their corresponding errors as input.The code uses the value to estimate the absolute magnitude M V using a distance value that comes from a simple parallax inversion.Here, we decided to adopt for our stars the geometric distances from Bailer-Jones et al. ( 2021), which have already been corrected for the parallax zero-point offsets.As the code needs parallaxes as input, we simply inverted the distances from Bailer-Jones et al. (2021).As parallax errors, those directly given in the Gaia catalogue were adopted; for our nearby stars, they are virtually equivalent to those obtained by inverting the 16th and 84th percentiles of the geometric distance posterior.For the GALAH sample, V magnitudes were corrected from extinction, using −3.1 × E(B − V), and adopting reddening values from the GALAH DR3 catalogue.These A18, page 3 of 16 values were estimated as described in Casagrande et al. (2021).The GALAH observations mostly avoid the Galactic plane (see Fig. 27 in Buder et al. 2018), which implies small reddening values for most of the sample.In fact, for 85% of the stars that we selected, the E(B − V) values are distributed between 0 and 0.1 mag, with a peak at 0.04 mag.In addition, for 10% of the stars, the E(B − V) values vary between 0.1 and 0.2 mag while the remaining 5% have E(B − V) > 0.2 mag.For these 5% (which correspond to eight stars), we compared our adopted reddening values with estimates made using the database Stilism6 (Capitanio et al. 2017).We found a systematic offset of about ∼0.1 mag that would in turn increase the age of these few stars by about 1 Gyr.We note here that we assumed ±0.08 mag as the typical combined value of extinction and magnitude uncertainty, thus the effect of such offsets is essentially already included in the error we estimated for the ages.We decided to still use E(B − V) from GALAH for these eight stars, to ensure consistency with the rest of the sample.
Turn-off stars were selected for the determination of age because the spacing between the isochrones in this region allows for more precise values to be estimated compared to any other place on the HR diagram (see Fig. 1).The most important factor that limits the precision in isochronal ages for single stars is the error in T eff ; in this work, the average accuracy in T eff is 50 K, and this is unmatched in any other stellar sample.Another important factor is the determination of log g or luminosity (closely related quantities).To that end, precise parallaxes are needed, such as those that have become available thanks to Gaia EDR3.All these conditions have allowed us to derive ages with an average precision of 6.5% (0.9 Gyr) for 45 of the 48 stars in the Titans sample and with an average precision of 10% (1 Gyr) for 183 stars of the GALAH catalogue.For stars in the Splash and Erebus populations, we accepted ages that are slightly worse to increase their numbers.In such cases, the errors can be up to 1.5 and 2.0 Gyr for Splash and Erebus, respectively.

Chemical abundances
To derive chemical abundances for the Titans, we adapted iSpec (Blanco-Cuaresma et al. 2014) Python subroutines to be able to perform a supervised line-by-line fitting.Synthetic spectra calculations were performed with the Turbospectrum radiative transfer code (Plez 2012) using the MARCS (Gustafsson et al. 2008) model atmospheres.Abundances were derived keeping all stellar parameters fixed to the values obtained in Giribaldi et al. (2021): T eff , log g, [Fe/H], resolution, projected rotational velocity (vsin i), microturbulence velocity (v mic ), and macroturbulence velocity (v mac ).
The abundances of Mg include corrections for atomic diffusion, which we obtained by interpolating the values of Korn et al. (2007) and Gruyters et al. (2013).These corrections only affect stars with metallicities below [Fe/H] = −1.65 dex.At higher values of [Fe/H], the corrections for Mg cancel out with those for Fe.For stars with [Fe/H] below −2.10, the correction is of +0.1 in [Mg/Fe] and [Si/Fe].The correction decreases linearly between [Fe/H] = −2.10 and −1.65,where they become zero.

Non-LTE corrections
Our abundances are computed using one-dimensional (1D) model atmospheres and assuming local thermodynamic equilibrium (LTE).However, in particular for metal-poor stars, deviations from LTE (the so-called non-LTE effects) have been shown to be of importance.Therefore, when possible, we used literature calculations of non-LTE effects to correct the abundance that we derived.
Regarding the abundances of Mg, it was found in Mashonkina (2013) that, for turnoff stars (T eff = 6000 K and log g = 4 dex) with metallicities of [Fe/H] = −1, −2, and −3 dex, a constant correction in the abundances of ∼+0.05 is needed for the Mg lines analysed here.This is consistent with the fact that we find agreement between the abundances derived from these lines, as we can see in Fig. 3.The most extreme outliers in this figure are among the coolest and most metal-rich Titans ([Fe/H] > −1.2 dex) with the highest log g values.Stars in this range of surface gravity were not considered for the main discussion in the Titans or GALAH samples.Therefore, we adopted the constant non-LTE correction of +0.05 from Mashonkina (2013) for the Mg abundances of the entire sample.
For titanium, the abundances obtained from line λ4911 Å have been found to be affected by negligible non-LTE effects (Bergemann 2011).Therefore, we used this line as the main indicator of the Ti abundances.For the other lines, small offsets as a function of [Fe/H] were applied to correct all other lines on the scale defined by the λ4911 Å lines; see For silicon, all the lines used here seem to be free from non-LTE effects (Shi et al. 2009(Shi et al. , 2011)).Therefore, the final abundances are given by averaging the line values.For calcium, we extracted the non-LTE corrections (Mashonkina et al. 2007) for lines λ6166 and λ6169 Å by running the NLTE MPIA online tool (Kovalev & Bergemann 2018).For line λ6166 Å, the non-LTE corrections are zero for all combination of atmospheric parameters.For line λ6169 Å, the corrections vary with T eff and [Fe/H].Therefore, we chose the abundances of λ6166 Å as the standard reference to examine the behaviour of the relative abundances of the other three lines.For Eu, the work done in Zhao et al. (2016) showed that the non-LTE effects on abundances computed using the line at λ4129 Å are negligible in metal-poor turn-off stars.Therefore, we used the abundances from this line as a reference and corrected for the relative offsets of abundances from λ4205 Å as function of [Fe/H].Figure A.5 compares the abundances of both lines.
For Ba, the work of Mashonkina et al. (1999) shows that the abundances determined from the line at λ5853 Å are not affected A18, page 4 of 16

Uncertainties of the abundances
The uncertainties of the abundances were computed by estimating the separated effect of each atmospheric parameter.The typical uncertainties of the Titans atmospheric parameters are ±40 K in T eff , ±0.04 dex in log g, ±0.05 dex in [Fe/H], and ±0.1 km s −1 in v mic .New abundances are computed for each spectral line by changing each parameter, one at a time, by its typical error.The influence of errors in v mic was found to be negligible.To estimate the final uncertainty, all changes in the abundances are then added quadratically.Finally, to that uncertainty estimate, we also add the standard deviation of the abundances of the multiple lines.For Mg, for example, this is of ±0.023 dex, as shown in Fig. 3.The typical uncertainties are 0.037, 0.023, 0.029, 0.031, 0.041, 0.080, and 0.041 for A(Mg), A(Si), A(Ca), A(Ti), A(Ni), A(Ba), and A(Eu), respectively.

Magnesium abundances of the GALAH stars
Corrections were applied to the Mg abundances of the GALAH stars to ensure consistency with the abundance scale of the Titans.For that, we used the offset between the new values of T eff and log g that we adopted here and the original values from GALAH DR3.As we mentioned before, the values of [Fe/H] did not change, so there is no offset in this parameter.The offsets in the parameters were mapped into abundance offsets using the same approach adopted to calculate the abundance uncertainties.We found that the effect of the log g offset in Mg abundance is negligible.Using the line λ5711 Å for this analysis, we found that the necessary change in the Mg abundances from GALAH is an offset of +0.07 dex, independent of atmospheric parameters.The fact that the correction is just a constant offset guarantees that the change in [Mg/Fe] values along the age-metallicity relation discussed in Sect.7 (Fig. 7) is not an artefact of our abundance scale correction.The GALAH Mg abundances already take into account non-LTE effects, so this type of correction does not need to be applied again.

Separation of stellar populations
Our goal here is to classify the combined (GALAH + Titans) sample of stars into those stars that have an accreted origin and those that most likely formed in situ.As a first step to separate the metal-poor stellar populations, we used the Toomre diagram to select and exclude stars with prograde motion and total velocity below 210 km s −1 (see Fig. 4).Most stars with total velocities lower than that are probably part of the Milky Way disc.In addition, the stars Wolf 1492 (G64-12), BPS CS 22166-30, and CD-71 1234 (the first two are the most metal-poor stars in the sample) are the only Titans to reach maximum distances from the Galactic plane above 20 kpc.They have been considered to be part of the outer halo population and, as such, they are not used in our main discussion.The remaining Titans make up the sample from which we have the highest probability of identifying stars originating from, in particular, the Gaia-Enceladus merger.
Although there is no certain way to identify stars of accreted origin, there are a few telltale signs that can be used.In particular, previous mergers are thought to be the main source of stars that are usually found in loosely bound orbits that can be retrograde, with high spatial velocities, and/or are shown to be metal-poor, but displaying a low [α/Fe] ratio (Majewski 1992;Gratton et al. 2003).To first order, our classification is based on chemical differences; that is, we consider stars with a lower [Mg/Fe] ratio, at a given metallicity, to have a higher chance of have been accreted during a merger.This separation is then supported and improved by verifying the orbital and kinematic signatures of the accreted populations.In this procedure, the  Titans and GALAH samples were combined and used as a single sample.The classification of the populations was implemented using a Gaussian mixture model clustering algorithm (scikit-learn, Pedregosa et al. 2011) and proceeded as follows.First, the space made of angular momentum (L Z ) and orbital energy (E), the so-called Lindblad diagram, was divided into many over-lapping sections of dimensions of 0.5 [10 3 km s −1 kpc −1 ] × 0.2 [10 5 km s −2 ] in ∆L Z × ∆Energy.Each region was examined separately to find the dominant pattern in the [Mg/Fe] vs. [Fe/H] space.Mixtures of two, three, and four clusters were tested by visual inspection.In general, a configuration with three clusters seemed to better capture the division between two groups with low and high [Mg/Fe] and a third group of stars with uncertain classification.
From the literature (e.g., Helmi et al. 2018;Feuillet et al. 2021), it is known that retrograde stars tend to have higher [Mg/Fe] ratios and are thus separated from Gaia-Enceladus stars in the chemical diagram.Cone-shaped areas in the Lindblad diagram were tested to determine the L Z limits that correspond to such a chemical division.In this exercise, it was found that retrograde stars with high [Mg/Fe] extend to very low [Fe/H] values; this is the population we call Erebus, a name chosen for the reasons explained in Sect.7. On the other hand, stars that have null net rotation and high [Mg/Fe] mostly have high values of [Fe/H]; this is the population that we identify with the Splash.The division between Erebus and Splash is further refined as follows.First, we removed the low [Mg/Fe] stars, identified with Gaia-Enceladus, from the sample.Then, we examined the area that the remaining stars occupy in the Lindblad diagram (Fig. 5), where a clustering of two groups was applied.The top and mid panels of the figure show the classification obtained and, in grey scale, the probability density distribution of the groups.The bottom panel of the figure shows that the distribution in terms of angular momentum is clearly bimodal and well-fit by two Gaussians.The peak at prograde orbits is dominated by Splash stars, while the peak at slightly retrograde orbits, L Z = −0.26× 10 3 km s −1 kpc −1 , is dominated by the stars that we associate with Erebus.This value is likely displaced from zero to slightly negative because of the initial cut of stars in the prograde area where mostly disc stars remain.This shows that the two populations extend beyond the limits set by the clustering algorithm.On the other hand, it is evident that each group is certainly dominated by stars with their own characteristic Galactic rotation pattern.The final selection of the Gaia-Enceladus, Erebus, and Splash stars is shown in red, blue, and grey contours, respectively, in the diagrams of Figs. 4 and 6.While the distribution of dynamic properties of these groups overlaps somewhat, in the next Section we show how the different stellar groups have well-separated relations in age-[Fe/H]-[Mg/Fe] diagrams.Separation in these age-chemical relations shows that our selection criteria successfully group stars of common origin, perhaps with only a small fraction of cross-contamination among them.
Four of the Titans, BD+26 2621 (11.9 Gyr), HD 74000 (9.5 Gyr), Ross 453 (14.6 Gyr), and HD 106038 (13.8 Gyr), are classified as stars from Gaia-Enceladus, but are located in the retrograde area of the Lindblad diagram.These four stars also have the lowest eccentricity values among the Titans (0.6-0.7, see the lower panel in Fig. 4).Although Gaia-Enceladus stars are usually assumed to have high eccentricities (Naidu et al. 2020), it has also been suggested that its stars can have a range of eccentricity values (Perottoni et al. 2022).In our analysis, what classifies these four stars as members of Gaia-Enceladus is their low [Mg/Fe] ratio.Except for HD 74000, they are old stars born within the first two billion years after the Big Bang.HD 74000, although it exhibits a relatively young age among this sample, it still has a low metallicity value ([Fe/H] = −1.85),which places it to the right of the main Gaia-Enceladus sequence in the agemetallicity plot of Fig. 8. Due to its age, it seems to have a low probability of belonging to the population that we identified as Erebus, which is mainly composed of stars located to the A18, page 6 of 16 [Fe/H] plot (Fig. 6).However, its age seems to be in better agreement with the path of the Erebus [Fe/H]-age sequence in Fig. 8. HD 106038 is known to be a chemically peculiar star that likely underwent the influence of a nearby hypernova (Smiljanic et al. 2008).Most likely, its position in the age-metallicity diagram cannot be trusted.In any case, we remark that removing the four stars from the Gaia-Enceladus population does not change our conclusions.The remaining stars in Fig. 8 would still define the same age-metallicity sequence for Gaia-Enceladus.This discussion just highlights that, although we cannot be sure how to assign the population to each individual star, the general properties of the populations that we defined seem to be robust.

Discussion and conclusions
Distinguishing between metal-poor stars that have been formed in accreted systems or in different in-situ populations is by no means straightforward (Jean-Baptiste et al. 2017;Pagnini et al. 2023).Stars considered to have a high probability of having an accreted origin are usually those found in loosely bound orbits that can be retrograde, found to have high spatial velocities, and/or those that are relatively metal-poor, but display a low [α/Fe] ratio (Majewski 1992;Gratton et al. 2003).However, there is no single criterion that can determine the stellar origin with full certainty.To assign our sample stars to a given stellar population or to a given accreted system, we make use of combined information coming from kinematics, orbits, and elemental abundances.We believe to have reached a final classification of high confidence, where even if some individual cases are incorrect, our conclusions are not strongly affected.Taking into account the 45 Titans, with precise ages, we assigned 17 stars to the Milky Way disc, 20 stars to Gaia-Enceladus, three stars to the outer halo, and five stars to an inner halo of probable in situ origin (Erebus, see below).The GALAH sample with precise ages has 183 stars.We assign 61 of them to Gaia-Enceladus, 92 to the A18, page 7 of 16  Figure 7 shows the chemical enrichment sequence of the Gaia-Enceladus galaxy.This clear chronological step-by-step enrichment of Gaia-Enceladus is revealed here for the first time, thanks to the high precision of the stellar ages.It can be seen that some dispersion remains on the age axis.This can be partially ascribed to uncertainties, which are lower than ±0.9 for the Titans and lower than ±1.2 Gyr for GALAH stars.In addition, part of the dispersion is likely intrinsic.It is caused by the stochastic character of chemical evolution in spatially confined groups, which keeps the interstellar medium inhomogeneous at early times (Parizot et al. 2004).
The [Mg/Fe] ratio remains elevated in stars of Gaia-Enceladus up to an age of 11.4 Gyr ([Fe/H] ∼ −1.5 dex).At this point, the ejecta from supernovae (SN) type Ia became important, and the [Mg/Fe] ratio started to decrease.The metallicity at which this so-called 'knee' occurs in Gaia-Enceladus is six to seven times smaller in comparison to the metallicity of Galactic thick disc stars where the same feature is seen ([Fe/H] ∼ −0.7, grey circles in the top panel of Fig. 7).The age we find for the 'knee' differs from that of the Gaia-Enceladus chemical evolution model presented by Vincenzo et al. (2019;∼13.5 Gyr).Our result indicates that a revision of this model is in order.
Figure 7 shows that the chemical evolution of Gaia-Enceladus was truncated 9.6 ± 0.2 Gyr ago.This moment marks the final stage of the Gaia-Enceladus merger.Values for this age between 8 and 10 Gyr have previously been found in the literature (Gallart et al. 2019;Montalbán et al. 2021;Grunblatt et al. 2021;Borre et al. 2022), but here this moment is dated with greater precision.No star in Gaia-Enceladus was found to continue the sequence to metallicities higher than about −0.7 dex; this is similar to the findings of Feuillet et al. (2021).Simulations have shown that after a gas-rich merger, a thin disk can form from gas that has the time to be polluted by type Ia supernovae (Brook et al. 2007).Although we cannot confirm the causal connection here, our dating for the Gaia-Enceladus merger is consistent with the idea that its gas will be used to form the low-α thin disc population (Bonaca et al. 2020), since these stars mostly have younger ages (Haywood et al. 2018).
In Fig. 8, we investigate the time evolution of the chemical enrichment in other metal-poor stellar populations identified in our sample.The plots focus on the so-called Splash (Belokurov et al. 2020) and on a population that we call Erebus.The Splash is made of old, somewhat metal-rich stars (mainly with [Fe/H] > −0.8), probably formed in the thick disc of the Milky Way.Its stars mostly have prograde motions and were thrown into high-eccentricity orbits during the accretion of Gaia-Enceladus (Bonaca et al. 2017;Di Matteo et al. 2019;Belokurov et al. 2020).Erebus, on the other hand, is a population of more metal-poor stars (mostly [Fe/H] < −0.8), mostly in slightly retrograde orbits, although also including a wide tail of stars with prograde orbits (see Fig. 5).
As seen in Fig. 8, the age-metallicity sequence of the Splash stars overlaps that of thick-disc stars.Our precise age dating clearly shows that most of the Splash stars formed before the final stages of the Gaia-Enceladus merger.The peak of their age distribution is around 11 Gyr, and the mean metallicity is around [Fe/H] ∼ −0.7.After the peak, the age distribution shows a decline for about 2 Gyr, with its minimum roughly coinciding with the age truncation of Gaia-Enceladus.Subsequently, a slight increase is observed.However, the significance of this increase is unclear.Splash stars were only found in the GALAH sample, for which we obtained ages with slightly higher mean uncertainty.
Regarding Erebus, and as several halo substructures have already been identified and named in the literature (e.g., Koppelman et al. 2019;Myeong et al. 2019;Naidu et al. 2020;Horta et al. 2021;Lövdal et al. 2022), a valid question is whether we have rediscovered a substructure that is already known.In the region occupied by Erebus in the mid-panel of Fig. 5, there are two main substructures identified in the literature, Sequoia (Myeong et al. 2019) and Thamnos (Koppelman et al. 2019) The size of the points change according to stellar age and the color according to the [Mg/Fe] ratio.The Gaia-Enceladus sequence is shown as the red line.For the Splash, only stars with ages that have precision better than ±1.5 Gyr are considered.For Erebus, we plot those with age precision below ±2 Gyr.Thick and thin disc stars from Haywood et al. (2018) are shown as green and black crosses, respectively.Bottom: [Mg/Fe] ratio as a function of metallicity for stars in the Splash (left) and Erebus (right) components.Points are color-and size-coded according to the stellar age.
diagram.In this region, we only have a handful of objects.Thamnos stars, as defined by Koppelman et al. (2019), would be made of two components with, among other properties, V φ centred at −150 and −200 km s −1 .Instead, what we call Erebus appears strongly concentrated at V φ > −100 km s −1 (see Fig. 4), although we deduce by the analysis in Fig. 5 that part of it remains masked by the Splash rotation distribution, composing a dynamic structure with virtually null net rotation L Z ≈ 0. We cannot exclude the possibility that a few Sequoia and/or Thamnos stars are part of our Erebus selection, but it seems clear that Erebus is not dominated by stars belonging to these substructures.
An important point to note here is that our sample is very small (228 stars with precise ages) compared to large samples, with thousands or tens of thousands of stars, used by works that search for halo substructures.Therefore, we are not in a position to identify overdensities of stars that distinguish themselves from a general background stellar population.In contrast, what we call Erebus is, most likely, exactly such background stellar population.As we discuss below, the age-[Fe/H]-[α/Fe] relations of these stars suggest that they are most likely of in situ origin and a significant fraction of them acquired retrograde orbits after being subject to the effects of several mergers in the history of the Milky Way (Jean-Baptiste et al. 2017).
Figure 8 shows that the stars of Erebus attained a higher metallicity earlier than the stars from Gaia-Enceladus.This fact suggests that the Erebus stars formed in a system or region that could sustain star formation with higher efficiency.Moreover, Erebus stars consistently show a high [Mg/Fe] up to [Fe/H] ∼ −0.5.Only around this metallicity does the SNIa contribution becomes apparent.Together, these facts suggest that Erebus originated in a system that was more massive than Gaia-Enceladus.
The most obvious candidate is the Milky Way itself, where the knee in [Mg/Fe] ratios also appears at [Fe/H] ∼ −0.5.
The age distribution of Erebus stars peaks at a value of ∼11.5 Gyr, slightly older but still similar to the Splash.Interestingly, essentially no star in this group was found to be younger than the age of the final stage of the Gaia-Enceladus merger.Because the Erebus stars are old and mostly retrograde, it is possible that their orbits show the combined effect of heating caused by several mergers.As discussed in Jean-Baptiste et al. ( 2017), the higher the number of mergers suffered by the host galaxy, the broader the area in the retrograde region of the Lindblad diagram (see e.g., Fig. 5) that the in situ stars will end up occupying.Alternatively, there are at least two other possibilities for the origin of the Erebus stars: (i) the core region of Gaia-Enceladus itself and (ii) a separated merger with a galaxy more massive than Gaia-Enceladus.
The simulations of Koppelman et al. (2020) demonstrated that debris from a merger with a disc galaxy (as was probably the Gaia-Enceladus progenitor, Helmi et al. 2018) have a final complex phase-space structure, a range of orbital properties, and also stars with a range of chemical abundances.In the core regions of disc galaxies, the star formation is more intense.Consequently, as seen in the Milky Way, stars formed in the inner regions can be old, metal rich, and have enhanced [Mg/Fe] (Trevisan et al. 2011).In that sense, what we call Erebus could be a population of stars from the core of Gaia-Enceladus.The decreasing age distribution of Erebus, after 11.5 Gyr, and potential cut-off around 9.5 Gyr, could be understood as being influenced by the on-going merger.However, the simulations of Koppelman et al. (2020) also show that for the assumed configuration of the Gaia-Enceladus merger, stars from the core region would end up with prograde orbits of low eccentricity, which is not the case for the A18, page 9 of 16 Erebus stars.Thus, it seems that this hypothesis for the origin of Erebus can be excluded.
Another possibility is that of a separated merger with a system that was more massive than the Gaia-Enceladus galaxy.However, as we discussed above, the kinematic and orbital properties of the Erebus stars do not seem to match the substructures identified in the same region of the Linblad diagram in other works (Koppelman et al. 2019;Naidu et al. 2020;Shank et al. 2022).On the other hand, one should also keep in mind that a single merger can give origin to several substructures (Jean-Baptiste et al. 2017).To investigate any connection between Erebus and other halo substructures, we will need to extend our precise analysis of abundances and ages to stars from these other substructures for comparison.What we do observe is that the age-metallicity relation for Erebus stars is distinct from that of Gaia-Enceladus.At this point, a connection between the two populations seems unlikely.
Here, we highlight a possible connection with the simulations of galaxy formation in a cosmological context considered in Belokurov et al. (2020).As those authors discussed, some of the Auriga simulations where a host galaxy had a proto-disc produced a stellar population with properties consistent with those of the Splash stars after a merger (where the heated stars managed to retain much of their angular momentum).However, in one of the simulations, the signature of a prograde Splash was less clear.The authors speculated that the original configuration of the heated stars was less of a disc and more spheroidal.This could be a possible explanation for the Erebus stars.This population could be a component that was already in a spheroid configuration but that suffered additional heating by the Gaia-Enceladus accretion and perhaps also by other previous merger events.This possibility lends additional support to our conclusion that Erebus is a stellar population of in-situ origin.Because of the complex nature of this population and its somewhat "chaotic" orbital characteristics, we got inspired to name it Erebus, a primordial deity in Greek mythology born out of Chaos.
A similar in situ stellar component, dubbed Aurora, has been recently identified by Belokurov & Kravtsov (2022).Aurora was defined with a metallicity cut of [Fe/H] > −1.5 and a strict selection of energy values (on a scale quite different from the one used here, making a comparison difficult).Aurora stars were found to have a modest net rotation.This is different from our selection, which includes more metal-poor stars without restriction in energy, although it also appears to be of null net rotation.We cannot exclude that these are two sides of the same in-situ population, with Erebus being a more general and extended population that encompasses Aurora.We also note that a suggestion has been made by Myeong et al. (2022) that Aurora is connected to the Heracles structure identified by Horta et al. (2021) in the inner 4 kpc of the Galaxy (which is a region not probed by our stellar sample).On the other hand, Horta et al. (2023) finds clear chemical differences between stars associated with Aurora and Heracles.As a connection between Erebus and all these populations is, at least for now, not clear, we believe that giving it a separated identification is justified.Precise ages for stars in all these groups will be needed to clarify whether their age-chemical relations are distinct or not.
Finally, we mention that the plots in Fig. 8 include stars with orbits and [α/Fe] ratios consistent with those in the thin disc but with atypically old ages (Haywood et al. 2013(Haywood et al. , 2019)).Judging by their values of age, metallicity, and [α/Fe], these stars seem to be part of a coeval population of high prominence in the inner disc (up to 10 kpc radius) identified by Ciucȃ et al. (2023).
On the basis of comparisons with the cosmological simulations from Grand et al. (2020), Ciucȃ et al. (2023) interpret this population as evidence of a starburst generated by gas accretion induced by the merging with Gaia-Enceladus.This new generation of stars then has ages roughly consistent with those of the Splash and a range of [Fe/H] values that are lower than typical values of thin-disc stars due to dilution with the accreted metal-poor gas.
The discussion above shows how much the precise ages and abundances derived in this work help to put together all the pieces of the old Milky Way populations in a consistent chronological order.We can clearly detect the halo in-situ stellar population that existed before the Gaia-Enceladus merger and was heated to more retrograde orbits.We can also detect the stellar population of the disc that was affected by the same merger.Finally, we can also precisely date the final stages of the Gaia-Enceladus merger, at 9.6 ± 0.2 Gyr.This dating is also consistent with the idea that the gas from Gaia-Enceladus helps to form the low-α population of the thin disc (Brook et al. 2007;Bonaca et al. 2020).The data produced in this paper used for the discussion in this section are provided as supplementary material at CDS.

Fig. 1 .
Fig. 1.Kiel diagram with the stars divided into the different stellar populations.The symbols are color-coded according to metallicity as indicated in the scale on the right side.The symbol size is proportional to the stellar age.Stars belonging to Gaia-Enceladus are highlighted by dark contours and have age precision better than 1.2 Gyr.Stars belonging to the Splash and Erebus have age precision better than 1.5 Gyr and 2.0 Gyr, respectively.Isochrones are overplotted for reference.The isochrone corresponding to the youngest Gaia-Enceladus population discussed in Sect.7 is highlighted in red.These stars are the same shown in Figs.7 and 8.

Fig. 2 .
Fig. 2. Offset between spectroscopic and IRFM T eff .The distribution includes 408 turnoff stars from the GALAH catalog, restricted according to the selection explained in Sect. 2. Only stars with halo-like orbits, i.e., those confined within the retrograde and null net rotation areas in Fig. 4, are shown.The thick continuous line connects the medians computed in equally spaced bins of IRFM T eff .Vertical bars indicate the 1σ dispersion.
Fig. A.1.The final abundance of each star was then computed by averaging the values.
Figure A.2 shows a comparison among the abundances of each line.For Ni, on the other hand, we did not find studies on non-LTE and 3D effects for the lines we used in this work.No correction could be applied in this case.The final Ni abundances are the averages of the line values.Figure A.3 compares the abundances of both lines used here.
Figure A.4 shows a comparison of the abundances obtained with these lines.Similarly to what was done above for Mg and Ti, the final Ca abundances were obtained by averaging the offsetcorrected values from each line.

Fig. 3 .
Fig. 3. Consistency of Mg abundance from the lines 5528 and 5711 Å.For the Titans, difference between Mg abundances from the lines λ5528 Å and λ5711 Å as functions of the atmospheric parameters.The red lines are the linear regressions and the shaded areas indicate 1σ and 2σ dispersions.The reference stars HD 22879, HD 140283, and HD 84937 are identified with bold contours.
Figure A.7  shows plots with the abundances of the elements above relative to hydrogen as a function of stellar age, including only the Titans.Only the Enceladus and Erebus populations (defined in Sect.6) are distinguished because stars from other populations were not found in this sample.The main discussion in Sect.7 uses mostly the abundances of magnesium.

Fig. 4 .
Fig. 4. Inner halo populations in kinematic and dynamical diagrams.Top panel: Lindblad diagram.The contours represent the 33, 66, 90, 95, 98, 99, and 99.9% cumulative distribution of the populations disentangled with the Titans and GALAH metal-poor stars.The regions and stars corresponding to the Gaia-Enceladus, Splash, and Erebus populations are highlighted.Here, the Titans that have disc dynamics are also shown for completeness.Mid panel: Toomre diagram.Symbols and colors are the same as in the top panel.The region of disc stars, selected following Helmi et al. (2018), is delimited by the dashed line.Bottom panel: eccentricity vs. maximum distance from the Galactic plane (Z max ).Symbols and colors are the same as in the other panels.

Fig. 5 .
Fig. 5. Lindblad diagram of the Splash and Erebus populations.Top and middle panels highlight the probability density of the Splash and Erebus stars, respectively.Bottom panel shows the histogram of angular momentum, L Z , for stars with binding energy between −1.7 and −1.9 × 10 5 km −2 , approximately where most stars are concentrated.The histogram is well fit by a bi-modal Gaussian distribution (in black).One peak is dominated by Splash stars (solid bright red line) and the other by Erebus (dashed dark red line).

Fig. 6 .
Fig. 6.Inner halo populations in the [Mg/Fe] vs. [Fe/H] diagram.Cumulative distribution plots of the populations selected with GALAH in the [Mg/Fe] vs. [Fe/H] space.Contour colors are consistent with those used in Fig. 4. The symbols display the Titans associated to Gaia-Enceladus and Erebus, as well as two unclassified stars.Top panel: regions occupied by the stars from Gaia-Enceladus and the Splash.Bottom panel: regions occupied by the stars from Gaia-Enceladus and Erebus.

Fig. 7 .
Fig. 7.Chemical evolution of Gaia-Enceladus.Top panel: time evolution of the chemical enrichment of stars from Gaia-Enceladus.The symbol size changes as a function of stellar age and the points are colorcoded according to [Mg/Fe].Only stars for which an age precision better than ±1.2 Gyr was obtained are shown.The gray line connects median ages computed in equally separated bins of [Fe/H].The bars correspond to the 25 and 75% quartiles around the medians.The Titans are distinguished from the GALAH stars by the dark contours.Thick disc stars from Haywood et al. (2018) are shown as gray circles.Bottom panel: [Mg/Fe] vs. [Fe/H] diagram.Circles are color-and size-coded by stellar age.The approximate metallicity value where Fe coming from supernovae type Ia becomes important is indicated at [Fe/H] = −1.5.
Splash, and 30 to Erebus.These stars classified into different populations are shown in kinematic, dynamical, and chemical diagrams in Figs. 4 and 6.

Fig. 8 .
Fig. 8.Chemical evolution of the Splash and Erebus.Top: age-metallicity relation for stars in the Splash (left) and Erebus (right) components.The size of the points change according to stellar age and the color according to the [Mg/Fe] ratio.The Gaia-Enceladus sequence is shown as the red line.For the Splash, only stars with ages that have precision better than ±1.5 Gyr are considered.For Erebus, we plot those with age precision below ±2 Gyr.Thick and thin disc stars fromHaywood et al. (2018) are shown as green and black crosses, respectively.Bottom: [Mg/Fe] ratio as a function of metallicity for stars in the Splash (left) and Erebus (right) components.Points are color-and size-coded according to the stellar age.